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)

color='black', ls='solid', alpha=0.9)
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


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?