update 17

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

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]

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)
 

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)

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):

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]