This jupyter notebook is part of a collection of notebooks on various topics of Digital Signal Processing. Please direct questions and suggestions to Sascha.Spors@uni-rostock.de.
Computing the output $y[k] = \mathcal{H} \{ x[k] \}$ of a linear time-invariant (LTI) system is of central importance in digital signal processing. This is often referred to as filtering of the input signal $x[k]$. We already have discussed the realization of non-recursive filters. This section focuses on the realization of recursive filters.
Linear difference equations with constant coefficients represent linear time-invariant (LTI) systems
\begin{equation} \sum_{n=0}^{N} a_n \; y[k-n] = \sum_{m=0}^{M} b_m \; x[k-m] \end{equation}where $y[k] = \mathcal{H} \{ x[k] \}$ denotes the response of the system to the input signal $x[k]$, $N$ the order, $a_n$ and $b_m$ constant coefficients, respectively. Above equation can be rearranged with respect to the output signal $y[k]$ by extracting the first element ($n=0$) of the left-hand sum
\begin{equation} y[k] = \frac{1}{a_0} \left( \sum_{m=0}^{M} b_m \; x[k-m] - \sum_{n=1}^{N} a_n \; y[k-n] \right) \end{equation}It is evident that the output signal $y[k]$ at time instant $k$ is given as a linear combination of past output samples $y[k-n]$ superimposed by a linear combination of the actual $x[k]$ and past $x[k-m]$ input samples. Hence, the actual output $y[k]$ is composed from the two contributions
The impulse response of the system is given as the response of the system to a Dirac impulse at the input $h[k] = \mathcal{H} \{ \delta[k] \}$. Using above result and the properties of the discrete Dirac impulse we get
\begin{equation} h[k] = \frac{1}{a_0} \left( b_k - \sum_{n=1}^{N} a_n \; h[k-n] \right) \end{equation}Due to the feedback, the impulse response will in general be of infinite length. The impulse response is termed as infinite impulse response (IIR) and the system as recursive system/filter.
Applying a $z$-transform to the left- and right-hand side of the difference equation and rearranging terms yields the transfer function $H(z)$ of the system
\begin{equation} H(z) = \frac{Y(z)}{X(z)} = \frac{\sum_{m=0}^{M} b_m \; z^{-m}}{\sum_{n=0}^{N} a_n \; z^{-n}} \end{equation}The transfer function is given as a rational function in $z$. The polynominals of the numerator and denominator can be expressed alternatively by their roots as
\begin{equation} H(z) = \frac{b_M}{a_N} \cdot \frac{\prod_{\mu=1}^{P} (z - z_{0\mu})^{m_\mu}}{\prod_{\nu=1}^{Q} (z - z_{\infty\nu})^{n_\nu}} \end{equation}where $z_{0\mu}$ and $z_{\infty\nu}$ denote the $\mu$-th zero and $\nu$-th pole of degree $m_\mu$ and $n_\nu$ of $H(z)$, respectively. The total number of zeros and poles is denoted by $P$ and $Q$. Due to the symmetries of the $z$-transform, the transfer function of a real-valued system $h[k] \in \mathbb{R}$ exhibits complex conjugate symmetry
\begin{equation} H(z) = H^*(z^*) \end{equation}Poles and zeros are either real valued or complex conjugate pairs for real-valued systems ($b_m\in\mathbb{R}$, $a_n\in\mathbb{R}$). For the poles of a causal and stable system $H(z)$ the following condition has to hold
\begin{equation} \max_{\nu} | z_{\infty\nu} | < 1 \end{equation}Hence, all poles have to be located inside the unit circle $|z| = 1$. Amongst others, this implies that $M \leq N$.
The following example shows the pole/zero diagram, the magnitude and phase response, and impulse response of a recursive filter with so-called Butterworth lowpass characteristic.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.markers import MarkerStyle
from matplotlib.patches import Circle
import scipy.signal as sig
N = 5 # order of recursive filter
L = 128 # number of computed samples
def zplane(z, p, title="Poles and Zeros"):
"Plots zero and pole locations in the complex z-plane"
ax = plt.gca()
ax.plot(np.real(z), np.imag(z), "bo", fillstyle="none", ms=10)
ax.plot(np.real(p), np.imag(p), "rx", fillstyle="none", ms=10)
unit_circle = Circle(
(0, 0), radius=1, fill=False, color="black", ls="solid", alpha=0.9
)
ax.add_patch(unit_circle)
ax.axvline(0, color="0.7")
ax.axhline(0, color="0.7")
plt.title(title)
plt.xlabel(r"Re{$z$}")
plt.ylabel(r"Im{$z$}")
plt.axis("equal")
plt.xlim((-2, 2))
plt.ylim((-2, 2))
plt.grid()
# compute coefficients of recursive filter
b, a = sig.butter(N, 0.2, "low")
# compute transfer function
Om, H = sig.freqz(b, a)
# compute impulse response
k = np.arange(L)
x = np.where(k == 0, 1.0, 0)
h = sig.lfilter(b, a, x)
# plot pole/zero-diagram
plt.figure(figsize=(5, 5))
zplane(np.roots(b), np.roots(a))
# plot magnitude response
plt.figure(figsize=(10, 3))
plt.plot(Om, 20 * np.log10(abs(H)))
plt.xlabel(r"$\Omega$")
plt.ylabel(r"$|H(e^{j \Omega})|$ in dB")
plt.grid()
plt.title("Magnitude response")
# plot phase response
plt.figure(figsize=(10, 3))
plt.plot(Om, np.unwrap(np.angle(H)))
plt.xlabel(r"$\Omega$")
plt.ylabel(r"$\varphi (\Omega)$ in rad")
plt.grid()
plt.title("Phase response")
# plot impulse response (magnitude)
plt.figure(figsize=(10, 3))
plt.stem(20 * np.log10(np.abs(np.squeeze(h))))
plt.xlabel(r"$k$")
plt.ylabel(r"$|h[k]|$ in dB")
plt.grid()
plt.title("Impulse response (magnitude)")
Text(0.5, 1.0, 'Impulse response (magnitude)')
Exercise
N
of the filter?Solution: It can be concluded from the last illustration, showing the magnitude of the impulse response $|h[k]|$ on a logarithmic scale, that the magnitude of the impulse response decays continuously for increasing $k$ but does not become zero at some point. This behavior continues with increasing $k$ as can be observed when increasing the number L
of computed samples in above example. The magnitude response $|H(e^{j \Omega})|$ of the filter decays faster with increasing order N
of the filter.
Copyright
This notebook is provided as Open Educational Resource. Feel free to use the notebook for your own purposes. The text is licensed under Creative Commons Attribution 4.0, the code of the IPython examples under the MIT license. Please attribute the work as follows: Sascha Spors, Digital Signal Processing - Lecture notes featuring computational examples.