DSP 102

This article shows how to implement second-order lowpass and highpass filters using the biquad filter structure.

In the previous article we discussed first-order digital filters. They are simple structures and a good starting point for doing experiments with digital signal processing. However, they have limitations compared to second-order filters. A second-order filter has a better roll-off response near the cutoff frequency and can have peaking or ‘resonance’.

The default structure for implementing second-order filters is the so-called biquad structure. More complex filters are generally built using multiple biquad sections, either in series or in parallel.

The biquad structure

This filter has one additional state compared to the first-order filter and two additonal coefficients. Note that the ‘a’ coefficients are now negative. It can be a lowpass or highpass filter, depending on the value of the coefficients. It really is a very nice generic structure.

The code for this structure is as follows:

tmp = input - a1*state1 - a2*state2
output = b0*input + b1*state1 + b2*state2
state2 = state1
state1 = tmp

For the mathematically inclined, the transfer function of this block is: $$ H(z) = \dfrac{b_0 + b_1\cdot z^{-1} + b_2\cdot z^{-2}}{1 + a_1\cdot z^{-1} + a_2\cdot z^{-2}} $$

Calculating the coefficients

The first-order filter has a nice approximation for the filter coefficient. The added complexity of the biquad doesn’t allow for such an approximation and software is used to calculate the coefficients. This is fine for static filters that don’t require real-time variability of the cutoff frequency. Multiple sets of coefficients can be pre-calculated and the desired set selected if some variability is required.

There exist many tools to calculate the biquad coefficients. I prefer using Python’s Scipy library.

import numpy as np
import scipy.signal as signal

fs = 48000 # the sample rate [Hz]
fc = 1000  # the desired cutoff frequency [Hz]

sos = scipy.butter(2, fc, output='sos', fs=fs)
print(sos)

The response of this script should be:

[[ 0.00391613  0.00783225  0.00391613  1.         -1.81534108  0.83100559]]

The first three numbers are the b-coefficients b0, b1 and b2. The last two coefficients are a1 and a2.

Python can also be used to plot the frequency response of this filter by using the Matplotlib library:

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

fs = 48000 # the sample rate [Hz]
fc = 1000  # the desired cutoff frequency [Hz]

sos = signal.butter(2, fc, output='sos', fs=fs)
print(sos)

# plot the frequency response
w,h = signal.sosfreqz(sos, fs=fs)
plt.semilogx(w, 20*np.log10(np.abs(h)))
plt.grid()
plt.ylim([-100,3])
plt.xlabel("Frequency [Hz]")
plt.ylabel("Magnitude [dB]")
plt.title("Frequency response")
plt.show()

More complex filters

More complex filters can be built using multiple biquads in series. Simply increase the filter order given to the signal.butter function in the script above. For example:

order = 6
sos = signal.butter(order, fc, output='sos', fs=fs)

The response is:

[[ 6.15535185e-08  1.23107037e-07  6.15535185e-08  1.00000000e+00 -1.76088036e+00  7.76074924e-01]
 [ 1.00000000e+00  2.00000000e+00  1.00000000e+00  1.00000000e+00 -1.81534108e+00  8.31005589e-01]
 [ 1.00000000e+00  2.00000000e+00  1.00000000e+00  1.00000000e+00 -1.91809148e+00  9.34642618e-01]]

Each line represents a biquad.

A highpass filter

The coefficients for a highpass filter are calculated by adding btype='high' to signal.butter`:

order = 6
sos = signal.butter(order, fc, output='sos', fs=fs, btype='high')

Similarly, other filters can be design using Python Scipy. Check out the documentation for more information.