update 1

.gitignore

@@ -0,0 +1,4 @@
+*.pyc
+__pycache__/*
+lib/__main__.py
+version

README.md

@@ -0,0 +1,51 @@
+# random dsp code
+
+it's a bunch of half-baked python code that's kinda handy.
+
+don't expect commits, docs, or comments to be any verbose.
+
+## the stuff
+
+* a basic BS.1770-3 normalization implementation
+— [bs.py](/lib/bs.py)
+
+* biquad butterworth/chebyshev filters [(via DSPFilters)][dspf]
+— [butterworth.py](/lib/butterworth.py)
+
+* s-plane to z-plane conversion
+— [planes.py](/lib/planes.py)
+
+* various functions for biquad filters
+— [bq.py](/lib/bq.py) [\_\_init\_\_.py](/lib/__init__.py)
+
+* some functions for state-variable filters [(via Raph Levien)][svf]
+— [svf.py](/lib/svf.py)
+
+* sine sweeps, and the Optimized Aoshima's Time-Stretched Pulse [(via here)][sweeps]
+— [sweeps.py](/lib/sweeps.py)
+
+* a lot of stuff for magnitude plotting like tilting and smoothing
+— [fft.py](/lib/fft.py) [smoothfft.py](/lib/smoothfft.py) [\_\_init\_\_.py](/lib/__init__.py)
+
+* rough experiments with psychoacoustic equalization ("neon pink" and other crap)
+— [data.py](/lib/data.py) [\_\_init\_\_.py](/lib/__init__.py)
+
+* miscellaneous matplotlib stuff
+— [plot.py](/lib/plot.py)
+
+* miscellaneous utility functions
+— [util.py](/lib/util.py) [wav.py](/lib/wav.py)
+
+[dspf]: https://github.com/vinniefalco/DSPFilters/
+[sweeps]: http://www.sound.sie.dendai.ac.jp/dsp/e-21.html
+[svf]: http://nbviewer.ipython.org/urls/music-synthesizer-for-android.googlecode.com/git/lab/Second%20order%20sections%20in%20matrix%20form.ipynb
+
+all wrapped up in a inconveniently generic "lib" module!
+
+## dependencies
+
+python 3.4+
+
+numpy scipy sympy matplotlib ewave
+
+usually run in an ipython environment.

autoupdate

@@ -0,0 +1,25 @@
+#!/usr/bin/env zsh
+set -e
+alias db="dropbox_uploader"
+
+rm -r lib
+db download py/lib >/dev/null
+
+if [ ! -d .git ]; then
+    git init >/dev/null
+    [ -e version ] && rm version
+    git add .gitignore '*' '**/*'
+else
+    git add -u || echo "nothing new?"
+fi
+
+if [ -s version ]; then
+    read -r version < version
+else
+    version=0
+fi
+
+let version++ || true # -e mode doesn't like this normally
+echo $version > version
+
+git commit -m "update $version" . >/dev/null

lib/__init__.py

