IIR design by the bilinear transform

An analog filter becomes a digital one by a change of variable that bends the s-plane. Below, the jω axis curls into the unit circle, and a filter’s four poles ride along. Watch where each tick on the axis lands.

The jω axis curls into the unit circle

Analog prototype filters (20.1) gave us a shelf of good analog filters: Butterworth, Chebyshev, elliptic and Bessel. Each one is an $H(s)$, made for circuits. On this page I turn one into an $H(z)$, a digital filter that runs as code.

Analog and digital frequency meet on this page, so first a word on letters. Analog frequency is $\omega$ (rad/s), digital frequency is $\Omega$ (rad/sample); Oppenheim and Schafer use the opposite letters.

An analog filter is a ratio of two polynomials in $s$. So here is the plan: replace every $s$ by an expression in $z$, and the ratio becomes one in $z$. The replacement I’ll use is the bilinear transform:

Here $T_s=1/f_s$ is the sample period. The name says that the top and the bottom are both of first degree, linear, in $z$.

Does the order of the filter change? Take one pole $p_k$ of the analog filter. Its factor turns into

So each analog pole gives one digital pole and one zero at $z=-1$. A filter of order $N$ stays a filter of order $N$.

Where does a point $s$ land? Solve the transform for $z$:

Put in a point of the jω axis, $s=j\omega$. The top is $1+j\omega T_s/2$ and the bottom is $1-j\omega T_s/2$. They are mirror images, so they have the same length, and $\lvert z\rvert=1$. The whole jω axis lands on the unit circle.

The z-transform (16.1) had another map, $z=e^{sT_s}$. It wraps the jω axis round the circle again and again, once every $2\pi/T_s$ rad/s, and that is aliasing. The bilinear map goes round once.

Now the left half-plane. The top, $\lvert1+sT_s/2\rvert$, is the distance from the point $sT_s/2$ to $-1$, and the bottom is its distance to $+1$. A point left of the vertical axis is closer to $-1$, so $\lvert z\rvert<1$.

So every pole of a stable analog filter (Poles, zeros and the s-plane, 9.3) lands inside the unit circle. That is what a stable digital filter needs (Stability and causality, 16.4).

The picture at the top of the page shows the map as a motion. It starts from the s-plane scaled by $T_s/2$ and moved right by 1, which is the top of the fraction, $1+sT_s/2$. There the jω axis is the vertical line through 1. Then the division by $1-sT_s/2$ is faded in, and the line bends into the circle.

The filter on show is a 4th-order Butterworth at $f_s=8$ kHz, with its analog cutoff at 1054.8 Hz. That odd number is on purpose, and the second instrument explains it.

Watch the tick at 1054.8 Hz: it lands on 1000 Hz. And 16 kHz, far above $f_s/2=4$ kHz, lands at 3598 Hz, still below 4 kHz. Nothing aliases, because each analog frequency gets its own place on the circle.

Where each frequency lands

Follow a tick up the jω axis. For a very large $\omega$ the 1s hardly matter, and $z$ approaches $j\omega/(-j\omega)=-1$. So the far end of the axis, infinity, lands at $z=-1$, which is $f_s/2$.

A Butterworth $H(s)$ has only a constant on top, so its gain falls toward 0 as $\omega$ grows without end. We say it has its zeros at infinity, one for each pole: four here. The map sends them to $z=-1$, the $z+1$ in each factor above. So the digital gain is exactly 0 at $f_s/2$.

Where does a frequency in between land? The top, $1+j\omega T_s/2$, has the angle $\arctan(\omega T_s/2)$. The bottom has the same angle with a minus sign. Dividing subtracts angles, so the angle of $z$ is twice that:

This is warping. For a small number, $\arctan$ hardly changes it, so at low frequencies $\Omega\approx\omega T_s$, the $\Omega=2\pi f/f_s$ of Frequency in discrete time (12.1). To get hertz back, multiply $\Omega$ by $f_s/(2\pi)$.

