Discrete-Time Transfer Functions and Stability

H(z) from difference equations, poles inside the unit circle, inverse transforms with residuez, FIR and IIR filters, and impulse responses from poles

1202 — Fibonacci's rabbits

In Liber Abaci, the Book of Calculation of 1202, Leonardo of Pisa — Fibonacci — poses a puzzle about an idealised population. A newly born pair of rabbits is put in a field; each pair mates at one month old and, from the end of its second month, produces a new pair every month; and rabbits never die. How many pairs will there be in a year (Fibonacci sequence)? Each month's count is last month's plus the month before's: 1, 1, 2, 3, 5, 8, 13, … The numbers had been described in Indian mathematics as early as 200 BC, by Pingala, in counting the patterns of Sanskrit poetry (Fibonacci sequence). In this unit's language the rabbits are a filter with feedback,

y[n] = y[n-1] + y[n-2] + x[n]

. Its transfer function has a pole at the golden ratio, 1.618 — outside the unit circle — which is why the population explodes.

A linear difference equation,

y[n] + a_1y[n-1] + \dots + a_Ny[n-N] = b_0x[n] + \dots + b_Mx[n-M]

, describes every digital filter. If no past outputs appear (

a_k = 0

), the impulse response ends after

M + 1

samples: a finite impulse response (FIR) filter, always stable. Feed outputs back and the response can go on for ever: an infinite impulse response (IIR) filter, as sharp as feedback allows — and as risky.

Derivation

From difference equation to H(z)

Transform both sides with the delay rule

y[n - k] \leftrightarrow z^{-k}Y(z)

, starting from rest:

Y(z)\left(1 + a_1z^{-1} + \dots + a_Nz^{-N}\right) = X(z)\left(b_0 + b_1z^{-1} + \dots + b_Mz^{-M}\right),

H(z) = \frac{Y(z)}{X(z)} = \frac{b_0 + b_1z^{-1} + \dots + b_Mz^{-M}}{1 + a_1z^{-1} + \dots + a_Nz^{-N}} .

Multiply numerator and denominator by

z^{\max(M,N)}

and the poles are the roots of

z^N + a_1z^{N-1} + \dots + a_N

. For the rabbits,

H(z) = 1/(1 - z^{-1} - z^{-2}) = z^2/(z^2 - z - 1)

, with poles at

(1 \pm \sqrt5)/2

1.618

and

-0.618

💡Stable means inside the unit circle

Each pole

p

contributes

p^n u[n]

to the impulse response, which decays only if

|p| < 1

. A causal filter is stable if and only if every pole lies strictly inside the unit circle. In this lesson's lab, after 40 samples:

$y[n] = 1.6y[n-1] - 0.81y[n-2] + x[n]$ , poles of radius 0.9: the response has fallen to 0.0066 — stable;
$y[n] = 1.6y[n-1] - y[n-2] + x[n]$ , poles of radius exactly 1: a sinusoid that never fades — marginally stable, a digital oscillator;
the rabbits, a pole at 1.618: 102 million pairs — unstable.

💡Inverse transforms by partial fractions

Expand

H(z)

in terms of

1/(1 - p_kz^{-1})

, each of which inverts to

p_k^n u[n]

H(z) = \sum_k \frac{r_k}{1 - p_kz^{-1}} \quad\Longrightarrow\quad h[n] = \sum_k r_k\,p_k^n\,u[n] .

For

H(z) = (1 + 0.5z^{-1})/\big((1 - 0.9z^{-1})(1 - 0.5z^{-1})\big)

, SciPy's residuez finds

r = 3.5

p = 0.9

and

r = -2.5

p = 0.5

, so

h[n] = 3.5(0.9)^n - 2.5(0.5)^n

: 1, 1.9, 2.21, 2.239, 2.140, … — exactly what lfilter computes sample by sample.

📐A 440 Hz resonator

A complex pole pair

re^{\pm j\Omega_0}

gives

h[n] \propto r^n\cos(\Omega_0 n + \phi)

: it rings at

\Omega_0

, and its envelope shrinks by a factor

r

each sample, falling to

1/e

after about

1/(1 - r)

samples. Put the poles at

0.99\,e^{\pm j2\pi\cdot 440/8000}

and, sampled at 8 kHz, the filter rings at 440 Hz — the A above middle C. Its envelope falls to

1/e

in about 96 samples (12 ms), and its peak is 25.6 Hz wide at −3 dB, close to the estimate

(1 - r)f_s/\pi = 25.5

Hz: two coefficients set both the pitch and the sharpness.

There from a marvelous instruction in the art of the nine Indian figures, the introduction and knowledge of the art pleased me so much above all else

– Leonardo of Pisa (Fibonacci), Liber Abaci (1202), on his schooling in Bugia, quoted on Wikiquote

🕰From Pingala to Pisa

The sequence was counted in India centuries before Fibonacci, in Pingala's work on Sanskrit prosody. Fibonacci himself learned arithmetic as a boy in Bugia, on the North African coast, where his father worked for Pisan merchants; he went on to study the methods of Egypt, Syria, Greece, Sicily and Provence, and judged the Indian numerals superior to them all (Fibonacci sequence; Fibonacci on Wikiquote).

⚠Rounding can move a pole outside

A filter whose poles sit at radius 0.999 is only a thousandth from instability. Store its coefficients with too few bits and the poles move — possibly across the circle. High-order IIR filters are therefore built as cascades of second-order sections, each with its own pair of poles (Lesson 9-3), rather than as one long difference equation.

💡Where you will use this

Frequency responses read from pole and zero positions are Lesson 8-3. FIR design is Lesson 9-1, IIR design Lesson 9-2, and the arithmetic of second-order sections Lesson 9-3. Digital oscillators with poles on the circle generate the carriers of Unit 10.

📚Further reading

Oppenheim and Willsky, Signals and Systems, sections 10.6–10.8 — the inverse z-transform and LTI systems described by difference equations.
Alan V. Oppenheim and Ronald W. Schafer, Discrete-Time Signal Processing (3rd ed., Pearson, 2010), chapters 3 and 5.
Fibonacci sequence on Wikipedia.

↓But there is a problem

Poles near the circle make peaks; zeros on it make silences. Every frequency response is the play of these points against the circle. And every ECG machine relies on it: the heart's signal is about a millivolt, and the building's wiring hums at 50 or 60 hertz.