@@ -0,0 +1,247 @@
+import numpy as np
+#from IPython.display import display
+from matplotlib.pylab import show
+
+from .util import *
+from .data import *
+gen_filters = lambda cascade, srate: [
+    s2z(*f[1], fc=f[0], srate=srate, gain=10**(f[2]/20)) for f in cascade
+]
+from .bq import *
+from .butterworth import *
+from .sweeps import *
+from .smoothfft import *
+from .plot import *
+from .wav import *
+from .planes import *
+from .fft import *
+from .bs import *
+
+def analog(b, a):
+    import sympy as sym
+    w,s = sym.symbols('w s')
+    filt_expr = sym.Poly(b, s)/sym.Poly(a, s)
+    mag_expr = abs(filt_expr.subs({s: w*sym.I}))**2
+    return sym.lambdify(w, mag_expr, 'numpy')
+
+def makemag(w0, ba, gain=0):
+    f = analog(*ba)
+    def magf(w):
+        a = f(w/w0)
+        a[0] = 1e-35
+        a = np.log10(a)*10 + gain
+        a[0] = a[1] # safety measure
+        return a
+    return magf
+
+def test_filter_raw(ba, fc=1000, gain=0, precision=4096):
+    fig, ax = new_response(ymin=-24, ymax=24)
+    xs = xsp(precision)
+    ax.semilogx(xs, makemag(fc, ba, gain)(xs))
+    show(fig)
+
+def test_filter(ff, A=toA(12), Q=toQ(1), **kwargs):
+    test_filter_raw(ff(A, Q), **kwargs)
+
+npc = [makemag(*f) for f in cascades['raw']]
+def neonpink(xs):
+    print("neonpink(): DEPRECATED")
+    combined = np.zeros(len(xs))
+    for f in npc:
+        combined += f(xs)
+    return combined
+
+def c_render(cascade, precision=4096):
+    xs = xsp(precision)
+    ys = np.zeros_like(xs)
+    c = [makemag(*f) for f in cascade]
+    for f in c:
+        ys += f(xs)
+    return xs, ys
+
+def firize(xs, ys, n=4096, srate=44100, plot=None):
+    import scipy.signal as sig
+    if plot:
+        plot.semilogx(xs, ys, label='desired')
+    xf = xs/srate*2
+    yg = 10**(ys/20)
+
+    xf = np.hstack((0, xf, 1))
+    yg = np.hstack((0, yg, yg[-1]))
+
+    b = sig.firwin2(n, xf, yg, antisymmetric=True)
+
+    if plot:
+        _, ys = sig.freqz(b, worN=xs/srate*tau)
+        ys = 20*np.log10(np.abs(ys))
+        plot.semilogx(xs, ys, label='FIR ({} taps)'.format(n))
+        plot.legend(loc=8)
+
+    return b
+
+def normalize_test(s, srate):
+    # FIXME: where is fir3 defined?
+    #rms_naive = np.sqrt(np.mean(s**2))
+    filters3 = gen_filters(cascades['raw3'], srate)
+    db_standard = BS1770_3(s, srate)
+    db = BS1770_3(s, srate, filters=filters3)
+    print('raw3 would be\t{:+6.2f} dB louder/quieter than RG2'.format(db_standard - db))
+    db = BS1770_3(s, srate, filters=[fir3])
+    print('fir3 would be\t{:+6.2f} dB louder/quieter than RG2'.format(db_standard - db))
+    rms = 10**(db/20)
+    return s/rms, rms
+
+def tilter(xs, ys, tilt):
+    """tilts a magnitude plot by some decibels, or by equalizer curve."""
+    # should really just do this instead:
+    # ys -= tilt(xs, 3)
+    print("tilter(): DEPRECATED")
+    if tilt == 'neon':
+        noise = neonpink(xs)
+    elif type(tilt) is str:
+        noise = np.zeros(len(xs))
+        c = [makemag(*f) for f in cascades[tilt]]
+        for f in c:
+            noise += f(xs)
+    elif isinstance(tilt, int) or isinstance(tilt, float):
+        noise = tilt*(np.log2(1000) - np.log2(xs))
+    else:
+        noise = np.zeros(xs.shape)
+    return xs, ys - noise
+
+def tilter2(xs, tilt):
+    if type(tilt) is str:
+        noise = np.zeros(len(xs))
+        c = [makemag(*f) for f in cascades[tilt]]
+        for f in c:
+            noise += f(xs)
+    elif isinstance(tilt, int) or isinstance(tilt, float):
+        noise = tilt*(np.log2(1000) - np.log2(xs + 1e-35))
+    else:
+        noise = np.zeros(xs.shape)
+    return noise
+
+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)
+
+# this is similar to default behaviour of having no __all__ variable at all,
+# but ours ignores modules as well. this allows for `import sys` and such
+# without clobbering `from our_module import *`.
+__all__ = [o for o in locals() if type(o) != 'module' and not o.startswith('_')]

lib/bq.py

@@ -0,0 +1,35 @@
+import numpy as np
+import scipy.signal as sig
+
+from .util import *
+
+bq_run = lambda bq, xs: sig.lfilter(*bq, x=xs, axis=0)
+
+nf = lambda t, f, g, bw, mg: (f, t(toA(g), toQ(bw)), mg)
+
+LP1 = lambda A, Q: ((0,1),(1,1))
+HP1 = lambda A, Q: ((1,0),(1,1))
+LS1 = lambda A, Q: ((1,A),(1,1/A))
+HS1 = lambda A, Q: ((A,1),(1/A,1))
+
+# a always gets divided by A instead of multiplied
+# b1 and a1 always /= Q
+
+LP2 = lambda A, Q: ((0,     0, 1),
+                    (1,   1/Q, 1))
+HP2 = lambda A, Q: ((1,     0, 0),
+                    (1,   1/Q, 1))
+PE2 = lambda A, Q: ((1,   A/Q, 1),
+                    (1, 1/A/Q, 1))
+AP2 = lambda A, Q: ((1,   1/Q, 1),
+                    (1,   1/Q, 1))
+BP2a= lambda A, Q: ((0,  -A/Q, 0),
+                    (1, 1/A/Q, 1))
+BP2b= lambda A, Q: ((0,-A*A/Q, 0),
+                    (1,   1/Q, 1))
+NO2 = lambda A, Q: ((1,     0, 1),
+                    (1,   1/Q, 1))
+LS2 = lambda A, Q: ((1,     np.sqrt(A)/Q,   A),
+                    (1,   1/np.sqrt(A)/Q, 1/A))
+HS2 = lambda A, Q: ((A,     np.sqrt(A)/Q,   1),
+                    (1/A, 1/np.sqrt(A)/Q,   1))

