Impulse invariance

One way to copy an analog filter is to sample its impulse response. Below, the impulse response of an analog filter is drawn and then sampled. Watch where the two digital poles land on the plane at the end.

Sample the analog impulse response

In IIR design by the bilinear transform (20.2) we borrowed an analog filter and bent its frequency axis onto the unit circle. Here is a second way to borrow one. It starts from time instead of frequency: take the analog filter’s impulse response and sample it.

As in 20.2, analog frequency is $\omega$ (rad/s), digital frequency is $\Omega$ (rad/sample); Oppenheim and Schafer use the opposite letters.

I write the analog filter’s impulse response $h_c(t)$ and its transfer function $H_c(s)$. The subscript c means continuous time, as it does for a sampled signal $x_c(t)$. The impulse-invariant filter takes one sample every $T_s$ seconds and scales it by $T_s$:

Why the factor $T_s$? The gain at 0 Hz is the area under $h_c(t)$ for the analog filter and the sum of $h[n]$ for the digital one. Samples spaced $T_s$ apart, each times $T_s$, add up to nearly that area, so the two gains stay close.

Now split $H_c(s)$ into one piece per pole, as in Properties and the inverse Laplace transform (9.2). This works when $H_c$ has more poles than zeros, as every Butterworth filter does. Each piece is one exponential:

Sampling at $t=nT_s$ turns $e^{p_knT_s}$ into $\big(e^{p_kT_s}\big)^n$, so every piece becomes a geometric sequence:

A geometric impulse response is one pole (Properties of LTI systems, 5.4). Summing its geometric series (Difference equations, 6.1) gave the pair $a^nu[n]\leftrightarrow1/(1-az^{-1})$ of The z-transform (16.1). One piece at a time, that gives

Each analog pole $p_k$ becomes the digital pole $e^{p_kT_s}$. That is 16.1’s map from the s-plane to the z-plane, $z=e^{sT_s}$.

Let’s try it on one filter: the order-2 Butterworth of Analog prototype filters (20.1), with $f_c=1$ kHz, sampled at $f_s=8$ kHz. So $\omega_c=2\pi\cdot1000=6283.2$ rad/s and $T_s=0.125$ ms. Its two poles sit on the half-circle of radius $\omega_c$, at $\pm135°$:

The transfer function is $H_c(s)=\omega_c^2/\big((s-p_1)(s-p_2)\big)$, which has gain 1 at $s=0$ because $p_1p_2=\omega_c^2$. The cover-up rule gives $r_1=\omega_c^2/(p_1-p_2)=\omega_c^2/(j\sqrt2\,\omega_c)$. That is $-j\omega_c/\sqrt2=-4442.9j$, and the lower pole has the mirror weight $r_2=+4442.9j$.

9.2’s mirror pair adds the two exponentials into a damped sine:

The picture at the top of the page samples this curve. Think of a flip-book of a falling ball. Every page shows the ball where it really was at that moment; the pages only skip the moments in between.

Notice the angle. The digital poles sit at 31.8°, and that is $\omega_cT_s/\sqrt2$ in degrees: the analog pole’s spin rate times one sample period. Their distance from 0, 0.574, is $e^{-\omega_cT_s/\sqrt2}$: the decay over one sample period.

The poles follow $z=e^{sT_s}$. The zeros do not follow a simple rule, because they come out of adding the pieces. The analog filter has no finite zeros. Adding the two digital pieces over a common denominator gives

with 6.1’s coefficient names: $b_1=0.336071$, $a_1=-0.975239$ and $a_2=0.329322$.

The weights cancel in the constant term, $r_1+r_2=0$, so $h[0]=0$, just as $h_c(0)=0$. In powers of $z$ it reads $b_1z/(z^2+a_1z+a_2)$, so it has a zero at $z=0$, which no analog zero maps to.

You have met this design once already. In Differential equations and analog systems (6.2), the RC circuit’s discrete copy was $y[n]=a\,{y[n-1]}+(1-a)\,x[n]$ with $a=e^{-T_s/RC}$. The circuit’s impulse response is $\tfrac1{RC}e^{-t/RC}$ (Frequency response and Bode plots, 8.4): one pole at $-1/RC$ with weight $1/RC$.

Impulse invariance turns it into

The pole is 6.2’s $a$, by the same rule. Only the factor on top differs: 6.2 has $1-a$ where impulse invariance has $T_s/RC$. 6.2 matched the circuit’s step response at the end of each interval, a variant called step invariance.

