The unilateral z-transform

One recursion, run from rest and then with numbers stored. Watch the stored values add a second term on top of the input’s part, which does not move.

Starting a recursion at n = 0

On Difference equations (6.1) you ran a recursion with a number already sitting in its delay box. You split the result into the zero-input response, made by the stored value alone, and the zero-state response, made by the input alone. There you did the algebra by hand. This page puts both through one transform.

The bilateral transform of The z-transform (16.1) adds up a signal over every $n$. A recursion switched on at $n=0$ does not fit that. The equation holds from $n=0$ onward, and before that I only know the stored numbers, $y[-1]$ and $y[-2]$. They are the past, not part of the equation.

So I use the unilateral transform, which sums from $n=0$ only:

It is the same letter as before, and this page says “unilateral” when it matters. For a signal that is zero before $n=0$ the two transforms are the same function.

What changes is the delay rule. Take $x[n-1]$ and split off its first term, at $n=0$, which is $x[-1]$. The rest is $z^{-1}$ times the transform of $x$. Doing the same for $x[n-2]$ gives two terms before the pattern starts:

The extra terms are the samples that slide in from before $n=0$. This is the same job that $y(0^-)$ did in Properties and the inverse Laplace transform (9.2): the stored values enter the algebra through the delay rule.

Now apply it to $y[n]=y[n-1]-0.5\,y[n-2]+x[n]$ with a step input, $x[n]=u[n]$, switched on at $n=0$. Transform every term and collect the $Y(z)$ terms on the left:

Write $D(z)=1-z^{-1}+0.5z^{-2}$. Dividing by $D(z)$ shows that $Y$ has two parts. The input’s part is $X(z)/D(z)$, which is the zero-state response $y_\text{zs}$. The stored values’ part is $(y[-1]-0.5\,y[-2]-0.5\,y[-1]z^{-1})/D(z)$, which is the zero-input response $y_\text{zi}$. The total is $y=y_\text{zi}+y_\text{zs}$, as on page 6.1.

Both parts share the denominator $D(z)$, so both are made by the same poles. They are the roots of $z^2-z+0.5$, which are $0.5\pm0.5j$. Each has radius 0.7071 and angle $\pi/4$ (45°), which is 8 samples per cycle.

SciPy turns stored values into the state that lfilter starts from: zi = lfiltic(b, a, [y[-1], y[-2]]), then lfilter(b, a, x, zi=zi). I use it below to check every number.

Stored values add their own term

The picture at the top of this page runs this recursion. From rest, with $y[-1]=y[-2]=0$, only the input acts: the output climbs to 2.5, rings, and settles at 2. The zero-state samples are 1, 2, 2.5, 2.5, 2.25, 2, 1.875, 1.875, 1.9375, 2 for $n=0$ to $9$.

Store $y[-1]=2$, and a zero-input part starts at 2 and rings away to 0, with the same poles, on top of the input’s part. The total starts at 3 and 3. The input’s part did not change; a second term came in.

Store 2 in both, and the zero-input part exactly cancels the climb: the total is 2 from the start. The stored values were already the final value.

When the clip ends, drag the stored squares to set the stored values yourself. Try doubling both stored values: the zero-input stems double, and the zero-state stems stay where they were.

What the algebra predicts

Each frame follows from the algebra above. With $y[-1]=2$ and $y[-2]=0$, the stored values’ part of $Y$ is

Its samples are 2, 1, 0, −0.5, −0.5, −0.25, 0, 0.125 for $n=0$ to $7$. The inverse transform of 16.2 gives a closed form, with the weight 1 at each of the two poles:

At $n=3$ it gives $2\cdot0.3536\cdot(-0.7071)=-0.5$, as the stems show.

The zero-state part is $X(z)/D(z)=1/\big((1-z^{-1})D(z)\big)$. Partial fractions (16.2) give three terms. One is at the step’s pole, 1. The other two are a pair at $0.5\pm0.5j$ that dies away:

The weight 2 is the final value, and the dashed line in the instrument. Written in real terms, $y_\text{zs}[n]=2+2^{-n/2}\big(\sin(n\pi/4)-\cos(n\pi/4)\big)$. The total for $(2,0)$ is 3, 3, 2.5, 2, 1.75, 1.75 for $n=0$ to $5$.

The last frame can be solved exactly. With $y[-1]=y[-2]=2$, the stored values’ numerator is $2-1-z^{-1}=1-z^{-1}$. The input’s term is $X(z)=1/(1-z^{-1})$. Adding them and using $2D(z)=2-2z^{-1}+z^{-2}$ gives

That is $y[n]=2$ for every $n\ge0$. You can check it without a transform: if every value is 2, the equation reads $2=2-1+1$. The two stored values were already the output the equation settles at.

Switched on: a part that dies, a part that stays

The second split is about where the poles come from. For a system $H(z)$ and an input with transform $X(z)$, the output is $Y(z)=H(z)X(z)$ from rest. Its partial fractions have terms at the poles of $H$ and terms at the poles of $X$.

Take a sine switched on at $n=0$, $x[n]=\sin(\Omega_1n)\,u[n]$. Its transform has poles on the unit circle, at $e^{j\Omega_1}$ and $e^{-j\Omega_1}$. The terms at those poles make the steady state, $y_\text{ss}[n]$. It is the sine that Frequency response of discrete-time systems (12.4) taught you to expect: the same frequency, scaled by $\lvert H\rvert$ and shifted by $\angle H$, both read at $e^{j\Omega_1}$:

The terms at the poles of $H$ make the transient, $y_\text{tr}[n]$. For a stable system those poles are inside the unit circle (Stability and causality, 16.4), so this part dies away. The output is the sum, $y=y_\text{tr}+y_\text{ss}$.

The system here has one pole: $y[n]=0.8\,y[n-1]+0.2\,x[n]$, so $H(z)=0.2/(1-0.8z^{-1})$. Its gain at $\Omega=0$ is 1. The transient is $y_\text{tr}[n]=c\cdot0.8^n$. The output starts from rest, so $y[0]=0$, and the constant must cancel the steady state at $n=0$: $c=-y_\text{ss}[0]$.

Watch the transient’s open circles shrink by 0.8 each sample, until only the steady sine is left.

At $n=0$ the input is 0, and so is the output. The steady sine would be $-0.222$ there, so the transient starts at $+0.222$ to cancel it.

By $n=6$ the transient has shrunk to 0.058, by 0.8 each sample, and the total is mostly the steady sine.

By $n=15$ the transient is 0.008. What is left is $0.280\sin(0.25\pi n-52.5°)$, 12.4’s gain and phase for this test. Here $\lvert H\rvert=0.280$ and $\angle H=-0.916$ rad, which is −52.5°.

By $n=30$ the transient is gone. The system’s pole made the part that dies; the input’s poles made the part that stays. The transient first falls below 0.01 at $n=14$: it is 0.0122 at $n=13$ and 0.0098 at $n=14$.

When the clip ends, drag across the output to choose $\Omega_1$ yourself, from $0.05\pi$ to $0.95\pi$. At $0.10\pi$ the steady sine has size 0.581 and phase −45.9°, and the transient starts at 0.418.

Try $0.95\pi$ and then $0.05\pi$. The transient still shrinks by 0.8 each sample, since that is the pole’s doing, but its starting size changes. At $0.95\pi$ it starts at 0.008, and at $0.05\pi$ at 0.419. The start is set by the phase: a steady sine with phase near zero already starts near 0, so little needs cancelling.

I chose a sine so that the transient is big enough to see. A switched-on cosine at $0.25\pi$ gives $c=0.029$ instead, because the output starts at 0.2 and the steady cosine starts near 0.171.

What the split tells you

A stable system’s transient dies away, because its poles are inside the circle. What remains is the steady state, and the steady state is all that a measurement of $\lvert H\rvert$ and $\angle H$ ever sees. Frequency response of discrete-time systems (12.4) said “after the sine has played for a long time”. This is what that phrase meant.

The same split gives a rule for a step input. A step has one pole, at 1, so the steady part is a constant. Its size is $H(1)$, and for a stable $H$

For the system of the first instrument, $H(1)=1/(1-1+0.5)=2$, the dashed line. For the system of the second, $H(1)=0.2/(1-0.8)=1$.

Worked example

Take the first system, $y[n]=y[n-1]-0.5\,y[n-2]+u[n]$, with $y[-1]=1$ and $y[-2]=0$.

1. Stored values’ numerator. $y[-1]-0.5\,y[-2]-0.5\,y[-1]z^{-1}=1-0.5z^{-1}$. In SciPy, lfiltic([1], [1, -1, 0.5], [1, 0]) returns zi = [1, -0.5].
2. Zero-input part. This is half of the $(2,0)$ case, by linearity: 1, 0.5, 0, −0.25, −0.25, −0.125, 0, 0.0625 for $n=0$ to $7$.
3. Zero-state part. The same as before: 1, 2, 2.5, 2.5, 2.25, 2.
4. Total. Add them: 2, 2.5, 2.5, 2.25, 2, 1.875 for $n=0$ to $5$.

Iterate the equation directly as a check. $y[0]=1-0+1=2$. $y[1]=2-0.5\cdot1+1=2.5$. $y[2]=2.5-0.5\cdot2+1=2.5$. $y[3]=2.5-0.5\cdot2.5+1=2.25$. These match the total.

For the second system, find the steady state at $\Omega_1=0.25\pi$ by hand. $H(e^{j\pi/4})=0.2/(1-0.8e^{-j\pi/4})$. Since $0.8\cos(\pi/4)=0.8\sin(\pi/4)=0.5657$, the denominator is $0.4343+0.5657j$, with magnitude 0.7132 and angle $0.9160$ rad. So $\lvert H\rvert=0.2/0.7132=0.2804$ and $\angle H=-0.9160$ rad $=-52.5°$.

The transient size is $c=-\lvert H\rvert\sin(\angle H)=0.2804\cdot\sin(52.5°)=0.2224$. At $n=6$ the transient is $0.2224\cdot0.8^6=0.0583$ and the steady state is $0.2804\sin(6\pi/4-0.9160)=-0.1708$, so $y[6]=-0.1125$. The output $y[1]$ to $y[4]$ from lfilter is 0.141, 0.313, 0.392, 0.314.

Where you’ll meet this

Every streaming filter carries its stored values from one block to the next. The code calls them the filter’s state, and lfilter_zi and zi are how SciPy sets them up. Real-time processing (21.4) picks this up. The transient is also what you hear and see when a notch or a resonator is switched on, which Resonators, notches and combs (17.4) returns to.

Reference card