lib/bs.py

@@ -0,0 +1,83 @@
+from . import blocks, convolve_each, gen_filters, cascades, bq_run, toLK
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+def BS1770_3(s, srate, filters=None, window=0.4, overlap=0.75, gate=10, absolute_gate=70, detail=False):
+    if filters is None:
+        filters = gen_filters(cascades['1770'], srate)
+
+    sf = np.copy(s)
+    for f in filters:
+        if len(f) is 2: # dumb but effective
+            sf = bq_run(f, sf)
+        else:
+            sf = convolve_each(sf, f, 'same')
+
+    stepsize = round(window*srate*(1 - overlap))
+    blocksize = int(stepsize/(1 - overlap))
+
+    means = np.array([
+        np.sum(np.mean(b**2, axis=0)) for b in blocks(sf, stepsize, blocksize)
+    ])
+    LKs = toLK(means)
+
+    truths = LKs > -absolute_gate
+    LKs_g70 = LKs[truths]
+    means_g70 = means[truths]
+    avg_g70 = np.average(means_g70)
+    threshold = toLK(avg_g70) - gate
+    means_g10 = means[LKs_g70 > threshold]
+    avg_g10 = np.average(means_g10)
+
+    if detail is False:
+        return toLK(avg_g10)
+    else:
+        return toLK(avg_g10), toLK(avg_g70), LKs, threshold
+
+def BS_plot(ys, g10=None, g70=None, threshold=None, fig=None, ax=None):
+    if g10:
+        center = np.round(g10)
+        bins = np.arange(center - 10, center + 10.01, 0.25)
+    else:
+        bins = np.arange(-70, 0.1, 1)
+
+    if fig is None:
+        fig = plt.figure()
+    if ax is None:
+        ax = fig.gca()
+
+    if False: # histogram
+        ax.hist(ys, bins=bins, normed=True, facecolor='g', alpha=0.5)
+        ax.xlim(bins[0], bins[-1])
+        ax.ylim(0, 1)
+        ax.grid(True, 'both')
+        ax.xlabel('loudness (LKFS)')
+        ax.ylabel('probability')
+        fig.set_size_inches(10,4)
+        show()
+
+    xs = np.arange(len(ys))
+    #ax.plot(xs, ys, color='#066ACF', linestyle=':', marker='d', markersize=2)
+    ax.plot(xs, ys, color='#1459E0')
+    ax.set_xlim(xs[0], xs[-1])
+    ax.set_ylim(-70, 0)
+    ax.grid(True, 'both', 'y')
+    ax.set_xlabel('bin')
+    ax.set_ylabel('loudness (LKFS)')
+    fig.set_size_inches(12,5)
+    #_, _, ymin, _ = ax.axis()
+    if threshold:
+        ax.axhspan(-70, threshold, facecolor='r', alpha=1/5)
+    if g10:
+        ax.axhline(g10, color='g')
+    if g70:
+        ax.axhline(g70, color='0.3')
+
+    return fig, ax
+
+def normalize(s, srate):
+    """performs BS.1770-3 normalization and returns inverted gain."""
+    db = BS1770_3(s, srate)
+    rms = 10**(db/20)
+    return s/rms, rms

lib/butterworth.py