For 6.2’s $RC=1$ s and $T_s=0.5$ s the two factors are 0.393 and 0.5. They agree when $T_s$ is small next to $RC$.

The analog gain, folded

Impulse invariance copies the impulse response at the sampling instants. What does that do to the gain? The impulse response is a signal like any other. By The sampling theorem (10.2), its samples have a spectrum made of copies of the analog one, one at every multiple of $f_s$.

10.2’s copies carried a factor $1/T_s$, and the $T_s$ in $h[n]=T_s\,h_c(nT_s)$ cancels it. With $\omega_s=2\pi f_s$,

So the digital response at a frequency is the analog response there, plus the analog response at every frequency a multiple of $f_s$ away. The copy centred at $f_s$ reaches down into 0 to $f_s/2$. It carries the analog gain from above $f_s/2$, mirrored about $f_s/2$.

I call this the folding of the frequency response. It is 10.2’s aliasing, applied to a frequency response.

Think of the analog gain written on a long strip of paper, folded back and forth into a short box. What was written past the fold lands on top of what is already there. The copies add as complex numbers, with their phases, not as plain gains.

The bilinear design of 20.2 has no copies. It squeezes the whole analog frequency axis onto 0 to $f_s/2$, so the analog gain at infinity lands at $f_s/2$. For a Butterworth that gain is 0: the order-2 bilinear design has a double zero at $z=-1$. Below, the analog filter and its two digital copies are drawn together; watch the gain at 4 kHz as the cutoff rises.

After the clip, drag the dotted $f_c$ line, or use the arrow keys, to move the cutoff yourself. Try 1500 Hz. The impulse-invariant gain at 4 kHz is −10.4 dB, while the analog gain there is −17.1 dB. Now try 3500 Hz: −4.9 dB against the analog −4.3 dB. Here the copies partly cancel, because they add with their phases.

Impulse invariance or the bilinear transform

Impulse invariance copies a time response at the sampling instants, and it keeps frequencies where they were, with no warping. That makes it useful when the impulse or step response is the specification. Control engineers use it that way, and so do people who model an analog circuit’s transient.

It works only for filters whose analog gain is already small above $f_s/2$. That means low-pass and band-pass filters, at a high enough $f_s$. A high-pass or band-stop filter is never small at $f_s/2$, so its copies overlap at full strength. For a specification on the gain, the bilinear transform is the default: its frequencies are warped, but nothing folds.

The gain at 0 Hz drifts too, by −0.45 dB for our 1 kHz filter. To fix it, divide $H(z)$ by its gain at 0 Hz, which is its value at $z=1$.

Worked example

Let’s redo the page’s numbers.

1. The order-2 Butterworth at 1 kHz, sampled at 8 kHz. The analog poles are $-4442.9\pm4442.9j$ rad/s, with weights $\mp4442.9j$, which is $\mp j\omega_c/\sqrt2$. Over one sample period the spin turns $4442.9\cdot0.000125=0.5554$ rad, which is 31.82°, and the decay shrinks the size by $e^{-0.5554}=0.5739$. So the digital poles are 0.5739∠±31.82°, or $0.4876\pm0.3026j$.

The samples are a damped sine, $h[n]=1.1107\,(0.5739)^n\sin(0.5554\,n)$, because $T_s\sqrt2\,\omega_c=1.1107$. 16.1’s damped-sine pair, with radius 0.5739 and angle 31.82°, gives the coefficients of $H(z)$:

Its gain is −0.453 dB at 0 Hz, −3.015 dB at 1 kHz and −16.72 dB at 4 kHz, where the analog filter has −24.10 dB.

2. The folding grows with the cutoff. I designed the impulse-invariant filter for eight cutoffs and read its gain at 0 Hz and at 4 kHz. The analog gain at 4 kHz is there for comparison.

Up to 1 kHz the folded gain at 4 kHz stays more than 16 dB down. At 2 kHz the gain at 0 Hz is already 1.9 dB off. From 2.5 kHz on, the gain at 4 kHz is at most 1.5 dB below the gain at 0 Hz, and from 3 kHz it is above it.

Where you’ll meet this

SciPy’s cont2discrete(system, dt, method='impulse') and MATLAB’s impinvar do this design. Control engineers use it to give a digital model the impulse or step response of an analog plant. Audio programmers use it to copy an analog circuit’s ringing, and to build sounds from sums of decaying sines, one pole pair per sine.

For shapes other than low-pass, see Frequency transformations (20.4). Whether an IIR filter is the right choice at all is Choosing FIR or IIR (20.6).

Reference card