diff --git a/lib/__init__.py b/lib/__init__.py index 6f406c0..2f89c84 100644 --- a/lib/__init__.py +++ b/lib/__init__.py @@ -10,6 +10,8 @@ from .planes import * from .fft import * from .bs import * from .cepstrum import * +from .windowing import * +from .piir import * import numpy as np from matplotlib.pylab import show diff --git a/lib/piir.py b/lib/piir.py new file mode 100644 index 0000000..fc4d1d9 --- /dev/null +++ b/lib/piir.py @@ -0,0 +1,120 @@ +import numpy as np + +# i'll be dumping coefficients here until i port the generator. +# coefficients via https://gist.github.com/notwa/3be345efb6c97d757398 +# which is a port of http://ldesoras.free.fr/prod.html#src_hiir + +halfband_c = {} + +halfband_c['16,0.1'] = [ + # reject: roughly -155 dB + 0.006185967461045014, + 0.024499027624721819, + 0.054230780876613788, + 0.094283481125726432, + 0.143280861566087270, + 0.199699579426327684, + 0.262004358403954640, + 0.328772348316831664, + 0.398796973552973666, + 0.471167216679969414, + 0.545323651071132232, + 0.621096845120503893, + 0.698736833646440347, + 0.778944517099529166, + 0.862917812650502936, + 0.952428157718303137, +] + +halfband_c['8,0.01'] = [ + # reject: -69 dB + 0.077115079832416222, + 0.265968526521094595, + 0.482070625061047198, + 0.665104153263495701, + 0.796820471331579738, + 0.884101508550615867, + 0.941251427774047134, + 0.982005414188607539, +] + +halfband_c['olli'] = [ + # via http://yehar.com/blog/?p=368 + # "Transition bandwidth is 0.002 times the width of passband, + # stopband is attenuated down to -44 dB + # and passband ripple is 0.0002 dB." + # roughly equivalent to ./halfband 8 0.0009074 + # reject: -44 dB + 0.4021921162426**2, + 0.6923878000000**2, + 0.8561710882420**2, + 0.9360654322959**2, + 0.9722909545651**2, + 0.9882295226860**2, + 0.9952884791278**2, + 0.9987488452737**2, +] + +class Halfband: + def __init__(self, c='olli'): + self.x = np.zeros(4) + + if isinstance(c, str): + c = halfband_c[c] + self.c = np.asfarray(c) + + self.len = len(c)//2 + self.a2 = np.zeros(self.len) + self.b2 = np.zeros(self.len) + self.a1 = np.zeros(self.len) + self.b1 = np.zeros(self.len) + + def process_halfband(self, xs): + real, imag = self.process_all(xs, mode='filter') + return (real + imag)*0.5 + + def process_power(self, xs): + real, imag = self.process_all(xs, mode='hilbert') + return np.sqrt(real**2 + imag**2) + + def process_all(self, xs, mode='hilbert'): + self.__init__() + real = np.zeros(len(xs)) + imag = np.zeros(len(xs)) + for i, x in enumerate(xs): + real[i], imag[i] = self.process(x, mode=mode) + return real, imag + + def process(self, x0, mode='hilbert'): + a = np.zeros(self.len) + b = np.zeros(self.len) + self.x[0] = x0 + + sign = 1 + if mode == 'hilbert': + #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] + sign = -1 + + in2 = self.x[2] + in0 = self.x[0] + for i in range(self.len): + in0 = a[i] = self.c[i*2+0]*(in0 + sign*self.a2[i]) - sign*in2 + in2 = self.a2[i] + + in2 = self.x[3] + in0 = self.x[1] + for i in range(self.len): + in0 = b[i] = self.c[i*2+1]*(in0 + sign*self.b2[i]) - sign*in2 + in2 = self.b2[i] + + self.x[3] = self.x[2] + self.x[2] = self.x[1] + self.x[1] = self.x[0] + + self.a2, self.a1 = self.a1, a + self.b2, self.b1 = self.b1, b + + return a[-1], b[-1] diff --git a/lib/plotwav.py b/lib/plotwav.py index 23b2b57..fc491cb 100644 --- a/lib/plotwav.py +++ b/lib/plotwav.py @@ -1,14 +1,14 @@ # this is a bunch of crap that should really be reduced to one or two functions from . import wav_read, normalize, averfft, tilter2, smoothfft2, firize -from . import new_response, convolve_each, monoize, count_channels +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): sm = monoize(s) - xs_raw = np.arange(0, srate/2, srate/2/size) + xs_raw = magnitude_x(srate, size) ys_raw = averfft(sm, size=size, mode=window) ys_raw -= tilter2(xs_raw, tilt) @@ -22,7 +22,7 @@ def plotfftsmooth(s, srate, ax=None, bw=1, tilt=None, size=8192, window=0, raw=F return xs, ys def plotwavinternal(sm, ss, srate, bw=1, size=8192, smoother=smoothfft2): - xs_raw = np.arange(0, srate/2, srate/2/size) + xs_raw = magnitude_x(srate, size) ys_raw_m = averfft(sm, size=size) ys_raw_s = averfft(ss, size=size) diff --git a/lib/smoothfft.py b/lib/smoothfft.py index a6a031e..380b2a5 100644 --- a/lib/smoothfft.py +++ b/lib/smoothfft.py @@ -21,6 +21,7 @@ def smoothfft(xs, ys, bw=1, precision=512): def smoothfft2(xs, ys, bw=1, precision=512, compensate=True): """performs log-lin smoothing on magnitude data, generally from the output of averfft.""" + # this is probably implementable with FFTs now that i think about it xs2 = xsp(precision) ys2 = np.zeros(precision) log2_xs2 = np.log2(xs2) diff --git a/lib/util.py b/lib/util.py index 2c71c80..7feb21c 100644 --- a/lib/util.py +++ b/lib/util.py @@ -19,6 +19,8 @@ pad2 = lambda x: np.hstack((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] +# x axis for plotting above magnitude +magnitude_x = lambda srate, size: np.arange(0, srate/2, srate/2/size) degrees_clamped = lambda x: ((x*180/np.pi + 180) % 360) - 180 @@ -39,7 +41,7 @@ def blocks(a, step, size=None): break yield a[start:end] -def convolve_each(s, fir, mode=None, axis=0): +def convolve_each(s, fir, mode='same', axis=0): return np.apply_along_axis(lambda s: sig.fftconvolve(s, fir, mode), axis, s) def count_channels(s): diff --git a/lib/windowing.py b/lib/windowing.py new file mode 100644 index 0000000..4ee2926 --- /dev/null +++ b/lib/windowing.py @@ -0,0 +1,41 @@ +from . import tau +import numpy as np + +def gen_hamm(a0, a1=0, a2=0, a3=0, a4=0): + form = lambda x: a0 - a1*np.cos(x) + a2*np.cos(2*x) - a3*np.cos(3*x) + a4*np.cos(4*x) + dc = form(np.pi) + def f(N): + a = tau*np.arange(N)/(N - 1) + return form(a)/dc + return f + +# TODO: rename or something +sinc = lambda N: np.sinc((np.arange(N) - (N - 1)/2)/2) + +triangular = lambda N: 1 - np.abs(2*np.arange(N)/(N - 1) - 1) +welch = lambda N: 1 - (2*np.arange(N)/(N - 1) - 1)**2 +cosine = lambda N: np.sin(np.pi*np.arange(N)/(N - 1)) + +hann = gen_hamm(0.5, 0.5) +hamming = gen_hamm(0.53836, 0.46164) +blackman = gen_hamm(0.42659, 0.49656, 0.076849) +blackman_harris = gen_hamm(0.35875, 0.48829, 0.14128, 0.01168) +blackman_nutall = gen_hamm(0.3635819, 0.4891775, 0.1365995, 0.0106411) +nutall = gen_hamm(0.355768, 0.487396, 0.144232, 0.012604) +# specifically the Stanford Research flat-top window (FTSRS) +flattop = gen_hamm(1, 1.93, 1.29, 0.388, 0.028) + +def parzen(N): + # TODO: verify this is accurate. the code here is kinda sketchy + n = np.arange(N/2) + center = 1 - 6*(2*n/N)**2*(1 - 2*n/N) + outer = 2*(1 - 2*n/N)**3 + center = center[:len(center)//2] + outer = outer[len(outer)//2:] + if N % 2 == 0: + return np.hstack((outer[::-1], center[::-1], center, outer)) + else: + return np.hstack((outer[::-1], center[::-1], center[1:], outer)) + +# some windows reach 0 at either side, so this can avoid that +winmod = lambda f: lambda N: f(N + 2)[1:-1]