@@ -0,0 +1,87 @@
+import numpy as np
+
+def LPB(n):
+    # crap ripped from https://github.com/vinniefalco/DSPFilters/blob/master/shared/DSPFilters/source
+    """n-th degree butterworth low-pass filter cascade
+
+    -3 dB at center frequency."""
+    series = []
+    pi2 = np.pi/2
+    pairs = int(n/2)
+    for i in range(pairs):
+        # magnitude = 1
+        theta = pi2*(1 + (2*i + 1)/n)
+        real = np.cos(theta)
+        imag = np.sin(theta)
+        num = (0, 0, 1)
+        den = (1, -2*real, real*real + imag*imag)
+        series += [(num, den)]
+    if n % 2:
+        num = (0, 1)
+        den = (1, 1)
+        series += [(num, den)]
+    return series
+
+def LPC(n, ripple, type=1):
+    # crap ripped from https://github.com/vinniefalco/DSPFilters/blob/master/shared/DSPFilters/source
+    # FIXME: type 2 has wrong center frequency?
+    """n-th degree chebyshev low-pass filter cascade
+
+    0 dB at center frequency for type 1.
+    -ripple dB at center frequency for type 2.
+    when ripple=0 and type=1, acts as butterworth."""
+    series = []
+
+    ripple = np.abs(ripple)
+    lin = np.exp(ripple*0.1*np.log(10))
+
+    if ripple != 0:
+        if type == 2:
+            eps = np.sqrt(1/(lin - 1))
+        else:
+            eps = np.sqrt(lin - 1)
+
+        v0 = np.arcsinh(1/eps)/n
+    else:
+        if type == 2:
+            v0 = 0 # allpass?
+        else:
+            v0 = 1 # butterworth
+
+    sinh_v0 = -np.sinh(v0)
+    cosh_v0 = np.cosh(v0)
+
+    fn = np.pi/(2*n)
+    pairs = int(n/2)
+    for i in range(pairs):
+        k = 2*i + 1
+        theta = (k - n)*fn
+        real = sinh_v0*np.cos(theta)
+        imag = cosh_v0*np.sin(theta)
+
+        if type == 2:
+            d2 = real*real + imag*imag
+            im = 1/np.cos(k*fn)
+            real = real/d2
+            imag = imag/d2
+            num = (1, 0, im*im)
+        else:
+            num = (0, 0, 1)
+
+        den = (1, -2*real, real*real + imag*imag)
+
+        # normalize such that 0 Hz is 0 dB
+        den = np.array(den)
+        gain = den[2]/num[2]
+        den /= gain
+
+        series += [(num, den)]
+    if n % 2:
+        real = sinh_v0
+        if type == 2:
+            real = 1/real
+        num = (0, 1)
+        den = (1/-real, 1)
+        series += [(num, den)]
+    return series
+

lib/data.py

@@ -0,0 +1,83 @@
+from .util import *
+from .bq import *
+
+import numpy as np
+
+_bq2a = 1/.76536686473017945
+_bq2b = 1/1.8477590650225735
+_bq2a_bw = isqrt2/_bq2a
+_bq2b_bw = isqrt2/_bq2b
+
+cascades = {
+    '1770': [
+        (1501,      HS2(toA(4), toQ(1)),    0),
+        (38.135457, HP2(0, 0.5003268), np.log10(1.004995)*20),
+    ],
+    # "neon pink"
+    'raw': [
+        nf(LP1,    20,    0,    1,   29),
+        nf(HS1,   800,   12,    1,    0),
+        (   45,   HP2(    0, 1.32), 0.5), # roughly estimates
+        (   45,   HP2(    0, 0.54), 0.5), # a 4-pole butterworth highpass
+        nf(LP2, 14000,    0, 1.33,    0),
+    ],
+    # like neon pink but for feeding into RMS
+    'raw2': [
+        (10000, HP1(0,0),                  26),
+        (  750, HS2(toA(-10), toQ(1.33)),   0),
+        (   45, HP2(0,             1.32), 0.5),
+        (   45, HP2(0,             0.54), 0.5),
+        (14000, LP2(0,        toQ(1.33)),   0),
+        (  250, PE2(toA(3),   toQ(1.33)),  -1),
+        ( 4000, PE2(toA(3),   toQ(1.33)),  -1),
+    ],
+    # loosely based on the equal loudness contour at 60dB or something
+    'raw_ELC': [
+        (   40, HP2(0,        toQ(1.33)),   0),
+        (  400, HP1(0,0),                   6),
+        ( 1400, PE2(toA(-3),  toQ(1.33)),   1),
+        ( 4000, PE2(toA(5),   toQ(1.00)),-1.5),
+        ( 4000, LP2(0,        toQ(1.33)), 1.5),
+    ],
+    # here's the ideas written out:
+    # low (<40) freqs dont contribute much to ears (feel doesnt count.)
+    # high (>14000) freqs are mostly unheard. 14000 is the top for someone with bad hearing.
+    # 750 Hz isn't too painful to the ears, but cutting them would give
+    # overly-produced songs not enough gain to hear vocals, so keep it flat.
+    # we're supposedly less sensitive to 1400 Hz, but i need to
+    # boost high freqs so the RMS even catches it, so keep that flat too.
+    'raw3': [
+        (   40, HP2(0,        toQ(1.33)), 0.0),
+        (  400, HP1(0,                0), 4.5),
+        ( 1400, PE2(toA(-2),  toQ(1.33)), 0.0),
+        ( 8000, PE2(toA(3),   toQ(1.00)), 0.0),
+        (10000, LP2(0,        toQ(0.50)),-0.5),
+    ],
+    'tilt_test': [
+        (10000, HP1(0,0),                  30),
+        ( 1000, HS1(toA(-16),         0), 1.5),
+
+        (   40, HP2(0,        toQ(1.00)), 0.0),
+        (10000, LP1(0,                0), 0.0),
+    ],
+    # tested against your 227 top-rated songs
+    'np2': [
+        nf(LP1,    20,    0,    1,   32),
+        nf(HS1,   800,    9,    1, -4.5),
+        nf(LP2, 14000,    0, 1.33,    0),
+        nf(LP2, 14000,    0, 0.90,    0),
+        nf(HP2,    77,    0, 1.00,    0),
+        nf(LS2,    38,   -9, 1.00,    0),
+        nf(PE2,    64,  4.5, 1.20,    0),
+    ],
+    # side channel
+    'np2s': [
+        nf(LP1,    20,    0,    1,   32),
+        nf(HS1,   800,    9,    1, -4.5),
+        nf(LP2, 14000,    0, 1.33,    0),
+        nf(HP2,    90,    0, 1.11,    0),
+        nf(PE2,    30, -9.5, 1.00,    0),
+        #(17500, LP2(0,            _bq2a),   0),
+        #(17500, LP2(0,            _bq2b),   0),
+    ],
+}

