IIR filters

import numpy as np
import matplotlib.pyplot as plt

import sdr

%config InlineBackend.print_figure_kwargs = {"facecolor" : "w"}
# %matplotlib widget
colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]

Create an IIR filter

The user creates an IIR filter with the sdr.IIR class by specifying the feedforward coefficients \(b_i\) and feedback coefficients \(a_j\). Alternatively, an IIR filter may be created by specifying the zeros and poles in sdr.IIR.ZerosPoles.

Below is an IIR filter with one real zero and two complex-conjugate poles.

zero = 0.6
pole = 0.8 * np.exp(1j * np.pi / 8)
iir = sdr.IIR.ZerosPoles([zero], [pole, pole.conj()])
print(iir)
<sdr._iir_filter.IIR object at 0x0000028CBDA77348>
print(f"Feedforward taps: {iir.b_taps}")
print(f"Feedback taps: {iir.a_taps}")
Feedforward taps: [ 1.  -0.6]
Feedback taps: [ 1.         -1.47820725  0.64      ]

Examine the impulse response, \(h[n]\)

The impulse response of the IIR filter is computed and returned from the sdr.IIR.impulse_response() method. The impulse response \(h[n]\) is the output of the filter when the input is an impulse \(\delta[n]\).

h = iir.impulse_response(30)
print(h)
[ 1.00000000e+00  8.78207252e-01  6.58172329e-01  4.10862468e-01
  1.86109590e-01  1.21565653e-02 -1.01140214e-01 -1.57286400e-01
 -1.67772160e-01 -1.47338728e-01 -1.10422993e-01 -6.89312837e-02
 -3.12240078e-02 -2.03953322e-03  1.69685122e-02  2.63882791e-02
  2.81474977e-02  2.47193366e-02  1.85259041e-02  1.15647504e-02
  5.23851924e-03  3.42176895e-04 -2.84684395e-03 -4.42721858e-03
 -4.72236648e-03 -4.14721649e-03 -3.10813095e-03 -1.94024315e-03
 -8.78877688e-04 -5.74077567e-05]

The impulse response is conveniently plotted using the sdr.IIR.plot_impulse_response() method.

plt.figure(figsize=(10, 5))
iir.plot_impulse_response(30)
plt.show()
../../_images/c84556fa0164f9afcd95bbd7ca02235bfd5e60f016b8bcfcc9e2fe487e936e82.png

Examine the step response, \(s[n]\)

The step response of the IIR filter is computed and returned from the sdr.IIR.step_response() method. The step response \(s[n]\) is the output of the filter when the input is a unit step \(u[n]\).

s = iir.step_response(30)
print(s)
[1.         1.87820725 2.53637958 2.94724205 3.13335164 3.1455082
 3.04436799 2.88708159 2.71930943 2.5719707  2.46154771 2.39261642
 2.36139242 2.35935288 2.3763214  2.40270968 2.43085717 2.45557651
 2.47410241 2.48566716 2.49090568 2.49124786 2.48840102 2.4839738
 2.47925143 2.47510421 2.47199608 2.47005584 2.46917696 2.46911955]
plt.figure(figsize=(10, 5))
iir.plot_step_response(30)
plt.show()
../../_images/6aceef0cad52d4901b5740ef3d5054e27520c368b1f5760a16b7c17d9aadbb36.png

Examine the zeros and poles

Zeros are \(z\) values that set the numerator of \(H(z)\) to zero.

print(iir.zeros)
[0.6]

Poles are \(z\) values that set the denominator of \(H(z)\) to zero. The poles define the stability of the IIR filter.

print(iir.poles)
[0.73910363+0.30614675j 0.73910363-0.30614675j]

The zeros and poles are conveniently plotted in the complex plane using the sdr.IIR.plot_zeros_poles() method.

plt.figure(figsize=(5, 5))
iir.plot_zeros_poles()
plt.show()
../../_images/9d4cf4b665dc256cd8b7987d2aad7fa0069fb4ce452bfdd807b8566d4acdc6c3.png

Examine the frequency response, \(H(\omega)\)

The frequency response is the transfer function \(H(z)\) evaluated at the complex exponential \(e^{j \omega}\), where \(\omega = 2 \pi f / f_s\).

The frequency response is conveniently plotted using the sdr.IIR.plot_frequency_response() method.

plt.figure(figsize=(10, 5))
iir.plot_frequency_response()
plt.show()
../../_images/dadfc6850f933a564e950a50f97e28112185fb682480753d522a70c083c22a4a.png

The frequency response may also be plotted with a logarithmic frequency scale using the sdr.IIR.plot_frequency_response_log() method.

plt.figure(figsize=(10, 5))
iir.plot_frequency_response_log()
plt.show()
../../_images/0f8238b2f83a60d4653151ddcf20739fd0064fcf1066bf4d6edefb30e7f69492.png

