DSP 106

Until now we’ve used Python to calculate the filter coefficients for the biquad sections of high-order filters. This is a quick way to implement a filter as it only requires the biquad code to actually use the filter. However, sometimes it’s necessary to have a filter with a variable cutoff and this method won’t allow that.

This acticle shows how to calculate the biquad filter coefficients from an analog prototype filter using the pole/zero mapping technique from DSP105.

The analog prototype filter will have a fixed cutoff of 1 radians per second. We will again use Python to obtain the (analog) poles and zeros for this filter and convert them into their digital equivalent using the so-called bilinear transform. This transform will also move the cutoff frequency of the prototype filter to the desired cutoff frequency. The biquad coefficients are easily derived from the digital poles and zeros.

Note that the analog prototype filter has a fixed cutoff and its poles and zeros can therefore be pre-calculated and easily stored in the signal processing code.

The prototype filter

In this example we will generate a Butterworth prototype filter. This filter has a maximally flat passband and it often the default choice if one needs a filter. This time, we will obtain the poles and zeros directly from signal.butter. Note that this time we supply analog=True to get the data for the analog filter, i.e. s-domain poles and zeros.

#!/usr/bin/python3

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

fs = 48000  # sample rate
fc = 1000   # desired cutoff of digital filter
N = 4       # filter order

# design an analog prototype filter
# zero locations are in 'z'
# pole locations are in 'p'
# k is the gain constant to get a 0 dB passband

z,p,k = signal.butter(N, 1, analog=True, output='zpk')

print("Zeros: ", f"{z}")
print("Poles: ", f"{p}")
print("Gain : ", f"{k}")

This will give the following output:

Zeros:  []
Poles:  [-0.38268343+0.92387953j -0.92387953+0.38268343j -0.92387953-0.38268343j
 -0.38268343-0.92387953j]
Gain :  1.0

Note that there are no zeros in the analog prototype filter and four poles. Like the digital poles from the previous article, the analog poles appear in complex conjugate pairs.

The bilinear transform

The bilinear transform maps the analog pole-zero space (the s-plane) onto the digital pole-zero space (the z-plane) or back. In addition it moves the analog cutoff frequency of 1 rad per second to our desired digital cutoff frequency.

For formula for the bilinear transform is: $$ a = 2 \pi \tan \left ( \pi \dfrac{f_c}{f_s} \right ) $$ $$ z = \dfrac{2+s\cdot a}{2 - s \cdot a} $$ where $s$ is the analog (s-domain) pole or zero to be transformed into a digital (z-domain) pole or zero $z$. Here, $f_c$ is the desired cutoff frequency in Hz and $f_s$ is the sample rate in samples per second.

In our example filter, we have the following pole pairs:

-0.38268343+0.92387953j and -0.38268343-0.92387953j

and

-0.92387953+0.38268343j and -0.92387953-0.38268343j

Let’s choose our cutoff to be $f_c = 1000$ Hz and $f_s = 48000$ samples per second. This means that: $$ 2 \pi \tan \left ( \pi \dfrac{1000}{48000} \right ) \approx 0.4118 $$ and the pole pairs transform to:

0.188968141+0.65155560j and 0.188968141-0.65155560j

and

0.136767685+0.19533078j and 0.136767685-0.19533078j

Even though it seems we have no analog zeros they do exist, namely at $s = 0.0\pm inf\cdot j $, If we look at the bilinear transform again, we see that we can write $z$ as $$ z = \dfrac{\frac{2}{s}+a}{\frac{2}{s} - a} $$ which gives $z=-1$ for these zero locations.

All that remains is to calculate the gain factor $k$. For a lowpass filter we want the gain at DC (0 Hz) to be 1. The gain of a biquad at DC is simply the sum of the b-coefficients divided by the a-coefficients. So, to compensate for that gain, we need $k = \dfrac{a_0 + a_1 + a_2}{b_0 + b_1 + b2}$ applied to each biquad.