lib/fft.py

@@ -0,0 +1,48 @@
+from . import rfft
+
+import numpy as np
+import scipy.signal as sig
+
+def magnitudes_window_setup(s, size=8192):
+    L = s.shape[0]
+    overlap = 0.661
+    step = np.ceil(size*(1 - overlap))
+    segs = np.ceil(L/step)
+    return step, segs
+
+def magnitudes(s, size=8192):
+    import scipy.linalg as linalg
+
+    step, segs = magnitudes_window_setup(s, size)
+
+    L = s.shape[0]
+
+    # blindly pad with zeros for friendlier ffts and overlapping
+    z = np.zeros(size)
+    s = np.hstack((s, z))
+
+    win_size = size
+
+    win = sig.blackmanharris(win_size)
+    win /= linalg.norm(win)
+
+    count = 0
+    for i in range(0, L - 1, int(step)):
+        windowed = s[i:i+win_size]*win
+        power = np.abs(rfft(windowed, size))**2
+        # this scraps the nyquist value to get exactly size outputs
+        yield power[0:size]
+        count += 1
+
+    #assert(segs == count)
+
+def averfft(s, size=8192):
+    """calculates frequency magnitudes by fft and averages them together."""
+    step, segs = magnitudes_window_setup(s, size)
+
+    avg = np.zeros(size)
+    for power in magnitudes(s, size):
+        avg += power/segs
+
+    avg_db = 10*np.log10(avg)
+    return avg_db

lib/planes.py

@@ -0,0 +1,83 @@
+from . import tau
+
+import numpy as np
+import sympy as sym
+
+def zcgen_py(n, d):
+    zcs = np.zeros(d + 1)
+    zcs[0] = 1
+    for _ in range(n):
+        for i in range(d, 0, -1):
+            zcs[i] -= zcs[i - 1]
+    for _ in range(d - n):
+        for i in range(d, 0, -1):
+            zcs[i] += zcs[i - 1]
+    return zcs
+
+def zcgen_sym(n, d):
+    z = sym.symbols('z')
+    expr = sym.expand((1 - z**-1)**n*(1 + z**-1)**(d - n))
+    coeffs = expr.equals(1) and [1] or expr.as_poly().all_coeffs()
+    return coeffs[::-1]
+
+def s2z_two(b,a,fc,srate,gain=1):
+    """
+    converts s-plane coefficients to z-plane for digital usage.
+    hard-coded for 3 coefficients.
+    """
+    if (len(b) < 3):
+        b = (b[0], b[1], 0)
+    if (len(a) < 3):
+        a = (a[0], a[1], 0)
+    w0 = tau*fc/srate
+    cw = np.cos(w0)
+    sw = np.sin(w0)
+    zb = np.array((
+           b[2]*(1 - cw) + b[0]*(1 + cw) + b[1]*sw,
+        2*(b[2]*(1 - cw) - b[0]*(1 + cw)),
+           b[2]*(1 - cw) + b[0]*(1 + cw) - b[1]*sw,
+    ))
+    za = np.array((
+           a[2]*(1 - cw) + a[0]*(1 + cw) + a[1]*sw,
+        2*(a[2]*(1 - cw) - a[0]*(1 + cw)),
+           a[2]*(1 - cw) + a[0]*(1 + cw) - a[1]*sw,
+    ))
+    return zb*gain, za
+
+def s2z1(w0, s, d):
+    """
+    s: array of s-plane coefficients (num OR den, not both)
+    d: degree (array length - 1)
+    returns output array of size d + 1
+    """
+    y = np.zeros(d + 1)
+    sw = np.sin(w0)
+    cw = np.cos(w0)
+    for n in range(d + 1):
+        zcs = zcgen(d - n, d)
+        trig = sw**n/(cw + 1)**(n - 1)
+        for i in range(d + 1):
+            y[i] += trig*zcs[i]*s[n]
+    return y
+
+def s2z_any(b, a, fc, srate, gain=1, d=-1):
+    """
+    converts s-plane coefficients to z-plane for digital usage.
+    supports any number of coefficients; b or a will be padded accordingly.
+    additional padding can be specified with d.
+    """
+    cs = max(len(b), len(a), d + 1)
+    sb = np.zeros(cs)
+    sa = np.zeros(cs)
+    sb[:len(b)] = b
+    sa[:len(a)] = a
+    w0 = tau*fc/srate
+    zb = s2z1(w0, sb, cs - 1)
+    za = s2z1(w0, sa, cs - 1)
+    return zb*gain, za
+
+# set our preference. zcgen_py is 1000+ times faster than zcgen_sym
+zcgen = zcgen_py
+
+# s2z_any is only ~2.4 times slower than s2z_two and allows for filters of any degree
+s2z = s2z_any

