Source code in C is available on GitHub.

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

$$$y_i = \sum_{k=0}^{N}{b_k x_{i-k}} - \sum_{k=1}^{N}{a_k y_{i-k}}$$$

Let's look at the 6-order low-pass Butterworth filter with sample rate 24000 Hz and cut-off requency 110 Hz.

ba
08.43345791e-12  1
15.06007475e-11 -5.88873409
21.26501869e-10 14.4498436
31.68669158e-10-18.91181708
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:

Stageb0b1b2a0a1a2
00.000201710.000403410.000201711-1.945072770.9458796
10.000203180.000406360.000203181-1.959279030.96009175
20.000205780.000411560.000205781-1.984382280.98520541

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:

$$$\omega_c = \frac{2}{T} tan \biggl( \pi f_c T \biggr)$$$

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:

$$$K = tan\frac{\omega T}{2}$$$ $$$s = K \frac{z - 1}{z + 1}$$$ $$$\omega_n = \frac{2n + 1}{4 \pi N}$$$ $$$Q_n = \frac{1}{2 cos(\omega_n)}$$$

Applying bilinear transform and replacing $$K$$ and $$K^2$$ with

$$$K = tan \frac{\omega T}{2} = \frac{sin \omega T}{1 + cos \omega T}$$$ $$$K^2 = tan^2 \frac{\omega T}{2} = \frac{1 - cos \omega T}{1 + cos \omega T}$$$

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 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]:

$$$S_i = \frac{1}{\gamma} \biggl( \frac{s_i^2 + \omega_0^2}{\omega_0 s_i} \biggr), \label{bptolp} \tag{2.1}$$$

where $$\omega_0$$ is the central frequency, geometric mean of $$\omega_l$$ and $$\omega_h$$,

$$$\omega_0 = \sqrt{\omega_l \omega_u}$$$

$$\gamma$$ is the relative bandwidth,

$$$\gamma = \frac{\omega_u - \omega_l}{\omega_0}$$$

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$$:

$$$s_i = \omega_0 \biggl( \frac{\gamma}{2} S_i \pm j \sqrt{1 - \biggl( \frac{\gamma S_i}{2} \biggr)^2} \biggr)$$$

This implementation converts one of these poles to z-plane using bilinear transform,

$$$z_i = \frac{1 + s_i \frac{T}{2}}{1 - s_i \frac{T}{2}} \label{bilinear} \tag{2.2}$$$

and creates a second order section with poles at $$z$$ and $$z^*$$.

\begin{aligned} a_{i, 0} & = & 1, \\ a_{i, 1} & = & -2 Re(z_i), \\ a_{i, 2} & = & |{z_i}|^2. \end{aligned}

In the analog filter prototype, $$\omega_l$$ and $$\omega_u$$ are pre-warped frequencies of discrete filter:

$$$\omega = \frac{2}{T} tan \biggl( \pi f T \biggr)$$$

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