update 18

lib/__init__.py

@@ -51,19 +51,31 @@ def c_render(cascade, precision=4096):
     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_render(xs, cascade, mode='phase')
+    else:
+        return c_render(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 not phase:
+
+    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:
-        fmt = 'arctan2(imag({0}), real({0}))' # gross
+        raise Exception("c_render(): unknown mode: {}".format(mode))
+
     ys = np.zeros(len(xs))
     for f in cascade:
-        w0, ba, gain = f
+        f, ba, gain = f
         b, a = ba
         if len(b) == 3 and len(a) == 3:
             eq = fmt.format(eq2)
@@ -75,9 +87,24 @@ def c_render2(xs, cascade, phase=False):
             a1, a0 = a
         else:
             raise Exception("incompatible cascade; consider using c_render instead")
-        s = xs/w0*j
+
-        ys += ne.evaluate(eq)
+        if mode == 'group delay':
-    if phase:
+            # approximate derivative of phase by slope of tangent line
+            step = 2**-8
+            fa = f - step
+            fb = f + 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/f*tau)
+        else:
+            s = xs/f*j
+            ys += ne.evaluate(eq)
+    if mode == 'phase':
         ys = degrees_clamped(ys)
     return ys
 

lib/bq.py

@@ -6,6 +6,7 @@ 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)
 
 LP1 = lambda A, Q: ((0,1),(1,1))

lib/plot.py

@@ -1,8 +1,6 @@
 import matplotlib.pyplot as plt
 from matplotlib import ticker
 
-# TODO: remove set_size_inches calls, move them inline as necessary
-
 def response_setup(ax, ymin=-24, ymax=24, yL=ticker.AutoMinorLocator(3)):
     ax.set_xlim(20, 20000)
     ax.set_ylim(ymin, ymax)
@@ -23,7 +21,6 @@ def phase_response_setup(ax, div=12, yL=ticker.AutoMinorLocator(2)):
 
 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
@@ -32,14 +29,12 @@ def new_response(*args, **kwargs):
     fig = plt.figure()
     ax = fig.gca()
     response_setup(ax, *args, **kwargs)
-    fig.set_size_inches(10, 6)
     return fig, ax
 
 def new_phase_response(*args, **kwargs):
     fig = plt.figure()
     ax = fig.gca()
     phase_response_setup(ax, *args, **kwargs)
-    fig.set_size_inches(10, 6)
     return fig, ax
 
 def new_bode(magnitude_offset=0):
@@ -50,8 +45,22 @@ def new_bode(magnitude_offset=0):
     ymax = 24 + magnitude_offset
     response_setup(ax1, ymin, ymax)
     phase_response_setup(ax2)
+
+    cc = plt.style.library['ggplot']['axes.color_cycle']
+    ax1.set_ylabel(ax1.get_ylabel(), color=cc[0])
+    ax2.set_ylabel(ax2.get_ylabel(), color=cc[1])
+    for tl in ax1.get_yticklabels():
+        tl.set_color(cc[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 will use the initial color, so tell ax2 to iterate to the second
+    next(ax2._get_lines.color_cycle)
+
     # ax1 and ax2 should have identical grids;
     # disable ax2's so it doesn't overlap ax1's plot lines.
     ax2.grid(False, which='both')
-    fig.set_size_inches(10, 6)
     return fig, ax1, ax2