DSP 105

Analog and digital filters are defined by their poles and zeros. These are locations in a two-dimensional space where the z-domain response goes to infinity (in the case of a pole) or zero (in the case of a zero).

This article shows how to calculate the biquad coefficients for a given set of poles and zeros. This is the first step in designing custom filters on-the-fly.

Obtaining poles and zeros

The pole or zero locations are specified by a two-dimensional number, for example (1,3), where the 1 is the value in one dimension and 3 is the value in the second dimension. To show a location in a plot, the first value is plotted on the x-axis, and the second on the y-axis.

In electrical engineering poles and zeros are represented by a complex number, itself being a two-dimensional quantity consisting of two parts: the real part and the imaginary part. Rather than (1,3) it is written as 1 + 3i or 1 + 3j, where the multiplication by i or j denotes the imaginary part.

We can use Python to calculate the poles and zeros for a fourth-order Butterworth filter and print them:

import numpy as np
import scypi.signal as signal

N  = 4          # filter order
fs = 48000      # sample rate
fc = 4000       # filter cutoff frequency
z,p,k = signal.butter(N, fc, fs=fs, output="zpk")

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

The script produces the following output:

Poles:  [0.72693283+0.3877475j  0.59238104+0.13088208j 0.59238104-0.13088208j
 0.72693283-0.3877475j ]
Zeros:  [-1. -1. -1. -1.]
Gain:   0.0025764344253226203

Note that for every pole (or zero) with a positive imaginary part there is also a pole with the same negative imaginary part. This is not accidental. In fact it is a necessary property to be able to get biquad coefficients! A pole (or zero) without an imarginary part can appear on its own, however.

In figure 2 the pole and zero locations have been plotted.

Biquad revisited

The biquad implements the following transfer function: $$ H(z) = k \cdot \dfrac{(z-zero_1)(z-zero_2)}{(z-pole_1)(z-pole_2)} $$

This looks a lot different than the biquad transfer function seen everywhere: $$ H(z) = \dfrac{b_0 + b_1\cdot z^{-1} + b_2\cdot z^{-2}}{a_0 + a_1\cdot z^{-1} + a_2\cdot z^{-2}}$$

However, one can be written as the other, which gives a way to calculate the coefficients $a_0,a_1,a_2$ and $b_0,b_1,b_2$.

We start by expanding the numerator and denominator of the first equation: $$ H(z) = k \cdot \dfrac{z^2 -zero_1 \cdot z -zero_2 \cdot z + zero_1 \cdot zero_2}{z^2 -pole_1 \cdot z -pole_2 \cdot z + pole_1 \cdot pole_2} $$

The divide numerator and denominator by $z^2$: $$ H(z) = k \cdot \dfrac{1 -zero_1 \cdot z^{-1} -zero_2 \cdot z^{-1} + zero_1 \cdot zero_2 \cdot z^{-2}}{1 -pole_1 \cdot z^{-1} -pole_2 \cdot z^{-1} + pole_1 \cdot pole_2 \cdot z^{-2}} $$

It’s already starting to look familiar! Now the pole (and zero) pairs come into play; $pole_1$ and $pole_2$ must have real parts that are equal and imaginary parts that are equal but opposite in sign.

In the example above, we have two pole such pole pairs:

0.72693283+0.3877475j and 0.72693283-0.3877475j

and

0.59238104+0.13088208j and 0.59238104-0.13088208j

Let’s call the pole with the positive imaginary part $pole$ and the one with the negative imaginary part $pole^*$. The star indicates a complex conjugate – a fancy way of saying the sign of the imaginary part has been changed.

To continue we need the following two equivalents: $$ pole+pole^* = 2\cdot\real e(pole) $$ $$ pole \cdot pole^* = |pole|^2 $$ Here $\real e(x)$ means taking the real value of $x$. And $|x|^2$ meants the magnitude squared of $x$, which is the real part squared and the imaginary part squared added together. It is important to note that both equivalents have a real part only!

For example: $$ (1+2j)+(1-2j) = 2\cdot\real e(1+2j) = 2 \cdot 1 = 2 $$ and $$ (1+2j) \cdot (1-2j) = |(1+2j)|^2 = 1^2 + 2^2 = 5 $$

Because we have a pole pair $pole_1 = pole, pole_2 = pole^* $, we can write: $$ H(z) = k \cdot \dfrac{1 -zero \cdot z^{-1} -zero^* \cdot z^{-1} + zero \cdot zero^* \cdot z^{-2}}{1 -pole \cdot z^{-1} -pole^* \cdot z^{-1} + pole \cdot pole^* \cdot z^{-2}} $$

And, using the equivalents above, we arrive at:

$$ H(z) = k \cdot \dfrac{1 -2 \cdot \real e(zero) \cdot z^{-1} + |zero|^2 \cdot z^{-2}}{1 -2 \cdot \real e(pole) \cdot z^{-1} + |pole|^2 \cdot z^{-2}} $$

By comparing the above formula to $$ H(z) = \dfrac{b_0 + b_1\cdot z^{-1} + b_2\cdot z^{-2}}{a_0 + a_1\cdot z^{-1} + a_2\cdot z^{-2}}$$ we can almost directly read off the coefficients: $$ b_0 = k $$ $$ b_1 = -2\cdot k\cdot \real e(zero) $$ $$ b_2 = k\cdot |zero|^2 $$ $$ a_0 = 1 $$ $$ a_1 = -2\cdot \real e(pole) $$ $$ a_2 = |pole|^2 $$

Python implementation

The following Python script provides a function toBiquad that implements the coefficient calculations derived above. It uses the complex Python data type to represent a complex number. Most modern programming languages have support for complex numbers. C++, for example, has std::complex<>. Rust has support through the num crate.

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

# returns a the biquad coefficients in SOS format
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

N  = 4          # filter order
fs = 48000      # sample rate
fc = 4000       # filter cutoff frequency

# get the poles and zeros for a Butterworth filter
z,p,k = signal.butter(N, fc, fs=fs, output="zpk")

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

# convert poles and zeroz to biquad coefficients
# p[0] and p[3] are a pole pair and so are
# p[1] and p[2]. We only need p[0] and p[1].
#
# in the same fashion, we only need z[0] and z[1]
#
# note that we only need to take care of the gain
# 'k' once!
#
sos = []
bq1 = toBiquad(z[0], p[0], k)
bq2 = toBiquad(z[1], p[1], 1)   # no need to take care of the gain here.

print("SOS =",sos)

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

The script above gives in the following results:

SOS =
[[ 0.00257643, 0.00515287, 0.00257643, 1. , -1.45386566, 0.67877946]), 
 [ 1.        , 2.        , 1.        , 1. , -1.18476209, 0.36804542])]