At $f_s=8$ kHz, analog 500 Hz lands at 493.7 Hz, only 1.3 % low. High frequencies are squeezed instead. Everything from analog 4 kHz up to infinity fits between 2556.4 Hz and 4 kHz.

Pre-warp, and the cutoff lands where you asked

Warping moves frequencies, but it keeps the gains. Put a point of the unit circle, $z=e^{j\Omega}$, into the transform. Multiplying the top and bottom by $e^{-j\Omega/2}$ turns them into $2j\sin(\Omega/2)$ and $2\cos(\Omega/2)$:

So the digital gain at $\Omega$ is the analog gain at $\omega=\frac{2}{T_s}\tan\frac{\Omega}{2}$. Every value on the analog gain curve shows up on the digital one, at a moved frequency.

That is a problem for the cutoff. Design a Butterworth with its −3 dB point at 1000 Hz, map it, and the −3 dB point lands at 952.9 Hz. The fix is to aim off: ask where an analog frequency must be to land on the $\Omega$ we want.

The answer is the formula we just found, the landing formula turned round:

Designing the analog filter at $\omega_\text{pw}$ is pre-warping. For 1000 Hz at 8 kHz, $\Omega=2\pi\cdot1000/8000=\pi/4$, and $\omega_\text{pw}=16\,000\tan(\pi/8)=6627.4$ rad/s, which is 1054.8 Hz. That is the odd cutoff of the first instrument.

The instrument below designs two filters for each wanted cutoff. One has its analog cutoff set to the wanted one; the other has it pre-warped. It is like aiming at a stone under water: you aim where the bent light says, not where the stone seems to be. Watch where the two curves cross −3 dB.

Notice how the miss grows with the cutoff: 47 Hz at 1000 Hz, but 792 Hz at 3000 Hz. The squeeze is strongest near $f_s/2$, so a high cutoff moves furthest.

After the clip, drag the wanted line, or use the arrow keys, to move the wanted cutoff yourself. At 2000 Hz the plain design lands at 1695.4 Hz, and the pre-warped analog cutoff is 2546.5 Hz. Try 3500 Hz: the plain design lands at 2398.3 Hz, and the analog filter has to sit at 12 802.0 Hz. Then try 100 Hz, where the two designs nearly agree.

Pre-warping puts one frequency where you want it, and the rest of the curve still warps around it. For a low-pass that frequency is the cutoff. For a spec with two edges, you pre-warp both, as the worked example does.

The complete design

Now let’s design a whole filter by hand: a 4th-order Butterworth low-pass, −3 dB at 1 kHz, with $f_s=8$ kHz. There are four steps.

1. Pre-warp. $\omega_\text{pw}=2\cdot8000\tan(\pi/8)=6627.4$ rad/s.
2. Place the analog poles. 20.1’s Butterworth poles are $p_k=\omega_ce^{j\pi(2k+N+1)/(2N)}$ for $k=0$ to $N-1$. With $\omega_c=\omega_\text{pw}$ and $N=4$, their angles are 112.5°, 157.5°, 202.5° and 247.5°.
3. Map. Send each pole to $(1+p_kT_s/2)/(1-p_kT_s/2)$, and the four zeros at infinity to $z=-1$.
4. Set the gain. Choose the constant in front so that $H(1)=1$. The point $z=1$ is $\Omega=0$, so the gain at 0 Hz is 1.

The poles land at $0.7577\angle\pm42.73°$ and $0.4579\angle\pm20.94°$. All four zeros sit at $-1$, so the top of $H(z)$ is a constant $K$ times $(1+z^{-1})^4$, and $b_k$ is $K$ times 1, 4, 6, 4, 1. Step 4 gives $K=0.01021$.

Multiplying out the bottom gives the $a_k$. Here are both, as the coefficients of $H(z)$ in powers of $z^{-1}$, with $a_0=1$:

These are SciPy’s butter(4, 1000, fs=8000) to every digit shown. The bottom has $a_k$ beyond $a_0$, so the output feeds back on itself and the impulse response never quite ends: this is an IIR filter. butter with fs does all four steps for you.