lib/plot.py

@@ -0,0 +1,27 @@
+import matplotlib.pyplot as plt
+from matplotlib import ticker
+
+def response_setup(ax, ymin=-24, ymax=24, yL=ticker.AutoMinorLocator(3)):
+    ax.set_xlim(20, 20000)
+    ax.set_ylim(ymin, ymax)
+    #ax.set_yticks(np.arange(ymin, ymax + 1, 6))
+    ax.set_yticks(tuple(range(ymin, ymax + 1, 6)))
+    ax.yaxis.set_minor_locator(yL)
+    ax.grid(True, 'both')
+    ax.set_xlabel('frequency (Hz)')
+    ax.set_ylabel('magnitude (dB)')
+
+def cleanplot():
+    fig, ax = plt.subplots()
+    fig.set_size_inches((16, 16))
+    ax.set_axis_off()
+    ax.set_position([0,0,1,1])
+    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

lib/smoothfft.py

@@ -0,0 +1,60 @@
+from . import xsp
+import numpy as np
+
+def smoothfft(xs, ys, bw=1, precision=512):
+    """performs log-lin smoothing on magnitude data,
+    generally from the output of averfft."""
+    # TODO: option to extrapolate (pad) fft data
+    xs2 = xsp(precision)
+    ys2 = np.zeros(precision)
+    log_xs = np.log(xs)
+    for i, x in enumerate(xs2):
+        dist = np.exp(np.abs(log_xs - np.log(x + 1e-35)))
+        window = np.maximum(0, 1 - (dist - bw))
+        # at this point you could probably
+        # normalize our *triangular* window to 0-1
+        # and transform it into *another* windowing function
+        wsum = np.sum(window)
+        ys2[i] = np.sum(ys*window/wsum)
+    return xs2, ys2
+
+def smoothfft2(xs, ys, bw=1, precision=512, compensate=True):
+    """performs log-lin smoothing on magnitude data,
+    generally from the output of averfft."""
+    xs2 = xsp(precision)
+    ys2 = np.zeros(precision)
+    log2_xs2 = np.log2(xs2)
+    for i, x in enumerate(xs):
+        #dist = np.abs(np.log2(xs2/(x + 1e-35)))/bw
+        dist = np.abs(log2_xs2 - np.log2(x + 1e-35))/bw
+        #window = np.maximum(0, 1 - dist) # triangle window
+        window = np.exp(-dist**2/(0.5/2)) # gaussian function (non-truncated)
+        ys2 += ys[i]*window
+    if compensate:
+        _, temp = smoothfft2(xs, np.ones(len(xs)), bw=bw, precision=precision, compensate=False)
+        ys2 /= temp
+    return xs2, ys2
+
+def smoothfft3(xs, ys, bw=1, precision=1024):
+    # actually this will never work...
+    # you need to go back to smoothfft2,
+    # which technically works as-designed,
+    # and fix the compensation to work with widely-spaced data.
+    raise Exception("smoothfft3 is broken.")
+    xs2 = xsp(precision)
+    ys2 = np.zeros(precision)
+    step = (xs[1] - xs[0])
+    if True:
+        for i, x in enumerate(xs):
+            dist = np.abs(xs2 - x)
+            bw2 = x*bw/2
+            window = np.maximum(0, 1 - dist/bw2)
+            #window = np.minimum(1, np.maximum(0, 1 - (dist - bw)))
+            ys2 += ys[i]*window
+    else:
+        for i, x2 in enumerate(xs2):
+            dist = np.abs(xs - x2)
+            window = np.maximum(0, 1 - (dist/step/bw))
+            wsum = np.sum(window)
+            ys2[i] = np.sum(ys*window/wsum)
+    return xs2, ys2

