7.3. Cascaded Structures

The realization of recursive filters with a high order may be subject to numerical issues. For instance, when the coefficients span a wide amplitude range, their quantization may require a small quantization step or may impose a large relative error for small coefficients. The basic concept of cascaded structures is to decompose a high order filter into a cascade of lower order filters, typically first and second order recursive filters.

7.3.1. Decomposition into Second-Order Sections

The rational transfer function \(H(z)\) of a linear time-invariant (LTI) recursive system can be expressed by its zeros and poles as

\[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}}\]

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\).

The poles and zeros of a real-valued filter \(h[k] \in \mathbb{R}\) are either single real valued or conjugate complex pairs. This motivates to split the transfer function into

  • first order filters constructed from a single pole and zero
  • second order filters constructed from a pair of conjugated complex poles and zeros

Decomposing the transfer function into these two types by grouping the poles and zeros into single poles/zeros and conjugate complex pairs of poles/zeros results in

\[ H(z) = K \cdot \prod_{\eta=1}^{S_1} \frac{(z - z_{0\eta})}{(z - z_{\infty\eta})} \cdot \prod_{\eta=1}^{S_2} \frac{(z - z_{0\eta}) (z - z_{0\eta}^*)} {(z - z_{\infty\eta})(z - z_{\infty\eta}^*)}\]

where \(K\) denotes a constant and \(S_1 + 2 S_2 = N\) with \(N\) denoting the order of the system. The cascade of two systems results in a multiplication of their transfer functions. Above decomposition represents a cascade of first- and second-order recursive systems. The former can be treated as a special case of second-order recursive systems. The decomposition is therefore known as decomposition into second-order sections (SOSs) or biquad filters. Using a cascade of SOSs the transfer function of the recursive system can be rewritten as

\[H(z) = \prod_{\mu=1}^{S} \frac{b_{0, \mu} + b_{1, \mu} \, z^{-1} + b_{2, \mu} \, z^{-2}}{1 + a_{1, \mu} \, z^{-1} + a_{2, \mu} \, z^{-2}}\]

where \(S = \lceil \frac{N}{2} \rceil\) denotes the total number of SOSs. These results state that any real valued system of order \(N > 2\) can be decomposed into SOSs. This has a number of benefits

  • quantization effects can be reduced by sensible grouping of poles/zeros, e.g. such that the spanned amplitude range of the filter coefficients is limited
  • A SOS may be extended by a gain factor to further reduce quantization effects by normalization of the coefficients
  • efficient and numerically stable SOSs serve as generic building blocks for higher-order recursive filters

7.3.2. Example

The following example illustrates the decomposition of a higher-order recursive Butterworth lowpass filter into a cascade of second-order sections.

In [1]:
%matplotlib inline
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 = 9  # order of recursive filter

def zplane(z, p):

    ax = plt.gca()
    plt.hold(True)

    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.axis('equal')
    plt.xlim((-2, 2))
    plt.ylim((-2, 2))
    plt.grid()
    plt.xlabel(r'Re{$z$}')
    plt.ylabel(r'Im{$z$}')

    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)

    plt.hold(False)


# design filter
b, a = sig.butter(N, 0.2)
# decomposition into SOS
sos = sig.tf2sos(b, a, pairing='nearest')


# print filter coefficients
print('Coefficients of the recursive part')
print(['%1.2f'%ai for ai in a])
print('\n')
print('Coefficients of the recursive part of the SOSs')
print('Section \t a1 \t\t a2')
for n in range(sos.shape[0]):
    print('%d \t\t %1.5f \t %1.5f'%(n, sos[n, 4], sos[n, 5]))

# plot poles and zeros
plt.figure(figsize=(5,5))
zplane(np.roots(b), np.roots(a))
plt.title('Overall')

for n in range(sos.shape[0]):
    if not n%3:
        plt.figure(figsize=(10, 3.5))
    plt.subplot(131+n%3)
    zplane(np.roots(sos[n, 0:3]), np.roots(sos[n, 3:6]))
    plt.title('Section %d'%n)
    plt.tight_layout()

# compute and plot frequency response of sections
plt.figure(figsize=(10,5))
for n in range(sos.shape[0]):
    Om, H = sig.freqz(sos[n, 0:3], sos[n, 3:6])
    plt.plot(Om, 20*np.log10(np.abs(H)), label=r'Section $n=%d$'%n)
    plt.hold(True)
plt.hold(False)
plt.xlabel(r'$\Omega$')
plt.ylabel(r'$|H_n(e^{j \Omega})|$ in dB')
plt.legend()
plt.grid();
Coefficients of the recursive part
['1.00', '-5.39', '13.38', '-19.96', '19.62', '-13.14', '5.97', '-1.78', '0.31', '-0.02']


Coefficients of the recursive part of the SOSs
Section          a1              a2
0                -0.50953        0.00000
1                -1.04232        0.28838
2                -1.11568        0.37905
3                -1.25052        0.54572
4                -1.46818        0.81477
../_images/recursive_filters_cascaded_structures_3_1.png
../_images/recursive_filters_cascaded_structures_3_2.png
../_images/recursive_filters_cascaded_structures_3_3.png
../_images/recursive_filters_cascaded_structures_3_4.png

Exercise

  • What amplitude range is spanned by the filter coefficients?
  • What amplitude range is spanned by the SOS coefficients?
  • Change the pole/zero grouping strategy from pairing='nearest' to pairing='keep_odd'. What changes?
  • Increase the order N of the filter. What changes?