Cascaded biquad filter implementationSource code in C is available on GitHub. The reason for cascade implementationIncreasing the Butterworth filter order makes faster roll-off around the cutoff frequency, while preserving flatness in the passband and stopband. 
Direct implementation of high order recursive filters would have coefficients that differ in many orders of magnitude, making a practical implementation difficult [Smith:CH20]. \[ \begin{equation} y_i = \sum_{k=0}^{N}{b_k x_{i-k}} - \sum_{k=1}^{N}{a_k y_{i-k}} \end{equation} \]Let's look at the 6-order low-pass Butterworth filter with sample rate 24000 Hz and cut-off requency 110 Hz.
Implementing this filter with single floating point precision requires changes in the algorithm1. The sum of input samples \(x_k\) multiplied by \(b_k\) will have order of magnitude much less than sum of output samples \(y_k\) multiplied by \(a_k\). Rounding errors make the filter unstable: 
Cascaded implementation has coefficients that are less different in orders of magnitude; it can be implemented with single precision:
Low-pass and high-pass filtersAnalog low-pass and high-pass filters have the same set of poles a zero-centered circle with radius \( \omega_c \). The frequency of analog prototype is pre-warped: \[ \begin{equation} \omega_c = \frac{2}{T} tan \biggl( \pi f_c T \biggr) \end{equation} \]For digital filter with cut-off frequency \(f_c\) = 880 and sample rate 8000 Hz, the analog prototype has frequency \(\omega_c\) = 5760.35445 
The poles of analog prototype are mapped to digital using bilinear transform: \[ \begin{equation} K = tan\frac{\omega T}{2} \end{equation} \] \[ \begin{equation} s = K \frac{z - 1}{z + 1} \end{equation} \] \[ \begin{equation} \omega_n = \frac{2n + 1}{4 \pi N} \end{equation} \] \[ \begin{equation} Q_n = \frac{1}{2 cos(\omega_n)} \end{equation} \]Applying bilinear transform and replacing \( K \) and \( K^2 \) with \[ \begin{equation} K = tan \frac{\omega T}{2} = \frac{sin \omega T}{1 + cos \omega T} \end{equation} \] \[ \begin{equation} K^2 = tan^2 \frac{\omega T}{2} = \frac{1 - cos \omega T}{1 + cos \omega T} \end{equation} \]gives the numerator and denominator terms used in the cookbook [RBJ]. The analog low-pass filter has zeros at \(s = \pm\infty\). Bi-linear transform maps them to \(z = -1\). The analog high-pass filter has zeros at \(s = 0 \). Bi-linear transform maps them to zeros at \(z = 1\). 
Notice that the middle stage is a second order Butterworth filter. Band-pass Butterworth filterLow-pass to band-pass transformation is performed in the frequency domain. The band-pass poles s are mapped to low-pass poles S by the following equation [BZ4]: \[ \begin{equation} S_i = \frac{1}{\gamma} \biggl( \frac{s_i^2 + \omega_0^2}{\omega_0 s_i} \biggr), \label{bptolp} \tag{2.1} \end{equation} \]where \( \omega_0 \) is the central frequency, geometric mean of \( \omega_l \) and \( \omega_h \), \[ \begin{equation} \omega_0 = \sqrt{\omega_l \omega_u} \end{equation} \]\( \gamma \) is the relative bandwidth, \[ \begin{equation} \gamma = \frac{\omega_u - \omega_l}{\omega_0} \end{equation} \]and \( T \) is the sampling period. Solving \ref{bptolp} maps every low-pass pole \(S_i\) to a pair of band-pass poles \(s_i\): \[ \begin{equation} s_i = \omega_0 \biggl( \frac{\gamma}{2} S_i \pm j \sqrt{1 - \biggl( \frac{\gamma S_i}{2} \biggr)^2} \biggr) \end{equation} \]
This implementation converts one of these poles to z-plane using bilinear transform, \[ \begin{equation} z_i = \frac{1 + s_i \frac{T}{2}}{1 - s_i \frac{T}{2}} \label{bilinear} \tag{2.2} \end{equation} \]and creates a second order section with poles at \( z \) and \( z^* \). \[ \begin{equation} \begin{aligned} a_{i, 0} & = & 1, \\ a_{i, 1} & = & -2 Re(z_i), \\ a_{i, 2} & = & |{z_i}|^2. \end{aligned} \end{equation} \]In the analog filter prototype, \(\omega_l\) and \(\omega_u\) are pre-warped frequencies of discrete filter: \[ \begin{equation} \omega = \frac{2}{T} tan \biggl( \pi f T \biggr) \end{equation} \]A pair of complex conjugate poles is not enough to define second-order polynom. For bandpass and bandstop filters, the result of bilinear transform has to be scaled to achieve unity gain in the passband [NR]. Every section of a band-pass filter has zeros at \( z = ±1 \). Band-stop filter zeros have to be mapped from continuous using bilinear transform \ref{bptolp}. Simply placing zeros at \( e^{\pm 2\pi f_0} \) would give wrong result. 
References[Smith] Smith, S.W. The Scientist and Engineer's Guide to Digital Signal Processing [RBJ] Robert Bristow-Johnson Cookbook formulae for audio EQ biquad filter coefficients [NR] Neil Robertson Design IIR Bandpass Filters [BZ4] Blinchikoff H.J., Zverev A.I.Filtering in the Time and Frequency Domains Chapter 4 © Ivan Gorinov, 2018 |