Now the running spec of Filter specifications (18.1): the pass band to 1 kHz within ±0.05, the stop band from 1.5 kHz at 40 dB. A Butterworth peaks at gain 1, so the pass band allows a loss of $-20\log_{10}0.95=0.4455$ dB.

buttord(1000, 1500, 0.4455, 40, fs=8000) pre-warps both edges and returns order 12, with −3 dB at 1086.2 Hz. That design has gain 0.950 at 1 kHz, right on the limit, and −40.18 dB at 1.5 kHz. Its twelve poles sit at radii from 0.381 to 0.906.

Twelve poles is a lot to run in one piece. Second-order sections (21.2) shows how to run it safely.

Worked example

Let’s check the page’s numbers, by hand where we can and with SciPy where we can’t. Everything is at $f_s=8$ kHz, so $2/T_s=16\,000$ per second.

1. Where frequencies land. Use $\Omega=2\arctan(\omega T_s/2)$, then multiply by $f_s/(2\pi)$ for hertz. For 1000 Hz, $\omega T_s/2=\pi\cdot1000/8000=0.3927$, and $2\arctan(0.3927)=0.7487$ rad/sample, which is 952.88 Hz.

2. One pole by hand. The first analog pole is $6627.4\,e^{j112.5°}=-2536.2+6122.9j$ rad/s. Times $T_s/2$ it is $-0.15851+0.38268j$, of length $\tan(\pi/8)=0.41421$.

The top is $1+(-0.15851+0.38268j)=0.84149+0.38268j$, of length 0.92442 and angle 24.45°. The bottom is $1.15851-0.38268j$, of length 1.22008 and angle −18.28°.

Dividing, the length is $0.92442/1.22008=0.75767$ and the angle is $24.45°+18.28°=42.73°$. That is SciPy’s pole.

3. Back from z to s. Which analog frequency lands on $z=j$, a quarter of the way round, $f_s/4=2000$ Hz? Put $z=j$ into the transform: $(j-1)/(j+1)=j$, so $s=16\,000j$. That is 16 000 rad/s, or 2546.5 Hz, the pre-warped value for 2000 Hz in the second instrument.

4. The gain of the 4th-order design. At 0, 500, 1000, 1500 and 2000 Hz it is 0, −0.012, −3.010, −16.707 and −30.626 dB. The −3.010 dB at 1000 Hz is where we asked for it.

5. The order for the running spec, by hand. Pre-warp both edges: 1000 Hz gives 6627.4 rad/s and 1500 Hz gives 10 690.9 rad/s, a ratio of 1.6131.

20.1’s Butterworth gain satisfies $1/\lvert H\rvert^2=1+(\omega/\omega_c)^{2N}$. The pass edge needs $1/\lvert H\rvert^2\le1/0.95^2$, so $(\omega/\omega_c)^{2N}\le0.10803$ there. The stop edge needs gain at most 0.01, so $(\omega/\omega_c)^{2N}\ge9999$ there.

Divide the two conditions: the ratio of the edges, 1.6131, raised to $2N$, must reach $9999/0.10803=92\,555$. Taking logarithms, $2N\ge4.9664/0.20767$, so $N\ge11.96$ and the order is 12, as buttord said.

Where you’ll meet this

The bilinear transform is the standard way to turn an analog design into code. SciPy’s butter, cheby1, cheby2 and ellip use it for every digital design, and bilinear and bilinear_zpk do the map on its own. MATLAB’s design functions work the same way.

It is also how analog circuits become software: guitar tone stacks, synthesiser filters and old equalisers are modelled by writing down the circuit’s $H(s)$ and mapping it. Audio equalisers and biquads (20.5) builds the standard equaliser filters this way.

There is another way to sample an analog filter, Impulse invariance (20.3). High-pass, band-pass and band-stop shapes come from the same low-pass prototypes in Frequency transformations (20.4). Whether an IIR filter is the right choice at all is Choosing FIR or IIR (20.6).

Reference card