Cascaded biquad filter implementation

Source code in C is available on GitHub.

The reason for cascade implementation

Increasing 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.

08.43345791e-12  1
15.06007475e-11 -5.88873409
21.26501869e-10 14.4498436
41.26501869e-10 13.92373555
55.06007475e-11 -5.46772375
68.43345791e-12  0.89469577

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 filters

Analog 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

Analog low-pass filter poles Digital low-pass filter poles

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

Low-pass filter frequency response High-pass filter frequency response

Notice that the middle stage is a second order Butterworth filter.

Band-pass Butterworth filter

Low-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} \] Analog band-pass filter poles Analog band-pass filter poles

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.

Band-pass filter frequency response Band-stop filter frequency response

Test program source


[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