update 5

autoupdate

@@ -10,7 +10,7 @@ if [ ! -d .git ]; then
     [ -e version ] && rm version
     git add .gitignore '*' '**/*'
 else
-    git add -u
+    git add -A
     changes="$(git status --porcelain | wc -l)"
     [ "$changes" -eq 0 ] && exit 0
 fi

lib/__init__.py

@@ -60,13 +60,16 @@ def c_render(cascade, precision=4096):
         ys += f(xs)
     return xs, ys
 
-def c_render2(xs, cascade):
+def c_render2(xs, cascade, phase=False):
     """c_render optimized and specifically for first/second-order filters"""
     import numexpr as ne
     j = np.complex(0, 1)
-    eq2 = 'abs((b0 + j*b1*ws - b2*ws**2)/(a0 + j*a1*ws - a2*ws**2))**2'
-    eq1 = 'abs((b0 + j*b1*ws)/(a0 + j*a1*ws))**2'
-    fmt = 'real(log10({})*10 + gain)'
+    eq2 = '(b0 + j*b1*ws - b2*ws**2)/(a0 + j*a1*ws - a2*ws**2)'
+    eq1 = '(b0 + j*b1*ws)/(a0 + j*a1*ws)'
+    if not phase:
+        fmt = 'real(log10(abs({})**2)*10 + gain)'
+    else:
+        fmt = 'arctan2(imag({0}), real({0}))' # gross
     ys = np.zeros(len(xs))
     for f in cascade:
         w0, ba, gain = f
@@ -83,6 +86,8 @@ def c_render2(xs, cascade):
             raise Exception("incompatible cascade; consider using c_render instead")
         ws = xs/w0
         ys += ne.evaluate(eq)
+    if phase:
+        ys = degrees_clamped(ys)
     return ys
 
 def firize(xs, ys, n=4096, srate=44100, plot=None):

lib/plot.py

@@ -1,6 +1,8 @@
 import matplotlib.pyplot as plt
 from matplotlib import ticker
 
+# TODO: remove set_size_inches calls, move them inline as necessary
+
 def response_setup(ax, ymin=-24, ymax=24, yL=ticker.AutoMinorLocator(3)):
     ax.set_xlim(20, 20000)
     ax.set_ylim(ymin, ymax)
@@ -10,6 +12,15 @@ def response_setup(ax, ymin=-24, ymax=24, yL=ticker.AutoMinorLocator(3)):
     ax.set_xlabel('frequency (Hz)')
     ax.set_ylabel('magnitude (dB)')
 
+def phase_response_setup(ax, div=12, yL=ticker.AutoMinorLocator(2)):
+    ax.set_xlim(20, 20000)
+    ax.set_ylim(-180, 180)
+    ax.set_yticks(tuple(range(-180, 180 + 1, int(360/div))))
+    ax.yaxis.set_minor_locator(yL)
+    ax.grid(True, 'both')
+    ax.set_xlabel('frequency (Hz)')
+    ax.set_ylabel('phase (degrees)')
+
 def cleanplot():
     fig, ax = plt.subplots()
     fig.set_size_inches((16, 16))
@@ -18,9 +29,29 @@ def cleanplot():
     return fig, ax
 
 def new_response(*args, **kwargs):
-    #fig, ax = plt.subplots()
     fig = plt.figure()
     ax = fig.gca()
     response_setup(ax, *args, **kwargs)
     fig.set_size_inches(10, 6)
     return fig, ax
+
+def new_phase_response(*args, **kwargs):
+    fig = plt.figure()
+    ax = fig.gca()
+    phase_response_setup(ax, *args, **kwargs)
+    fig.set_size_inches(10, 6)
+    return fig, ax
+
+def new_bode(magnitude_offset=0):
+    fig = plt.figure()
+    ax1 = fig.gca()
+    ax2 = ax1.twinx()
+    ymin = -24 + magnitude_offset
+    ymax = 24 + magnitude_offset
+    response_setup(ax1, ymin, ymax)
+    phase_response_setup(ax2)
+    # ax1 and ax2 should have identical grids;
+    # disable ax2's so it doesn't overlap ax1's plot lines.
+    ax2.grid(False, which='both')
+    fig.set_size_inches(10, 6)
+    return fig, ax1, ax2

lib/plotwav.py