Generating the biquad coefficients

Once transformed, the z-domain poles and zeros can be mapped using the formulas from DSP105: $$ b_0 = 1 $$ $$ b_1 = -2\cdot \real e(zero) $$ $$ b_2 = |zero|^2 $$ $$ a_0 = 1 $$ $$ a_1 = -2\cdot \real e(pole) $$ $$ a_2 = |pole|^2 $$

and then to set the gain right: $$ k = \dfrac{1 + a_1 + a_2}{1 + b_1 + b_2} $$ $$ b_0 = b_0 \cdot k $$ $$ b_1 = b_1 \cdot k $$ $$ b_2 = b_2 \cdot k $$

This procedure results in the following biquad coefficients:

[ 0.00407407,  0.00814814,  0.00407407,  1. , -1.88855595, 0.90485223]
[ 0.00381725,  0.00763449,  0.00381725,  1. , -1.76950435, 0.78477333]

Conclusion

The method above shows how to go from an analog filter prototype to a digital filter that can be implemented by biquads. In a practical application, a mathematical package such as Scipy is used to pre-calculate the poles and zeros of the analog prototype filter with a cutoff at 1 rad per second. These poles and zeros are then stored on the microcontroller or digital signal processor. A bilinear transform is applied to transform the poles and zeros into the z-domain using the bilinear transform. This also moves the filter cutoff to the desired frequency. The z-domain poles and zeros are then converted to biquad coefficients to implement the filter.

This method is flexible in that it can handle any kind of analog filter type, Butterworth, Bessel, Chebyshev and elliptical filters. The only thing that the implementer needs to be aware of is the gain normalization. In this article only the DC gain normalization was presented. Bandpass and highpass filters required a different normalization procedure.

Appendix

A python script that generates a 4th order Butterworth filter from an analog prototype.

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

def toBiquad(zero: complex, pole: complex, k: float) -> np.array:
    biquad = np.zeros(6)
    biquad[0] = k
    biquad[1] = -2*k*np.real(zero)
    biquad[2] = k*(np.real(zero)*np.real(zero)+np.imag(zero)*np.imag(zero))
    biquad[3] = 1
    biquad[4] = -2*np.real(pole)
    biquad[5] = np.real(pole)*np.real(pole)+np.imag(pole)*np.imag(pole)
    return biquad

# calculate the bilinear transform of location 's' 
def bilinear(s : complex, fc : float, fs : float) -> complex:
    a = 2.0*np.tan(np.pi*fc/fs)
    return (2+s*a)/(2-s*a)

fs = 48000  # sample rate
fc = 1000  # desired cutoff of digital filter
N = 4       # filter order

# design an analog prototype filter
# zero locations are in 'z'
# pole locations are in 'p'
# k is the gain constant to get a 0 dB passband

z,p,k = signal.butter(N, 1, analog=True, output='zpk')

# transform the poles
p_d = np.zeros(len(p), dtype=complex)
for idx, pole in enumerate(p):
    p_d[idx] = bilinear(pole, fc, fs)
    print(p_d[idx])

# map the poles to the biquads
# note: there are four zeros at -1 in the 
# z-domain
# we handle them directly in the
# toBiquad call

sos = []
bq1 = toBiquad(-1, p_d[0], 1)
bq2 = toBiquad(-1, p_d[1], 1)

# gain normalization for DC
k1 = (bq1[3]+bq1[4]+bq1[5]) / (bq1[0]+bq1[1]+bq1[2])
bq1[0] *= k1
bq1[1] *= k1
bq1[2] *= k1

k2 = (bq2[3]+bq2[4]+bq2[5]) / (bq2[0]+bq2[1]+bq2[2])
bq2[0] *= k2
bq2[1] *= k2
bq2[2] *= k2

sos.append(bq1)
sos.append(bq2)
print(sos)

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