Second-order sections

A long IIR filter is safest run as a chain of small filters, each with two poles and two zeros. Below, an order-5 filter is split into its sections. Watch each section’s gain peak on its own, while the whole filter never rises above 0 dB.

Pair each pole pair with its nearest zeros

In Filter structures (21.1) we wired up one biquad in several ways, and saw that a long filter can be a cascade of them. This page asks how to split a long filter into biquads, and why it is worth doing. I’ll use the elliptic filter of Choosing FIR or IIR (20.6), which has order 5: five poles and five zeros.

Transfer functions, poles & zeros (16.3) wrote such a filter as a gain $K$ times one factor for each root:

Multiply the factors back together in small groups, and each group is a small filter. Run the small filters one after another, and by 21.1 the cascade is $H(z)$ again.

Which groups? A complex pole $p=re^{j\theta}$ comes with its mirror $p^*=re^{-j\theta}$, and their two factors multiply to

Both coefficients are real. In the same way, a pair of zeros $e^{\pm j\theta}$ on the unit circle gives $1-2\cos\theta\,z^{-1}+z^{-2}$.

One pole pair over one zero pair is the biquad of Audio equalisers and biquads (20.5). Used as one piece of a longer filter, I call it a second-order section. Each section has its own five coefficients:

An odd order leaves one real pole over. It goes with a real zero into a first-order section, $(b_0+b_1z^{-1})/(1+a_1z^{-1})$. So our order-5 filter becomes one first-order section and two biquads.

Which zeros go with which poles? That choice is the pairing. The arrows of 16.3 show why it matters. Where the walk round the circle passes close to a pole, that pole’s arrow is short and the gain climbs. A zero close by has a short arrow too, and it pulls the gain back down.

So give each pole pair the zeros nearest it, and no section climbs very high. It is like seating dinner guests next to the people they know, so that no table gets loud. SciPy’s zpk2sos(z, p, k, pairing='nearest') does it in this order:

1. Take the pole pair nearest the unit circle, and give it the two zeros nearest to it.
2. Repeat with the poles and zeros that are left.
3. A real pole takes a real zero, and the gain $K$ goes into the first section’s $b$‘s. (SciPy’s k is our $K$.)

The cascade gives the same $H(z)$ in any order, but SciPy still picks one. This is the ordering: the sections with poles nearest the circle go last. I’ll show why further down.

The picture at the top of the page splits 20.6’s filter this way, one section at a time, and draws each section’s gain. Notice that both biquads peak far above 0 dB on their own. Their peaks sit at different frequencies, and section 1 stays below −17 dB everywhere. Gains in a cascade multiply, so their dB add (20.5), and the sum never rises above 0 dB.

After the clip, choose the swapped pairing with its button: the loops redraw around the new groups, and the section curves change. The pole pair nearest the circle now gets the far zeros, and its section peaks at 22.99 dB; the other drops to 10.17 dB. The whole filter is unchanged.

The arrows of 16.3 explain the jump. At 1000 Hz the walk is at 45° on the circle, next to the poles of section 3. Their own zeros are 0.217 away, and the far zeros 0.541, two and a half times as far. Swapped, nothing near the poles pulls the gain down, and the peak rises by 9.0 dB.

Now the ordering. It does not change $H(z)$, but it changes the signals between the sections. In SciPy’s order the largest gain from the input to the end of section 1 is −17.43 dB, to the end of section 2 it is −4.11 dB, and to the output 0 dB. The sharp sections come last, so the signal reaches them already made smaller by the others.

Run the same three sections in the opposite order. After two of them, a 981 Hz sine comes out 24.47 dB louder than it went in, about 17 times. In floating point that does no harm. Finite word-length effects (21.3) shows what a stored value that large does in fixed point.

Raise the order until the direct form breaks

Why not skip the sections and run the whole filter as one difference equation, a direct form of 21.1? On paper the answer is the same. But a computer rounds every coefficient, and the direct form is far more sensitive to that rounding.

Most computers store a number in floating point: a sign, the first 24 binary digits of the number, and a power of 2 that says where the point goes. This is single precision, float32 in NumPy, and many DSP chips and audio plug-ins use it. 24 binary digits are about 7 decimal digits.

Rounding to single precision is the rounding to a grid of levels of Quantization & noise (11.1), with one difference: the step grows with the number. The error is at most $2^{-24}\approx6\times10^{-8}$ of the number’s size. NumPy’s usual float64, double precision, keeps 53 binary digits, about 16 decimal digits.

My test filter is a Butterworth low-pass from Analog prototype filters (20.1), with cutoff 100 Hz at $f_s=8$ kHz. 100 Hz is a small slice of the circle, 4.5°, so all its poles crowd near $z=1$.

I write the direct form’s denominator as $A(z)=1+a_1z^{-1}+\dots+a_Nz^{-N}$, with $N$ the order. For order 6 the coefficients are 1, −5.696561, 13.528499, −17.144063, 12.227073, −4.653138 and 0.738190. Put $z=1$, and $A(1)$ is their sum. By the factor form, it is also the product of the six small numbers $1-p_k$.

So the coefficients, up to 17 in size, add up to $A(1)=2.0\times10^{-7}$. Rounding −17.14406324 to single precision gives −17.14406395, a change of $7.1\times10^{-7}$. That one change is larger than $A(1)$ itself, so every digit of the coefficients matters.

The instrument below rounds the coefficients of this filter to single precision in two ways, as one direct form and as second-order sections. Then it finds the poles of what was rounded. The filtering itself runs in double precision, so the rounded coefficients are the only cause of what you see.

