commit bf22e5a326fad4c37a7c4007aca5bb88c31b67c0 Author: Connor Olding Date: Sun Oct 18 23:06:39 2015 -0700 update 1 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e89fb22 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +*.pyc +__pycache__/* +lib/__main__.py +version diff --git a/README.md b/README.md new file mode 100644 index 0000000..b3b79b8 --- /dev/null +++ b/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. diff --git a/autoupdate b/autoupdate new file mode 100755 index 0000000..595ed63 --- /dev/null +++ b/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 diff --git a/lib/__init__.py b/lib/__init__.py new file mode 100644 index 0000000..8e2edb8 --- /dev/null +++ b/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('_')] diff --git a/lib/bq.py b/lib/bq.py new file mode 100644 index 0000000..3e48b31 --- /dev/null +++ b/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)) diff --git a/lib/bs.py b/lib/bs.py new file mode 100644 index 0000000..d8bf299 --- /dev/null +++ b/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 diff --git a/lib/butterworth.py b/lib/butterworth.py new file mode 100644 index 0000000..8472f76 --- /dev/null +++ b/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 + diff --git a/lib/data.py b/lib/data.py new file mode 100644 index 0000000..403fb58 --- /dev/null +++ b/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), + ], +} diff --git a/lib/fft.py b/lib/fft.py new file mode 100644 index 0000000..c2a7781 --- /dev/null +++ b/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 diff --git a/lib/planes.py b/lib/planes.py new file mode 100644 index 0000000..8bcd2f4 --- /dev/null +++ b/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 diff --git a/lib/plot.py b/lib/plot.py new file mode 100644 index 0000000..813a8f7 --- /dev/null +++ b/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 diff --git a/lib/smoothfft.py b/lib/smoothfft.py new file mode 100644 index 0000000..2eec514 --- /dev/null +++ b/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 diff --git a/lib/svf.py b/lib/svf.py new file mode 100644 index 0000000..43bd252 --- /dev/null +++ b/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 diff --git a/lib/sweeps.py b/lib/sweeps.py new file mode 100644 index 0000000..4881c85 --- /dev/null +++ b/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 diff --git a/lib/util.py b/lib/util.py new file mode 100644 index 0000000..6b40460 --- /dev/null +++ b/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 diff --git a/lib/wav.py b/lib/wav.py new file mode 100644 index 0000000..4cc7ca1 --- /dev/null +++ b/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