@@ -0,0 +1,126 @@
+# this is a bunch of crap that should really be reduced to one or two functions
+
+from . import wav_read, normalize, averfft, tilter2, smoothfft2, firize
+from . import new_response, show, convolve_each, monoize, count_channels
+
+import numpy as np
+
+def plotwavsmooth(fn, ax, tilt=None, bw=1, size=8192, raw=False, fix=False, smoother=smoothfft2, **kwargs):
+    s, srate = wav_read(fn)
+
+    s, rms = normalize(s, srate)
+    sm = monoize(s)
+    ss = monoize(s*np.array((1, -1)))
+
+    xs_raw = np.arange(0, srate/2, srate/2/size)
+    ys_raw = averfft(sm, size=size)
+
+    # tilting beforehand is negligible besides lowest frequencies, but eh
+    if tilt is not None:
+        ys_raw -= tilter2(xs_raw, tilt)
+
+    xs, ys = smoother(xs_raw, ys_raw, bw=bw)
+
+    if not 'label' in kwargs:
+        kwargs['label'] = fn
+
+    if raw:
+        ax.semilogx(xs_raw, ys_raw, **kwargs)
+    ax.semilogx(xs, ys, **kwargs)
+
+    if not fix: return
+
+    fno = fn[:-4]+"-proc.wav"
+    fir = firize(xs, -ys, srate=srate)
+    sf = convolve_each(s/8, fir, mode='same')
+
+    import ewave
+    with ewave.open(fno, 'w', sampling_rate=srate, nchannels=count_channels(sf)) as f:
+        f.write(sf)
+    print('wrote '+fno)
+
+def plotfftsmooth(s, srate, ax, bw=1, tilt=None, size=8192, window=0, raw=False, **kwargs):
+    sm = monoize(s)
+
+    xs_raw = np.arange(0, srate/2, srate/2/size)
+    ys_raw = averfft(sm, size=size, mode=window)
+
+    ys_raw -= tilter2(xs_raw, tilt)
+
+    xs, ys = smoothfft(xs_raw, ys_raw, bw=bw)
+
+    if raw: ax.semilogx(xs_raw, ys_raw, **kwargs)
+    ax.semilogx(xs, ys, **kwargs)
+
+    return xs, ys
+
+def plotwav2(fn, ax, bw=1, size=8192, raw=False, fix=False,
+             smoother=smoothfft2, side_compensate=9, **kwargs):
+    s, srate = wav_read(fn)
+
+    s, rms = normalize(s, srate)
+    sm = monoize(s)
+    ss = monoize(s*np.array((1, -1)))
+
+    xs_raw = np.arange(0, srate/2, srate/2/size)
+    ys_raw = averfft(sm, size=size)
+    ys_raw_side = averfft(ss, size=size)
+
+    # tilting beforehand is negligible besides lowest frequencies, but eh
+    ys_raw -= tilter2(xs_raw, 'np2')
+    ys_raw_side -= tilter2(xs_raw, 'np2s')
+
+    xs, ys = smoother(xs_raw, ys_raw, bw=bw)
+    xs, ys_side = smoother(xs_raw, ys_raw_side, bw=bw)
+
+    if not 'label' in kwargs:
+        kwargs['label'] = fn
+
+    if raw:
+        ax.semilogx(xs_raw, ys_raw, **kwargs)
+        ax.semilogx(xs_raw, ys_raw_side + side_compensate, **kwargs)
+    ax.semilogx(xs, ys, **kwargs)
+    ax.semilogx(xs, ys_side + side_compensate, **kwargs)
+
+    side_gain = np.average(ys_raw_side) - np.average(ys_raw)
+    print("side gain:", side_gain)
+
+    if not fix: return
+    fno = fn[:-4]+"-proc.wav"
+
+    fir = firize(xs, -ys, srate=srate)
+    smf = convolve_each(sm/8, fir, mode='same')
+    fir = firize(xs, -ys_side, srate=srate)
+    ssf = convolve_each(ss/8, fir, mode='same')
+    ssf *= 10**(side_gain/20)
+    sf = np.array((smf + ssf, smf - ssf)).T
+
+    import ewave
+    with ewave.open(fno, 'w', sampling_rate=srate, nchannels=count_channels(sf)) as f:
+        f.write(sf)
+    print('wrote '+fno)
+
+def pw(fn, ax, **kwargs):
+    plotwavsmooth(fn, ax, tilt='np2', bw=1/6, **kwargs)
+
+def pwc(fn, **kwargs):
+    fig, ax = new_response(-18, 18)
+    ax.set_title('averaged magnitudes of normalized songs with tilt and smoothing')
+
+    pw(fn, ax, fix=True, **kwargs)
+    fno = fn[:-4]+"-proc.wav"
+    pw(fno, ax, fix=False, **kwargs)
+
+    ax.legend(loc=8)
+    show(fig)
+
+def pw2(fn, **kwargs):
+    fig, ax = new_response(-18, 18)
+    ax.set_title('averaged magnitudes of normalized songs with tilt and smoothing')
+
+    plotwav2(fn, ax, fix=True, bw=1/6, **kwargs)
+    fno = fn[:-4]+"-proc.wav"
+    plotwav2(fno, ax, fix=False, bw=1/6, **kwargs)
+
+    ax.legend(loc=8)
+    show(fig)

lib/util.py

@@ -17,6 +17,8 @@ warp = lambda w: np.arctan(w)*2
 rfft = lambda src, size: np.fft.rfft(src, size*2)
 magnitude = lambda src, size: 10*np.log10(np.abs(rfft(src, size))**2)[0:size]
 
+degrees_clamped = lambda x: ((x*180/np.pi + 180) % 360) - 180
+
 def xsp(precision=4096):
     """create #precision log-spaced points from 20 to 20480 Hz"""
     # i opt not to use steps or linspace here,