diff --git a/lib/__init__.py b/lib/__init__.py index 352aeb5..f733fb3 100644 --- a/lib/__init__.py +++ b/lib/__init__.py @@ -12,141 +12,13 @@ from .bs import * from .cepstrum import * from .windowing import * from .piir import * - -import numpy as np - -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)) - -def test_filter(ff, A=toA(12), Q=toQ(1), **kwargs): - test_filter_raw(ff(A, Q), **kwargs) - -def neonpink(xs): - lament("neonpink(): DEPRECATED; use tilter2(xs, 'raw') instead.") - return tilter2(xs, 'raw') - -def c_render(cascade, precision=4096): - # TODO: deprecate in favor of tilter2 - xs = xsp(precision) - return xs, tilter2(xs, cascade) - -def c_render2(xs, cascade, phase=False): - """c_render optimized and specifically for first/second-order filters""" - if phase: - return c_render3(xs, cascade, mode='phase') - else: - return c_render3(xs, cascade, mode='magnitude') - -def c_render3(xs, cascade, mode='magnitude'): - """c_render optimized and specifically for first/second-order filters""" - import numexpr as ne - j = np.complex(0, 1) - - # obviously this could be extended to higher orders - eq2 = '(b0 + b1*s + b2*s**2)/(a0 + a1*s + a2*s**2)' - eq1 = '(b0 + b1*s)/(a0 + a1*s)' - - if mode == 'magnitude': - fmt = 'real(log10(abs({})**2)*10 + gain)' - elif mode == 'phase' or mode == 'group delay': - fmt = '-arctan2(imag({0}), real({0}))' # gross - else: - raise Exception("c_render3(): unknown mode: {}".format(mode)) - - ys = np.zeros(len(xs)) - for f in cascade: - freq, ba, gain = f - b, a = ba - if len(b) == 3 and len(a) == 3: - eq = fmt.format(eq2) - b2, b1, b0 = b - a2, a1, a0 = a - elif len(b) == 2 and len(a) == 2: - eq = fmt.format(eq1) - b1, b0 = b - a1, a0 = a - else: - raise Exception("incompatible cascade; consider using c_render instead") - - if mode == 'group delay': - # approximate derivative of phase by slope of tangent line - step = 2**-8 - fa = freq - step - fb = freq + step - - s = xs/fa*j - ya = ne.evaluate(eq) - s = xs/fb*j - yb = ne.evaluate(eq) - - slope = (yb - ya)/(2*step) - ys += -slope/(xs/freq*tau) - else: - s = xs/freq*j - ys += ne.evaluate(eq) - if mode == 'phase': - ys = degrees_clamped(ys) - return ys - -def firize(xs, ys, n=4096, srate=44100, ax=None): - import scipy.signal as sig - if ax: - ax.semilogx(xs, ys, label='desired') - xf = xs/srate*2 - yg = 10**(ys/20) - - xf = np.r_[0, xf, 1] - yg = np.r_[0, yg, yg[-1]] - - b = sig.firwin2(n, xf, yg, antisymmetric=True) - - if ax: - _, ys = sig.freqz(b, worN=xs/srate*tau) - ys = 20*np.log10(np.abs(ys)) - ax.semilogx(xs, ys, label='FIR ({} taps)'.format(n)) - ax.legend(loc=8) - - return b - -def tilter(xs, ys, tilt): - """tilts a magnitude plot by some decibels, or by equalizer curve.""" - lament("tilter(): DEPRECATED; use ys -= tilter2(xs, tilt) instead.") - return xs, ys - tilter2(xs, tilt) - -def tilter2(xs, tilt): # TODO: rename - noise = np.zeros(xs.shape) - if isinstance(tilt, str) and tilt in cascades: - tilt = cascades[tilt] - if isinstance(tilt, list): - c = [makemag(*f) for f in 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)) - return noise - +from .mag import * from .plotwav import * + # 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('_')] +__all__ = [ + o for o in locals() + if type(o) != 'module' and not o.startswith('_')] diff --git a/lib/bq.py b/lib/bq.py index ef29c69..3c81a32 100644 --- a/lib/bq.py +++ b/lib/bq.py @@ -4,39 +4,79 @@ import scipy.signal as sig from .util import * from .planes import s2z -bq_run = lambda bq, xs: sig.lfilter(*bq, x=xs, axis=0) -nfba = lambda b, a: (1/tau, (b, a), 0) -nf = lambda t, f, g, bw, mg: (f, t(toA(g), toQ(bw)), mg) +# PEP 8 fucking destroyed this file. I'm sorry. + + +def bq_run(bq, xs): + return sig.lfilter(*bq, x=xs, axis=0) + + +def nfba(b, a): + return (1/tau, (b, a), 0) + + +def nf(t, f, g, bw, mg): + return (f, t(toA(g), toQ(bw)), mg) + + +def LP1(A, Q): + return ((0, 1), (1, 1)) + + +def HP1(A, Q): + return ((1, 0), (1, 1)) + + +def LS1(A, Q): + return ((1, A), (1, 1/A)) + + +def HS1(A, Q): + return ((A, 1), (1/A, 1)) -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)) # patterns observed, in case some simplification could be done: # 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)) +def LP2(A, Q): + return ((0, 0, 1), (1, 1/Q, 1)) -gen_filters = lambda cascade, srate: [ - s2z(*f[1], fc=f[0], srate=srate, gain=10**(f[2]/20)) for f in cascade -] + +def HP2(A, Q): + return ((1, 0, 0), (1, 1/Q, 1)) + + +def PE2(A, Q): + return ((1, A/Q, 1), (1, 1/A/Q, 1)) + + +def AP2(A, Q): + return ((1, -1/Q, 1), (1, 1/Q, 1)) + + +def BP2a(A, Q): + return ((0, -A/Q, 0), (1, 1/A/Q, 1)) + + +def BP2b(A, Q): + return ((0, -A*A/Q, 0), (1, 1/Q, 1)) + + +def NO2(A, Q): + return ((1, 0, 1), (1, 1/Q, 1)) + + +def LS2(A, Q): + return ((1, np.sqrt(A)/Q, A), (1, 1/np.sqrt(A)/Q, 1/A)) + + +def HS2(A, Q): + return ((A, np.sqrt(A)/Q, 1), (1/A, 1/np.sqrt(A)/Q, 1)) + + +def gen_filters(cascade, srate): + return [ + s2z(*f[1], fc=f[0], srate=srate, gain=10**(f[2]/20)) for f in cascade + ] diff --git a/lib/bs.py b/lib/bs.py index 1cf1ede..1f55f96 100644 --- a/lib/bs.py +++ b/lib/bs.py @@ -3,6 +3,7 @@ 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: @@ -10,7 +11,7 @@ def BS1770_3(s, srate, filters=None, window=0.4, overlap=0.75, sf = np.copy(s) for f in filters: - if len(f) is 2: # dumb but effective + if len(f) is 2: # dumb way to tell what type we're given. sf = bq_run(f, sf) else: sf = convolve_each(sf, f, 'same') @@ -35,6 +36,7 @@ def BS1770_3(s, srate, filters=None, window=0.4, overlap=0.75, 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) @@ -47,25 +49,25 @@ def BS_plot(ys, g10=None, g70=None, threshold=None, fig=None, ax=None): if ax is None: ax = fig.gca() - if False: # histogram + 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) + fig.set_size_inches(10, 4) xs = np.arange(len(ys)) - #ax.plot(xs, ys, color='#066ACF', linestyle=':', marker='d', markersize=2) + # 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() + fig.set_size_inches(12, 5) + # _, _, ymin, _ = ax.axis() if threshold: ax.axhspan(-70, threshold, facecolor='r', alpha=1/5) if g10: @@ -75,6 +77,7 @@ def BS_plot(ys, g10=None, g70=None, threshold=None, fig=None, ax=None): return fig, ax + def normalize(s, srate): """performs BS.1770-3 normalization and returns inverted gain.""" db = BS1770_3(s, srate) diff --git a/lib/cepstrum.py b/lib/cepstrum.py index 3fb1095..54eb129 100644 --- a/lib/cepstrum.py +++ b/lib/cepstrum.py @@ -1,12 +1,19 @@ import numpy as np from .util import pad2 -# fast cepstrum and inverted fast cepstrum -fcs = lambda s: np.fft.ifft(np.log(np.fft.fft(s))) -ifcs = lambda s: np.fft.fft(np.exp(np.fft.ifft(s))) -# magnitude -mcs = lambda s: (np.abs(np.fft.ifft(np.log(np.abs(np.fft.fft(s))**2)))**2)[:len(s)//2] +def fcs(s): # fast cepstrum + return np.fft.ifft(np.log(np.fft.fft(s))) + + +def ifcs(s): # inverted fast cepstrum + return np.fft.fft(np.exp(np.fft.ifft(s))) + + +def mcs(s): # magnitude + return (np.abs(np.fft.ifft(np.log(np.abs(np.fft.fft(s))**2)))**2 + )[:len(s)//2] + def clipdb(s, cutoff=-100): as_ = np.abs(s) @@ -16,6 +23,7 @@ def clipdb(s, cutoff=-100): thresh = mas*10**(cutoff/20) return np.where(as_ < thresh, thresh, s) + def fold(r): # via https://ccrma.stanford.edu/~jos/fp/Matlab_listing_fold_m.html # Fold left wing of vector in "FFT buffer format" onto right wing @@ -33,6 +41,7 @@ def fold(r): rw = np.r_[r[0], rf, np.zeros(n-nt-1)] return rw + def minphase(s, pad=True, os=False): # via https://ccrma.stanford.edu/~jos/fp/Matlab_listing_mps_m.html # TODO: actual oversampling diff --git a/lib/data.py b/lib/data.py index 45e7cc3..21037f5 100644 --- a/lib/data.py +++ b/lib/data.py @@ -14,18 +14,20 @@ cascades = { (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), # i don't use the exact _bq2 coeffecients here for legacy reasons - ( 45, HP2( 0, 1.32), 0.5), # roughly estimates - ( 45, HP2( 0, 0.54), 0.5), # a 4-pole butterworth highpass + ( 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), + (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), @@ -33,14 +35,16 @@ cascades = { ( 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), + ( 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. @@ -55,6 +59,7 @@ cascades = { ( 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), @@ -62,6 +67,7 @@ cascades = { ( 40, HP2(0, toQ(1.00)), 0.0), (10000, LP1(0, 0), 0.0), ], + # average curve of my 227 favorite songs 'np2': [ nf(LP1, 20, 0, 1, 32), @@ -72,6 +78,7 @@ cascades = { nf(LS2, 38, -9, 1.00, 0), nf(PE2, 64, 4.5, 1.20, 0), ], + # same but for the side channel 'np2s': [ nf(LP1, 20, 0, 1, 32), @@ -79,7 +86,5 @@ cascades = { 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 index 90fd06b..bb472bd 100644 --- a/lib/fft.py +++ b/lib/fft.py @@ -3,6 +3,7 @@ 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 @@ -10,6 +11,7 @@ def magnitudes_window_setup(s, size=8192): segs = np.ceil(L/step) return step, segs + def magnitudes(s, size=8192): import scipy.linalg as linalg @@ -34,7 +36,8 @@ def magnitudes(s, size=8192): yield power[0:size] count += 1 - #assert(segs == count) # this is probably no good in a generator + # assert(segs == count) # this is probably no good in a generator + def averfft(s, size=8192): """calculates frequency magnitudes by fft and averages them together.""" diff --git a/lib/mag.py b/lib/mag.py new file mode 100644 index 0000000..d39a284 --- /dev/null +++ b/lib/mag.py @@ -0,0 +1,145 @@ +from . import toA, toQ, cascades + +import numpy as np + + +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)) + + +def test_filter(ff, A=toA(12), Q=toQ(1), **kwargs): + test_filter_raw(ff(A, Q), **kwargs) + + +def neonpink(xs): + lament("neonpink(): DEPRECATED; use tilter2(xs, 'raw') instead.") + return tilter2(xs, 'raw') + + +def c_render(cascade, precision=4096): + # TODO: deprecate in favor of tilter2 + xs = xsp(precision) + return xs, tilter2(xs, cascade) + + +def c_render2(xs, cascade, phase=False): + """c_render optimized and specifically for first/second-order filters""" + if phase: + return c_render3(xs, cascade, mode='phase') + else: + return c_render3(xs, cascade, mode='magnitude') + + +def c_render3(xs, cascade, mode='magnitude'): + """c_render optimized and specifically for first/second-order filters""" + import numexpr as ne + j = np.complex(0, 1) + + # obviously this could be extended to higher orders + eq2 = '(b0 + b1*s + b2*s**2)/(a0 + a1*s + a2*s**2)' + eq1 = '(b0 + b1*s)/(a0 + a1*s)' + + if mode == 'magnitude': + fmt = 'real(log10(abs({})**2)*10 + gain)' + elif mode == 'phase' or mode == 'group delay': + fmt = '-arctan2(imag({0}), real({0}))' # gross + else: + raise Exception("c_render3(): unknown mode: {}".format(mode)) + + ys = np.zeros(len(xs)) + for f in cascade: + freq, ba, gain = f + b, a = ba + if len(b) == 3 and len(a) == 3: + eq = fmt.format(eq2) + b2, b1, b0 = b + a2, a1, a0 = a + elif len(b) == 2 and len(a) == 2: + eq = fmt.format(eq1) + b1, b0 = b + a1, a0 = a + else: + raise Exception( + "incompatible cascade; consider using c_render instead") + + if mode == 'group delay': + # approximate derivative of phase by slope of tangent line + step = 2**-8 + fa = freq - step + fb = freq + step + + s = xs/fa*j + ya = ne.evaluate(eq) + s = xs/fb*j + yb = ne.evaluate(eq) + + slope = (yb - ya)/(2*step) + ys += -slope/(xs/freq*tau) + else: + s = xs/freq*j + ys += ne.evaluate(eq) + if mode == 'phase': + ys = degrees_clamped(ys) + return ys + + +def firize(xs, ys, n=4096, srate=44100, ax=None): + import scipy.signal as sig + if ax: + ax.semilogx(xs, ys, label='desired') + xf = xs/srate*2 + yg = 10**(ys/20) + + xf = np.r_[0, xf, 1] + yg = np.r_[0, yg, yg[-1]] + + b = sig.firwin2(n, xf, yg, antisymmetric=True) + + if ax: + _, ys = sig.freqz(b, worN=xs/srate*tau) + ys = 20*np.log10(np.abs(ys)) + ax.semilogx(xs, ys, label='FIR ({} taps)'.format(n)) + ax.legend(loc=8) + + return b + + +def tilter(xs, ys, tilt): + """tilts a magnitude plot by some decibels, or by equalizer curve.""" + lament("tilter(): DEPRECATED; use ys -= tilter2(xs, tilt) instead.") + return xs, ys - tilter2(xs, tilt) + + +def tilter2(xs, tilt): # TODO: rename + noise = np.zeros(xs.shape) + if isinstance(tilt, str) and tilt in cascades: + tilt = cascades[tilt] + if isinstance(tilt, list): + c = [makemag(*f) for f in 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)) + return noise diff --git a/lib/nsf.py b/lib/nsf.py index e0b02b5..c37759e 100644 --- a/lib/nsf.py +++ b/lib/nsf.py @@ -2,8 +2,10 @@ import numpy as np + def LPB(n): - # via https://github.com/vinniefalco/DSPFilters/blob/master/shared/DSPFilters/source + # via: + # https://github.com/vinniefalco/DSPFilters/blob/master/shared/DSPFilters/source """n-th order butterworth low-pass filter cascade -3 dB at center frequency.""" @@ -24,8 +26,10 @@ def LPB(n): series += [(num, den)] return series + def LPC(n, ripple, type=1): - # via https://github.com/vinniefalco/DSPFilters/blob/master/shared/DSPFilters/source + # via: + # https://github.com/vinniefalco/DSPFilters/blob/master/shared/DSPFilters/source # FIXME: type 2 has wrong center frequency? """n-th order chebyshev low-pass filter cascade @@ -46,9 +50,9 @@ def LPC(n, ripple, type=1): v0 = np.arcsinh(1/eps)/n else: if type == 2: - v0 = 0 # allpass? + v0 = 0 # allpass? else: - v0 = 1 # butterworth + v0 = 1 # butterworth sinh_v0 = -np.sinh(v0) cosh_v0 = np.cosh(v0) diff --git a/lib/piir.py b/lib/piir.py index f6d8195..dcca681 100644 --- a/lib/piir.py +++ b/lib/piir.py @@ -57,6 +57,7 @@ halfband_c['olli'] = [ 0.9987488452737**2, ] + class Halfband: def __init__(self, c='olli'): self.x = np.zeros(4) @@ -94,10 +95,10 @@ class Halfband: sign = 1 if mode == 'hilbert': - #y[n] = c*(x[n] + y[n-2]) - x[n-2] + # y[n] = c*(x[n] + y[n-2]) - x[n-2] pass elif mode == 'filter': - #y[n] = c*(x[n] - y[n-2]) + x[n-2] + # y[n] = c*(x[n] - y[n-2]) + x[n-2] sign = -1 in2 = self.x[2] diff --git a/lib/planes.py b/lib/planes.py index 343e698..e2e2b8e 100644 --- a/lib/planes.py +++ b/lib/planes.py @@ -2,6 +2,7 @@ from . import tau import numpy as np + # implements the modified bilinear transform: # s <- 1/tan(w0/2)*(1 - z^-1)/(1 + z^-1) # this requires the s-plane coefficients to be frequency-normalized, @@ -20,6 +21,7 @@ def zcgen_py(n, d): zcs[i] += zcs[i - 1] return zcs + def zcgen_sym(n, d): import sympy as sym z = sym.symbols('z') @@ -27,6 +29,7 @@ def zcgen_sym(n, d): 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. @@ -40,17 +43,18 @@ def s2z_two(b, a, fc, srate, gain=1): 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, + (b[2]*(1 - cw) + b[0]*(1 + cw) + b[1]*sw), + (b[2]*(1 - cw) - b[0]*(1 + cw)) * 2, + (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, + (a[2]*(1 - cw) + a[0]*(1 + cw) + a[1]*sw), + (a[2]*(1 - cw) - a[0]*(1 + cw)) * 2, + (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) @@ -67,6 +71,7 @@ def s2z1(w0, s, d): 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. @@ -83,8 +88,10 @@ def s2z_any(b, a, fc, srate, gain=1, d=-1): za = s2z1(w0, sa, cs - 1) return zb*gain, za -# set our preference. zcgen_py is 1000+ times faster than zcgen_sym + +# 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_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 index 3b1d889..d3a12c8 100644 --- a/lib/plot.py +++ b/lib/plot.py @@ -1,6 +1,7 @@ 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) @@ -10,6 +11,7 @@ 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) @@ -19,22 +21,26 @@ def phase_response_setup(ax, div=12, yL=ticker.AutoMinorLocator(2)): ax.set_xlabel('frequency (Hz)') ax.set_ylabel('phase (degrees)') + def cleanplot(): fig, ax = plt.subplots() ax.set_axis_off() - ax.set_position([0,0,1,1]) + ax.set_position([0, 0, 1, 1]) return fig, ax + def new_response(*args, **kwargs): fig, ax = plt.subplots() response_setup(ax, *args, **kwargs) return fig, ax + def new_phase_response(*args, **kwargs): fig, ax = plt.subplots() phase_response_setup(ax, *args, **kwargs) return fig, ax + def new_bode(magnitude_offset=0): fig, ax1 = plt.subplots() ax2 = ax1.twinx() @@ -51,8 +57,8 @@ def new_bode(magnitude_offset=0): for tl in ax2.get_yticklabels(): tl.set_color(cc[1]) - #ax1.hlines(0, 20, 40, linewidth=0.5, color=cc[0]) - #ax2.hlines(0, 10000, 20000, linewidth=0.5, color=cc[1]) + # ax1.hlines(0, 20, 40, linewidth=0.5, color=cc[0]) + # ax2.hlines(0, 10000, 20000, linewidth=0.5, color=cc[1]) # share color cycles to prevent color re-use ax2._get_lines.prop_cycler = ax1._get_lines.prop_cycler diff --git a/lib/plotwav.py b/lib/plotwav.py index fc491cb..48e8e5f 100644 --- a/lib/plotwav.py +++ b/lib/plotwav.py @@ -5,7 +5,9 @@ from . import new_response, magnitude_x, convolve_each, monoize, count_channels import numpy as np -def plotfftsmooth(s, srate, ax=None, bw=1, tilt=None, size=8192, window=0, raw=False, **kwargs): + +def plotfftsmooth(s, srate, ax=None, bw=1, tilt=None, size=8192, + window=0, raw=False, **kwargs): sm = monoize(s) xs_raw = magnitude_x(srate, size) @@ -16,11 +18,13 @@ def plotfftsmooth(s, srate, ax=None, bw=1, tilt=None, size=8192, window=0, raw=F xs, ys = smoothfft(xs_raw, ys_raw, bw=bw) if ax: - if raw: ax.semilogx(xs_raw, ys_raw, **kwargs) + if raw: + ax.semilogx(xs_raw, ys_raw, **kwargs) ax.semilogx(xs, ys, **kwargs) return xs, ys + def plotwavinternal(sm, ss, srate, bw=1, size=8192, smoother=smoothfft2): xs_raw = magnitude_x(srate, size) ys_raw_m = averfft(sm, size=size) @@ -38,6 +42,7 @@ def plotwavinternal(sm, ss, srate, bw=1, size=8192, smoother=smoothfft2): return xs, ys_m, ys_s + def plotwav2(fn, bw=1, size=8192, fix=False, smoother=smoothfft2, **kwargs): s, srate = wav_read(fn) @@ -64,19 +69,22 @@ def plotwav2(fn, bw=1, size=8192, fix=False, sf = np.array((smf + ssf, smf - ssf)).T import ewave - with ewave.open(fno, 'w', sampling_rate=srate, nchannels=count_channels(sf)) as f: + with ewave.open(fno, 'w', sampling_rate=srate, + nchannels=count_channels(sf)) as f: f.write(sf) print('wrote '+fno) return xs, ys_m, ys_s + def pw2(fn, label=None, bw=1/6, **kwargs): fno = fn[:-4]+"-proc.wav" xs, ys_m, ys_s = plotwav2(fn, fix=True, bw=bw, **kwargs) xs, ys_m, ys_s = plotwav2(fno, fix=False, bw=bw, **kwargs) fig, ax = new_response(-18, 18) - ax.set_title('averaged magnitudes of normalized songs with tilt and smoothing') + ax.set_title( + 'averaged magnitudes of normalized songs with tilt and smoothing') label = label or fn ax.semilogx(xs, ys_m + 0, label=label+' (mid)') ax.semilogx(xs, ys_s + 9, label=label+' (side)') diff --git a/lib/smoothfft.py b/lib/smoothfft.py index 380b2a5..e670cf0 100644 --- a/lib/smoothfft.py +++ b/lib/smoothfft.py @@ -1,6 +1,7 @@ from . import xsp, lament 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.""" @@ -18,6 +19,7 @@ def smoothfft(xs, ys, bw=1, precision=512): 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.""" @@ -28,10 +30,11 @@ def smoothfft2(xs, ys, bw=1, precision=512, compensate=True): for i, x in enumerate(xs): # before optimizations: 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) + # 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) + _, temp = smoothfft2(xs, np.ones(len(xs)), + bw=bw, precision=precision, compensate=False) ys2 /= temp return xs2, ys2 diff --git a/lib/svf.py b/lib/svf.py index 43bd252..8a1e865 100644 --- a/lib/svf.py +++ b/lib/svf.py @@ -2,8 +2,10 @@ 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): + # via: + # http://nbviewer.ipython.org/urls/music-synthesizer-for-android.googlecode.com/git/lab/Second%20order%20sections%20in%20matrix%20form.ipynb # TODO: implement constant gain parameter g = unwarp(w0)*shelfA a1 = 1/(1 + g*(g + 1/Q)) @@ -17,52 +19,79 @@ def svf_core(w0, Q, m, shelfA=1, gain=1): 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)) +def LP2S(A, Q): + return (Q, [0, 0, 1], 1) + + +def BP2S(A, Q): + return (Q, [0, 1, 0], 1) + + +def HP2S(A, Q): + return (Q, [1, -1/Q, -1], 1) + + +# TODO: AP2S +# TODO: BP2aS +# TODO: BP2bS + + +def NO2S(A, Q): + return (Q, [1, -1/Q, 0], 1) + + +def PE2S(A, Q): + return (Q*A, [1, (A*A - 1)/A/Q, 0], 1) + + +def LS2S(A, Q): + return (Q, [1, (A - 1)/Q, A*A - 1], 1/np.sqrt(A)) + + +def HS2S(A, Q): + return (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)) + +def gen_filters_svf(cascade, srate): + return [ + svf_core(tau*f[0]/srate, *f[1], gain=10**(f[2]/20)) for f in cascade + ] -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 + 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]], + 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 + 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) diff --git a/lib/sweeps.py b/lib/sweeps.py index dbd2ce4..01a1671 100644 --- a/lib/sweeps.py +++ b/lib/sweeps.py @@ -2,18 +2,19 @@ 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)) + domain = np.log((tau * end)/(tau * begin)) + ys = amp*np.sin((tau * begin)/domain*(np.exp(xs*domain) - 1)) return ys + def tsp(N, m=0.5): """ OATSP(Optimized Aoshima's Time-Stretched Pulse) generator @@ -36,7 +37,7 @@ def tsp(N, m=0.5): if N < 0: raise Exception("The number of length must be the positive number") - NN = 2**np.floor(np.log2(N)) # nearest + NN = 2**np.floor(np.log2(N)) # nearest NN2 = NN//2 M = np.round(NN2*m) diff --git a/lib/util.py b/lib/util.py index 333f3d4..c9424a7 100644 --- a/lib/util.py +++ b/lib/util.py @@ -2,33 +2,73 @@ 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 -ceil2 = lambda x: np.power(2, np.ceil(np.log2(x))) -pad2 = lambda x: np.r_[x, np.zeros(ceil2(len(x)) - len(x))] -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 dummy(*args, **kwargs): + return None + + +def lament(*args, **kwargs): + return print(*args, file=sys.stderr, **kwargs) + + +def toLK(x): + return -0.691 + 10*np.log10(x) + + +def toQ(bw): + return isqrt2/bw + + +def toA(db): + return 10**(db/40) + + +def unwarp(w): + return np.tan(w/2) + + +def warp(w): + return np.arctan(w)*2 + + +def ceil2(x): + return np.power(2, np.ceil(np.log2(x))) + + +def pad2(x): + return np.r_[x, np.zeros(ceil2(len(x)) - len(x))] + + +def rfft(src, size): + return np.fft.rfft(src, size*2) + + +def magnitude(src, size): + return 10*np.log10(np.abs(rfft(src, size))**2)[0:size] + + # x axis for plotting above magnitude -magnitude_x = lambda srate, size: np.arange(0, srate/2, srate/2/size) +def magnitude_x(srate, size): + return np.arange(0, srate/2, srate/2/size) + + +def degrees_clamped(x): + return ((x*180/np.pi + 180) % 360) - 180 -degrees_clamped = lambda x: ((x*180/np.pi + 180) % 360) - 180 def xsp(precision=4096): - """create #precision log-spaced points from 20 Hz (inclusive) to 20480 Hz (exclusive)""" - xs = np.arange(0,precision)/precision + """ + create #precision log-spaced points from + 20 Hz (inclusive) to 20480 Hz (exclusive) + """ + 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: @@ -39,14 +79,18 @@ def blocks(a, step, size=None): break yield a[start:end] + def convolve_each(s, fir, mode='same', axis=0): - return np.apply_along_axis(lambda s: sig.fftconvolve(s, fir, mode), axis, s) + return np.apply_along_axis( + lambda s: sig.fftconvolve(s, fir, mode), axis, s) + def count_channels(s): if s.ndim < 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. @@ -56,6 +100,7 @@ def monoize(s): s = np.average(s, axis=1) return s + def div0(a, b): """division, whereby division by zero equals zero""" # http://stackoverflow.com/a/35696047 @@ -63,5 +108,5 @@ def div0(a, b): b = np.asanyarray(b) with np.errstate(divide='ignore', invalid='ignore'): c = np.true_divide(a, b) - c[~np.isfinite(c)] = 0 # -inf inf NaN + c[~np.isfinite(c)] = 0 # -inf inf NaN return c diff --git a/lib/wav.py b/lib/wav.py index b729f94..617d3fe 100644 --- a/lib/wav.py +++ b/lib/wav.py @@ -1,15 +1,18 @@ import numpy as np from .util import lament, count_channels + def wav_smart_read(fn): lament('wav_smart_read(): DEPRECATED; use wav_read instead.') - import scipy.io.wavfile as wav # don't use this, it fails to load good files + # don't use this, it fails to load good files. + import scipy.io.wavfile as wav 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('wav_smart_write(): DEPRECATED; use wav_write instead.') import scipy.io.wavfile as wav @@ -18,6 +21,7 @@ def wav_smart_write(fn, srate, s): si += np.clip(s*2**(bits - 1), -32768, 32767) wav.write(fn, srate, si) + def wav_read(fn): import ewave with ewave.open(fn) as f: @@ -30,6 +34,7 @@ def wav_read(fn): s = np.asfarray(s)/2**(bits - 1) return s, srate + def wav_write(fn, s, srate, dtype='h'): import ewave if dtype in ('b', 'h', 'i', 'l') and np.max(np.abs(s)) > 1: diff --git a/lib/windowing.py b/lib/windowing.py index acd9873..3847154 100644 --- a/lib/windowing.py +++ b/lib/windowing.py @@ -1,5 +1,6 @@ import numpy as np + def _deco_win(f): # gives scipy compatibility def deco(N, *args, sym=True, **kwargs): @@ -18,37 +19,46 @@ def _deco_win(f): return w return deco + def _gen_hamming(*a): L = len(a) - a += (0, 0, 0, 0, 0) # pad so orders definition doesn't error + a += (0, 0, 0, 0, 0) # pad so orders definition doesn't error orders = [ lambda fac: 0, lambda fac: a[0], - lambda fac: a[0] - a[1]*np.cos(fac), - lambda fac: a[0] - a[1]*np.cos(fac) + a[2]*np.cos(2*fac), - lambda fac: a[0] - a[1]*np.cos(fac) + a[2]*np.cos(2*fac) - a[3]*np.cos(3*fac), - lambda fac: a[0] - a[1]*np.cos(fac) + a[2]*np.cos(2*fac) - a[3]*np.cos(3*fac) + a[4]*np.cos(4*fac), + lambda fac: a[0] - a[1]*np.cos(1*fac), + lambda fac: a[0] - a[1]*np.cos(1*fac) + a[2]*np.cos(2*fac), + lambda fac: a[0] - a[1]*np.cos(1*fac) + a[2]*np.cos(2*fac) + - a[3]*np.cos(3*fac), + lambda fac: a[0] - a[1]*np.cos(1*fac) + a[2]*np.cos(2*fac) + - a[3]*np.cos(3*fac) + a[4]*np.cos(4*fac), ] f = orders[L] return lambda N: f(np.arange(0, N)*2*np.pi/(N - 1)) + def _normalize(*args): a = np.asfarray(args) return a/np.sum(a) -_h = lambda *args: _deco_win(_gen_hamming(*args)) + +def _h(*args): + return _deco_win(_gen_hamming(*args)) + + blackman_inexact = _h(0.42, 0.5, 0.08) blackman = _h(0.42659, 0.49656, 0.076849) blackman_nuttall = _h(0.3635819, 0.4891775, 0.1365995, 0.0106411) blackman_harris = _h(0.35875, 0.48829, 0.14128, 0.01168) nuttall = _h(0.355768, 0.487396, 0.144232, 0.012604) -flattop = _h(*_normalize(1, 1.93, 1.29, 0.388, 0.028)) # FTSRS -#flattop_weird = _h(*_normalize(1, 1.93, 1.29, 0.388, 0.032)) # ??? wtf -flattop_weird = _h(0.2156, 0.4160, 0.2781, 0.0836, 0.0069) # ??? scipy crap +flattop = _h(*_normalize(1, 1.93, 1.29, 0.388, 0.028)) # FTSRS +# flattop_weird = _h(*_normalize(1, 1.93, 1.29, 0.388, 0.032)) # ??? wtf +flattop_weird = _h(0.2156, 0.4160, 0.2781, 0.0836, 0.0069) # ??? scipy crap hann = _h(0.5, 0.5) hamming_inexact = _h(0.54, 0.46) hamming = _h(0.53836, 0.46164) + @_deco_win def triangular(N): if N % 2 == 0: @@ -56,6 +66,7 @@ def triangular(N): else: return 1 - np.abs(2*(np.arange(N) + 1)/(N + 1) - 1) + @_deco_win def parzen(N): odd = N % 2 @@ -71,17 +82,22 @@ def parzen(N): else: return np.r_[outer[::-1], center[::-1], center[1:], outer] + @_deco_win def cosine(N): return np.sin(np.pi*(np.arange(N) + 0.5)/N) + @_deco_win def welch(N): return 1 - (2*np.arange(N)/(N - 1) - 1)**2 + # TODO: rename or something @_deco_win def sinc(N): return np.sinc((np.arange(N) - (N - 1)/2)/2) -winmod = lambda f: lambda N: f(N + 2)[1:-1] + +def winmod(f): + return lambda N: f(N + 2)[1:-1]