Add plot_iir_filter_analysis() based on LTF-IIR-allgemein.ipynbcode from...

Add plot_iir_filter_analysis() based on LTF-IIR-allgemein.ipynbcode from https://github.com/chipmuenk .

applications/lofar2/model/pfb_os/dsp.py

@@ -22,11 +22,15 @@
 # Author: Eric Kooistra
 # Purpose: Utilities and functions for DSP
 # Description:
-
+#
+# References:
+# [1] dsp_study_erko.txt
+# [2] https://github.com/chipmuenk/dsp/blob/main/notebooks/02_LTF/LTF-IIR-allgemein.ipynb
 
 import numpy as np
+from scipy import signal
 import matplotlib.pyplot as plt
-
+import dsp_fpga_lib  # from [2]
 
 c_interpolate = 10
 c_atol = 1e-15
@@ -335,6 +339,58 @@ def plot_time_response(h, name='', markers=False):
     plt.grid(True)
 
 
+def plot_iir_filter_analysis(b, a, fs=1, Ntime=100, step=False):
+    """Plot iir filter analysis results.
+
+    Input:
+    . b, a: IIR filter coefficients in same format as for scipy.signal.freqz
+    . fs: sample frequency
+    . Ntime: number of timesamples for impulse response
+    . step: False for impulse response, True for step response
+    """
+
+    # Plot poles / zeros diagram in z-plane
+    fig1, ax1 = plt.subplots(1)
+    z, p, k = dsp_fpga_lib.zplane(b, a, plt_ax=ax1)  # uses np.roots(a), np.roots(b)
+    print('Zeros, poles and gain from b, a coefficients:')
+    if len(z) > 0:
+        print('. zeros:')
+        for zero in z:
+            print('  z = %s' % str(zero))
+    if len(p) > 0:
+        print('. poles:')
+        for pole in p:
+            print('  p = %s' % str(pole))
+    print('. gain: k = %.3f' % k)
+
+    # Derive b, a coefficients back from z, p, k
+    print('Coefficients back from z, p, k:')
+    print('  b = %s' % str(np.poly(z)))
+    print('  a = %s' % str(np.poly(p) / k))
+
+    # Plot transfer function H(f), is H(z) for z = exp(j w), so along the unit circle
+    # . 0 Hz at 1 + 0j, fs / 4 at 0 + 1j, fNyquist = fs / 2 at -1 + 0j
+    # . use whole=False to have only positve frequencies, so 0 to fs / 2
+    # . use Nfreq frequency points
+    fig2, ax2 = plt.subplots(1)
+    Nfreq = 1024
+    [f, HF] = signal.freqz(b, a, Nfreq, whole=False, fs=fs)
+    ax2.plot(f, pow_db(HF))
+    ax2.set_xlabel('frequency [fs = %f]' % fs)
+    ax2.set_ylabel('HF power [dB]')
+
+    # Plot impulse response
+    Ts = 1 / fs
+    fig3, ax3 = plt.subplots(1)
+    step = True  # step response (makes impz use np.cumsum(h))
+    step = False
+    [h, t] = dsp_fpga_lib.impz(b, a, FS=fs, N=Ntime, step=step)  # uses signal.lfilter()
+    (ml, sl, bl) = ax3.stem(t, h, linefmt='b-', markerfmt='ro', basefmt='k')
+    ax3.set_xlabel('time [Ts = %f]' % Ts)
+    ax3.set_ylabel('h[n]')
+    return z, p, k
+
+
 def plot_spectra(f, HF, fs=1.0, fLim=None, dbLim=None):
     """Plot spectra for power, magnitude, phase, real, imag