Examine the group delay, \(\tau_g(\omega)\)

The group delay \(\tau_g(\omega)\) is the time shift of the envelope of a signal passed through the filter as a function of its frequency \(\omega\).

The group delay is conveniently plotted using the sdr.IIR.plot_group_delay() method.

plt.figure(figsize=(10, 5))
iir.plot_group_delay()
plt.show()
../../_images/1b8c5b867d1a307f9d5aeba29e5d6f561530c90690cec0dceafc50033c9b7b29.png

Fully analyze an IIR filter

The user can easily analyze the perform of a given IIR filter using the sdr.IIR.plot_all() method.

Here is an IIR filter with one real zero and 8 complex poles.

zeros = np.array([0.8])
poles = 0.6 * np.exp(1j * np.linspace(np.pi / 8, np.pi / 4, 4, endpoint=False))
poles = np.concatenate((poles, poles.conj()))
iir = sdr.IIR.ZerosPoles(zeros, poles)
print(iir)
<sdr._iir_filter.IIR object at 0x0000028CC010FD48>
plt.figure(figsize=(10, 8))
iir.plot_all(N_time=30)
plt.show()
../../_images/87a637c0a30af5b5070b55c8d066392c13d05bbe9ef86fc8e89e96d8439687fa.png

Poles and digital filter stability

Reference:

  • R. Lyons, Understanding Digital Signal Processing 3rd Edition, Section 6.3.1.

When the pole is real and inside the unit circle, the impulse response \(h[n]\) is an exponential decay.

zeros = []
poles = [0.8]
iir = sdr.IIR.ZerosPoles(zeros, poles)

plt.figure(figsize=(10, 8))
iir.plot_all(N_time=30)
plt.show()
../../_images/08b269852db7cd05ce5f61123bac874e90d2fc6978c631dacebf7aa3d7b655b6.png

When the poles are complex conjugates and inside the unit circle, the impulse response \(h[n]\) is a decaying sinusoid.

zeros = []
pole = 0.8 * np.exp(1j * np.pi / 8)
poles = [pole, pole.conj()]
iir = sdr.IIR.ZerosPoles(zeros, poles)

plt.figure(figsize=(10, 8))
iir.plot_all(N_time=30)
plt.show()
../../_images/c5cea2bf91816669ac6892621219929d50d4ef49678ae44d4e10202d8ade3774.png

When the pole is real and on the unit circle, the impulse response \(h[n]\) is constant. This filter is an integrator.

zeros = []
poles = [1]
iir = sdr.IIR.ZerosPoles(zeros, poles)

plt.figure(figsize=(10, 8))
iir.plot_all(N_time=30)
plt.show()
C:\Users\matth\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.7_qbz5n2kfra8p0\LocalCache\local-packages\Python37\site-packages\scipy\signal\filter_design.py:476: RuntimeWarning: divide by zero encountered in true_divide
  npp_polyval(zm1, a, tensor=False))
C:\Users\matth\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.7_qbz5n2kfra8p0\LocalCache\local-packages\Python37\site-packages\scipy\signal\filter_design.py:476: RuntimeWarning: invalid value encountered in true_divide
  npp_polyval(zm1, a, tensor=False))
../../_images/00de1f6f9f23e126ff3e7d4ee5815652d6a72df243466bed8aa6ae7a9b8980bb.png

When the poles are complex conjugates and on the unit circle, the impulse response \(h[n]\) is a sinusoid.

zeros = []
pole = 1 * np.exp(1j * np.pi / 8)
poles = [pole, pole.conj()]
iir = sdr.IIR.ZerosPoles(zeros, poles)

plt.figure(figsize=(10, 8))
iir.plot_all(N_time=30)
plt.show()
../../_images/4eb27c2b3e331b18498afc7288dfff3a6b5862819459ba6b217726b84542828b.png

When the pole is real and outside the unit circle, the impulse response \(h[n]\) is an exponential. This filter is unstable.

zeros = []
poles = [1.2]
iir = sdr.IIR.ZerosPoles(zeros, poles)

plt.figure(figsize=(10, 8))
iir.plot_all(N_time=30)
plt.show()
../../_images/bd0a47003a9a3a9ac5591f3b7459779cc79b7a5e9eb086fe9506bef81e4fe86a.png

When the poles are complex conjugates and outside the unit circle, the impulse response \(h[n]\) is an exponentially-increasing sinusoid. This filter is unstable.

zeros = []
pole = 1.2 * np.exp(1j * np.pi / 8)
poles = [pole, pole.conj()]
iir = sdr.IIR.ZerosPoles(zeros, poles)

plt.figure(figsize=(10, 8))
iir.plot_all(N_time=30)
plt.show()
../../_images/0ac2de37f952d333bb100f1fb1c2d83af84365ab74a9a696e9e0889b9ea24a47.png

Last update: Jul 09, 2023