lib/svf.py

@@ -0,0 +1,70 @@
+from . import tau, unwarp
+
+import numpy as np
+
+# via http://nbviewer.ipython.org/urls/music-synthesizer-for-android.googlecode.com/git/lab/Second%20order%20sections%20in%20matrix%20form.ipynb
+def svf_core(w0, Q, m, shelfA=1, gain=1):
+    # TODO: implement constant gain parameter
+    g = unwarp(w0)*shelfA
+    a1 = 1/(1 + g*(g + 1/Q))
+    a2 = g*a1
+    a3 = g*a2
+    A = np.array([[2*a1 - 1, -2*a2], [2*a2, 1 - 2*a3]])
+    B = np.array([2*a2, 2*a3])
+    v0 = np.array([1, 0, 0])
+    v1 = np.array([a2, a1, -a2])
+    v2 = np.array([a3, a2, 1 - a3])
+    C = v0*m[0] + v1*m[1] + v2*m[2]
+    return A, B, C
+
+LP2S = lambda A, Q: (Q, [0,    0,  1], 1)
+BP2S = lambda A, Q: (Q, [0,    1,  0], 1)
+HP2S = lambda A, Q: (Q, [1, -1/Q, -1], 1)
+#AP2S = lambda A, Q:
+#BP2aS = lambda A, Q:
+#BP2bS = lambda A, Q:
+NO2S = lambda A, Q: (Q, [1, -1/Q,  0], 1)
+PE2S = lambda A, Q: (Q*A, [1, (A*A - 1)/A/Q, 0], 1)
+LS2S = lambda A, Q: (Q, [1,     (A - 1)/Q, A*A - 1], 1/np.sqrt(A))
+HS2S = lambda A, Q: (Q, [A*A, (1 - A)*A/Q, 1 - A*A],   np.sqrt(A))
+
+# actual peaking filter (not a bell?)
+#PE2S = lambda A, Q: ([1, -1/Q, -2], 1)
+# original uncompensated
+#PE2S = lambda A, Q: (Q, [1,   (A*A - 1)/Q,       0], 1)
+#LS2S = lambda A, Q: (Q, [1,     (A - 1)/Q, A*A - 1], 1/np.sqrt(A))
+#HS2S = lambda A, Q: (Q, [A*A, (A - A*A)/Q, 1 - A*A], 1/np.sqrt(A))
+
+gen_filters_svf = lambda cascade, srate: [
+    svf_core(tau*f[0]/srate, *f[1], gain=10**(f[2]/20)) for f in cascade
+]
+
+def svf_run(svf, xs):
+    A, B, C = svf
+    result = []
+    y = np.zeros(2) # filter memory
+    for x in xs:
+        result.append(np.dot(C, np.concatenate([[x], y])))
+        y = B*x + np.dot(A, y)
+    return np.array(result)
+
+def svf_mat(svf):
+    A, B, C = svf
+    AA = np.dot(A, A)
+    AB = np.dot(A, B)
+    CA = np.dot(C[1:], A)
+    cb = np.dot(C[1:], B)
+    return np.array([[ C[0],    0,     C[1],     C[2]],
+                     [   cb, C[0],    CA[0],    CA[1]],
+                     [AB[0], B[0], AA[0][0], AA[0][1]],
+                     [AB[1], B[1], AA[1][0], AA[1][1]]])
+
+def svf_run4(mat, xs):
+    assert(len(xs) % 2 == 0)
+    out = np.zeros(len(xs))
+    y = np.zeros(4) # filter memory
+    for i in range(0, len(xs), 2):
+        y[0:2] = xs[i:i+2]
+        y = np.dot(mat, y)
+        out[i:i+2] = y[0:2]
+    return out

lib/sweeps.py