Notice that the rings stay on the true poles at every order. The direct form’s crosses drift at order 4 and leave the circle at order 6.

After the clip, set the order yourself with the slider. Try 12: the direct form’s largest pole is 1.5037, while the sections’ is 0.9898.

Order 6 breaks in a way you can see in the numbers. After rounding, the coefficients add up to $A(1)=-1.19\times10^{-6}$: the sign has flipped. Multiply $A(z)$ by $z^6$, and you have a polynomial in $z$. Far out on the real axis it is positive, and at $z=1$ it is now negative.

So it crosses zero somewhere beyond $z=1$. That crossing is a real pole outside the circle, the 1.0417 of the readout. By Stability and causality (16.4), one pole outside the circle is enough to make the filter unstable.

Why the direct form is so sensitive

Write $P(z)=z^NA(z)$, the same polynomial in powers of $z$, so that $P(z)=\prod_k(z-p_k)$. Rounding the coefficients adds a tiny polynomial to $P$. Near one root $p_i$, $P$ is close to a straight line, and its slope there is the product of the distances from $p_i$ to the other $N-1$ roots.

The root moves by about the change in $P$ at $p_i$ divided by that slope. When the roots crowd together, that slope is tiny. For order 6, the pole nearest the circle is 0.039, 0.075, 0.106, 0.131 and 0.149 away from the other five, and the product is only $6.0\times10^{-6}$.

The rounding changes $P$ there by $1.3\times10^{-6}$, so the estimate is a move of about 0.2. That pole is 0.020 inside the circle, ten times closer than that. The estimate only holds for small moves, but it says clearly that the pole will not stay put.

A biquad’s two roots are moved only by its own two coefficients, and the slope at each root is the single distance to its mirror. In the order-6 design that distance is 0.149 for the section nearest the circle, so its poles move by $2.1\times10^{-7}$.

Double precision has 29 more binary digits, so it lasts longer, but it breaks too. The same Butterworth design as one direct form keeps every pole inside the circle up to order 11. At order 12 its largest pole is at 1.0241, with no rounding to single precision at all.

So I keep an IIR filter above order 2 or 3 as sections, from the design to the filtering. SciPy’s butter documentation says the same: the [b, a] form can have numerical problems “even for N >= 4”, and it recommends second-order sections, output='sos'.

Worked example

Let’s redo the page’s numbers, by hand where we can and with SciPy where we can’t.

1. The elliptic filter in sections. ellip(5, 0.4455, 40, 1000, fs=8000, output='sos') returns one row per section, with $b_0,b_1,b_2,1,a_1,a_2$. Each row is one section of the form above, and the filter is their product, $H(z)=H_1(z)H_2(z)H_3(z)$.

Check section 3 against its roots. The poles have radius 0.952677 and angles ±45.636°, so $2r\cos\theta=1.332241$ and $r^2=0.907593$. The zeros are at ±57.439°, so $2\cos\theta=1.076389$.

2. The dB add up. At 0 Hz, $z=1$, so each section’s gain is the sum of its $b$‘s over the sum of its $a$‘s. Section 1 gives $0.044892/0.334084=0.134374$, which is −17.43 dB. Section 2 gives $1.529560/0.329942=4.6358$, 13.32 dB, and section 3 gives $0.923611/0.575352=1.6053$, 4.11 dB.

The sum is $-17.43+13.32+4.11=0.00$ dB, the gain of the whole filter at 0 Hz.

3. The swapped pairing at 0 Hz. Section 2 now has section 3’s zeros: $0.923611/0.329942=2.7993$, 8.94 dB. Section 3 has section 2’s zeros: $1.529560/0.575352=2.6585$, 8.49 dB. Together they give 17.43 dB, the same as $13.32+4.11$.

So at 0 Hz the swap is hard to see. The cost is near 1 kHz, where the swapped section 3 peaks at 22.99 dB instead of 14.00.

4. Coefficient sizes. The order-6 Butterworth denominator at 100 Hz has coefficients up to 17.14 in size, and the order-12 one up to 681. Rounding to single precision changes each by at most $2^{-24}\approx6\times10^{-8}$ of its size. Here the largest such change is $4.8\times10^{-8}$ at order 6 and $5.3\times10^{-8}$ at order 12.

5. Where double precision breaks. With no rounding to single precision, the direct form’s largest pole is 0.9876 at order 10 (true 0.9878) and 0.9944 at order 11 (true 0.9889). At order 12 it is 1.0241, and the impulse response first passes 1 at $n=504$.

6. In SciPy. Design straight into sections with butter(6, 100, fs=8000, output='sos'), and filter with sosfilt(sos, x), which runs each section as the transposed direct form of 21.1. If you have poles and zeros, zpk2sos(z, p, k) pairs and orders them.

tf2sos(b, a) also exists, but it starts from the polynomial, so it cannot undo rounding already done to $b$ and $a$. Applied to the single-precision coefficients of order 6, it gives sections with the same pole at 1.0417.

Where you’ll meet this

Running IIR filters as cascades of biquads is the usual practice. SciPy has sosfilt and sosfiltfilt, and MATLAB has zp2sos and sosfilt. Arm’s CMSIS-DSP library for microcontrollers has the arm_biquad_cascade functions. The Web Audio API’s BiquadFilterNode is one section, and you chain several for a higher order.

The parametric equalisers of 20.5 are already cascades, one biquad for each band. How many bits each coefficient and stored value needs is Finite word-length effects (21.3). Running the sections sample by sample, against a clock, is Real-time processing (21.4).

Reference card