Add iir_biquad_butter_coefficients(), for biquads used in TBB.

applications/lofar2/model/pfb_os/dsp.py

@@ -117,7 +117,7 @@ def impulse_at_zero_crossing(x):
 
 
 ###############################################################################
-# Filter design
+# FIR Filter design
 ###############################################################################
 
 def nof_taps_kaiser_window(fs, fpass, fstop, atten_db):
@@ -157,7 +157,7 @@ def ideal_low_pass_filter(Npoints, Npass, bandEdgeGain=1.0):
     # Magnitude frequency reponse
     HF = np.zeros([Npoints])
     HF[0] = bandEdgeGain
-    HF[1 : Npass - 1] = 1
+    HF[1: Npass - 1] = 1
     HF[Npass - 1] = bandEdgeGain
     # Zero center HF to make it even
     HF = np.roll(HF, -(Npass // 2))
@@ -222,6 +222,39 @@ def fourier_interpolate(HFfilter, Ncoefs):
     return hInterpolated.real
 
 
+###############################################################################
+# IIR Filter design
+###############################################################################
+
+def iir_biquad_butter_coefficients(f0, BW, band='bandstop', fs=1):
+    """Calculate b0, b1, b2, a0, a1, a2 for biquad IIR filter section.
+
+    Bilinear transform of Butterworth filters. From
+    https://webaudio.github.io/Audio-EQ-Cookbook/Audio-EQ-Cookbook.txt
+
+    For single biquad bandpass alternatively use signal.iirpeak()
+    For single biquad bandstop alternatively use signal.iirnotch()
+    """
+    w0 = 2 * np.pi * f0 / fs
+    Q = f0 / BW
+    si = np.sin(w0)
+    cs = np.cos(w0)
+    alpha = np.sin(w0) / (2 * Q)
+    a1 = -2 * cs
+    if band == 'lowpass':
+        b0 = 1 - cs
+        b = [b0 / 2, b0, b0 / 2]
+    elif band == 'highpass':
+        b0 = 1 + cs
+        b = [b0 / 2, -b0, b0 / 2]
+    elif band == 'bandstop':
+        b = [1, a1, 1]
+    elif band == 'bandpass':
+        b = [si / 2, 0, -si / 2]
+    a = [1 + alpha, a1, 1 - alpha]
+    return b, a
+
+
 ###############################################################################
 # Hilbert transform filter
 ###############################################################################
@@ -475,7 +508,7 @@ def plot_spectra(f, HF, fs=1.0, fLim=None, dbLim=None):
     plt.grid(True)
 
 
-def plot_power_spectrum(f, HF, fmt='r', fs=1.0, fLim=None, dbLim=None, showRoll=False):
+def plot_power_spectrum(f, HF, fmt='r', fs=1.0, fLim=None, dbLim=None):
     """Plot power spectrum
 
     Input:
@@ -485,7 +518,7 @@ def plot_power_spectrum(f, HF, fmt='r', fs=1.0, fLim=None, dbLim=None, showRoll=
     . fs: sample frequency in Hz, scale f by fs, fs >= 1
     """
     if fs > 1:
-        flabel = 'Frequency [fs / %d]' % fs
+        flabel = 'Frequency'
     else:
         flabel = 'Frequency [fs]'