@@ -0,0 +1,54 @@
+from . import tau
+
+import numpy as np
+
+def sweep(amp, length, begin=20, end=20480, method='linear'):
+    method = method or 'linear'
+    xs = np.arange(length)/length
+    if method in ('linear', 'quadratic', 'logarithmic', 'hyperbolic'):
+        ys = amp*sig.chirp(xs, begin, 1, end, method=method)
+    elif method is 'sinesweep':
+        ang = lambda f: tau*f
+        # because xs ranges from 0:1, length is 1 and has been simplified out
+        domain = np.log(ang(end)/ang(begin))
+        ys = amp*np.sin(ang(begin)/domain*(np.exp(xs*domain) - 1))
+    return ys
+
+def tsp(N, m=0.5):
+    """
+        OATSP(Optimized Aoshima's Time-Stretched Pulse) generator
+        x = tsp( N, m )
+        N : length of the time-stretched pulse
+        m : ratio of the swept sine  (0 < m < 1)
+
+        Author(s): Seigo UTO 8-23-95
+
+        Reference:
+           Yoiti SUZUKI, Futoshi ASANO, Hack-Yoon KIM and Toshio SONE,
+           "Considerations on the Design of Time-Stretched Pulses,"
+               Techical Report of IEICE, EA92-86(1992-12)
+    """
+    # http://www.sound.sie.dendai.ac.jp/dsp/e-21.html
+
+    if m < 0 or m > 1:
+        raise Exception("sdfgsdfgsdg")
+
+    if N < 0:
+        raise Exception("The number of length must be the positive number")
+
+    NN = 2**np.floor(np.log2(N)) # nearest
+    NN2 = NN//2
+    M = np.round(NN2*m)
+
+    nn2 = np.arange(NN2 + 1)**2
+
+    j = np.complex(0, 1)
+
+    H = np.exp(j*4*M*np.pi*nn2/NN**2)
+    H2 = np.hstack([H, np.conj(H[1:NN2][::-1])])
+
+    x = np.fft.ifft(H2)
+    x = np.hstack([x[NN2 - M:NN + 1], x[0:NN2 - M + 1]])
+    x = np.hstack([x.real, np.zeros(1, N - NN)])
+
+    return x

lib/util.py

@@ -0,0 +1,53 @@
+import sys
+import numpy as np
+import scipy.signal as sig
+
+dummy = lambda *args, **kwargs: None
+lament = lambda *args, **kwargs: print(*args, file=sys.stderr, **kwargs)
+
+toLK = lambda x: -0.691 + 10*np.log10(x)
+isqrt2 = 1/np.sqrt(2)
+toQ = lambda bw: isqrt2/bw
+toA = lambda db: 10**(db/40)
+
+tau = 2*np.pi
+unwarp = lambda w: np.tan(w/2)
+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]
+
+def xsp(precision=4096):
+    """create #precision log-spaced points from 20 to 20480 Hz"""
+    # i opt not to use steps or linspace here,
+    # as the current method is less error-prone for me.
+    xs = np.arange(0,precision)/precision
+    return 20*1024**xs
+
+def blocks(a, step, size=None):
+    """break an iterable into chunks"""
+    if size is None:
+        size = step
+    for start in range(0, len(a), step):
+        end = start + size
+        if end > len(a):
+            break
+        yield a[start:end]
+
+def convolve_each(s, fir, mode=None, axis=0):
+    return np.apply_along_axis(lambda s: sig.fftconvolve(s, fir, mode), axis, s)
+
+def count_channels(s):
+    if len(s.shape) < 2:
+        return 1
+    return s.shape[1]
+
+def monoize(s):
+    """mixes an n-channel signal down to one channel.
+    technically, it averages a 2D array to be 1D.
+    existing mono signals are passed through unmodified."""
+    channels = count_channels(s)
+    if channels != 1:
+        s = np.sum(s, 1)
+        s /= channels
+    return s

lib/wav.py

@@ -0,0 +1,32 @@
+import numpy as np
+from . import lament
+
+# TODO: don't use wavfile, it breaks on perfectly good files
+import scipy.io.wavfile as wav
+import ewave
+
+def wav_smart_read(fn):
+    lament('DEPRECATED: wav_smart_read; use wav_read instead')
+    srate, s = wav.read(fn)
+    if s.dtype != np.float64:
+        bits = s.dtype.itemsize*8
+        s = np.asfarray(s)/2**(bits - 1)
+    return srate, s
+
+def wav_smart_write(fn, srate, s):
+    lament('DEPRECATED: wav_smart_write')
+    si = np.zeros_like(s, dtype='int16')
+    bits = si.dtype.itemsize*8
+    si += np.clip(s*2**(bits - 1), -32768, 32767)
+    wav.write(fn, srate, si)
+
+def wav_read(fn):
+    with ewave.open(fn) as f:
+        s = f.read()
+        srate = f.sampling_rate
+    if s.dtype == np.float32:
+        s = np.asfarray(s)
+    elif s.dtype != np.float64:
+        bits = s.dtype.itemsize*8
+        s = np.asfarray(s)/2**(bits - 1)
+    return s, srate