Difference equations

The echo loop of What is a system? (4.1), fed a single 1. Watch each stem come from the one before.

Each answer comes from the last one

Remember the feedback loop in What is a system? (4.1). Set its gain to one half. A sound goes out and comes back half as loud, and that echo goes round again and comes back half as loud again. Each echo is made from the one before it.

I can write that rule as an equation. Let $n$ be the sample number, $x[n]$ the input at sample $n$, and $y[n]$ the output. The symbol ${y[n-1]}$ is the output one step earlier, the shift from Shifting, reversing and scaling time (2.1). It is the number held in the delay block from 4.1, which I call the delay box; the picture at the top of the page labels it “one step back”. The echo loop says

In words: take the number sitting in the delay box, halve it, and add the new input. An equation that gives each output sample from the input and from earlier output samples is a difference equation. The input here is a single 1 at $n=0$ and 0 for ever after, the unit sample $\delta[n]$ from Impulse, step and ramp (3.1), and the box starts at 0.

Make a guess: what is $y[3]$? Then go back to the loop at the top of the page, with Step back and Step forward if you like, and watch the stem for $n=3$ appear: it is $0.125$, half of the $0.25$ the box held.

The input arrived once, and the output goes on. That is the loop remembering. Each stem is half the one before, so the stems are $1,\ 0.5,\ 0.25,\ 0.125,\dots$, which is $y[n]=0.5^n$.

The answer to a difference equation is not one number. It is the whole row of stems, one output for every $n$. Running the rule forward, one sample at a time, is how I solve it.

This row is also something you already know. The input was the unit sample and the box started empty, so the row is the impulse response $h[n]$ from The impulse response (5.1). The next sections explain why a box that starts empty matters.

With a loop and without one

Smoothing is a good job for a difference equation. Suppose a thermometer gives a noisy reading each second, and I want a steadier number. There are two ways to do it.

The first way is to average the last three readings. This is the 3-point moving average:

The second way is to keep a running estimate and nudge it halfway toward each new reading. The nudge gives $y[n]={y[n-1]}+0.5\,(x[n]-{y[n-1]})$, which simplifies to

I call this the leaky integrator. An integrator adds up its input as it goes, as the running sum in Operations on amplitude (2.2) does, and keeps all of its old total. This one keeps only half of the old total at each step, so the total leaks.

The two equations differ in a way you can see in a diagram. The moving average uses only inputs, so its diagram has no wire from the output back to the input side. We call it non-recursive. The leaky integrator uses its own past output, so its diagram has a loop. We call it recursive.

Feed both smoothers the same single 1 at $n=0$, and look at $n=3$.

The moving average has three stems of $\tfrac13$, then a stem of 0 at $n=3$, and stays at 0. The leaky integrator is at $0.0625$, still not 0. By $n=7$ its stem is $0.00390625$, tiny but not zero.

This is the reason behind the FIR and IIR of Properties of LTI systems (5.4). Without a loop, the output has no way to remember the input beyond the three samples it uses, so once the input has passed, the output stops. With this loop, the output feeds itself, so it never quite stops.

What is stored, and what is pushed in

Think of a bath. The water in it at any moment comes from two places: the water that was already there, and the water the tap has added since. The output of a difference equation has the same two sources.

Take $y[n]=0.8\,{y[n-1]}+x[n]$. Each step keeps $0.8$ of the old output and adds the new input. Before the first sample, the box already holds some number. I write it $y[-1]$, the output one step before the start, and I call it the stored value.

Here is a question to answer before the next instrument. There is a 2 in the box, so $y[-1]=2$, and then a 1 arrives. What is $y[0]$?

I want to pull that answer apart into the part from the 2 and the part from the 1. The way to do it is two experiments. First, switch the input off and let the stored 2 run on its own. Second, empty the box and switch the input on.

The first experiment gives the zero-input response, the output caused by the stored value alone. The second gives the zero-state response, the output caused by the input alone when the box starts empty. A system that starts with every stored value equal to 0 is said to be at initial rest, which for this equation means $y[-1]=0$.

Watch the two experiments run, then stack.

With the input off, the stored 2 shrinks: $1.6,\ 1.28,\ 1.024,\dots$. With the box emptied, the input is the unit step $u[n]$ from Impulse, step and ramp (3.1), switched on and left on, and the stems climb: $1,\ 1.8,\ 2.44,\dots$. Last, the two sets of stems stack.

Look at $n=0$ in the stacked picture. The open-circle stem is $1.6$ tall. The stem stacked on top of it is the zero-state part, $1$ long. The diamond at its top is their sum, $2.6$. That is $0.8\cdot2+1$, the answer to the question above. The same holds at every $n$: the total is the two stems added.

Why can the two be added? Each output sample is made by multiplying by fixed numbers and adding. This is superposition, from System properties (4.2): the response to a sum of causes is the sum of the responses to each. Here the two causes are the stored value and the input.

There is one more reason to care about initial rest. Scale the input by 0 and a linear system must give 0, which is the homogeneity half of the linearity test. With a 2 in the box, the output is not 0, because the stored 2 keeps going. So a system that starts with something stored fails the linearity test, and we usually start from rest. From rest the equation describes a linear, time-invariant system, and everything in The impulse response (5.1) applies.

Both parts have exact formulas. The stored part is multiplied by $0.8$ at every step, so after $n+1$ steps from $y[-1]$,

For the input part, run the rule step by step from rest with $x=u[n]$. Then $y[0]=1$, $y[1]=0.8+1$, $y[2]=0.8^2+0.8+1$, and so on, so $y_\text{zs}[n]=1+0.8+\dots+0.8^n$. Here I use the finite geometric sum, which I take from algebra rather than derive:

With $a=0.8$ this gives $y_\text{zs}[n]=5\,(1-{0.8}^{n+1})$. The total is $y=y_\text{zi}+y_\text{zs}$, and the worked example checks the first four values against direct iteration.

Every term is a wire

So far I have written equations and drawn diagrams as if they were two things. They are one thing, and this section shows how. Take the equation

It adds the input, twice the previous input, and half the previous output. There are three terms on the right.

Each term becomes one wire into the adder. A term with ${x[n-1]}$ needs a delay box on its wire, one box for each step back, and a gain triangle for its number. A gain of 1 is not drawn. A term with ${y[n-1]}$ sends a wire back from the output.

Watch each term light up and draw its wire into the diagram.

The wire back from the output is the loop. A ${y[n-1]}$ wire starts at the output and runs back through a delay box. A ${y[n-k]}$ term always sends such a wire, which is exactly what recursive means. If no term has a $y$, there is no loop and the equation is non-recursive.

Read the other way, a diagram gives its equation: write one term for each wire into the adder. A wire that starts at the input is an $x$ term, and one that starts at the output is a $y$ term. The number of delay boxes on a wire is how many steps back it reaches, and its gain triangle is the coefficient. A wire with no gain triangle has a gain of 1. The worked example below does this for a diagram.

The general form

Now I can write every such equation at once. On this page $a_k$ are the numbers that multiply past outputs and $b_k$ are the numbers that multiply inputs. (The page Fourier series coefficients (7.2) uses the same letters for a different job.) The general difference equation is

with the first coefficient fixed at

Here $N_a$ and $N_b$ are the largest steps back in the outputs and the inputs. The output terms are on the left, so $y[n]=x[n]+0.5\,{y[n-1]}$ has $a_1=-0.5$: moving $0.5\,{y[n-1]}$ to the left changes its sign. This is the convention of the SciPy function lfilter(b, a, x).

Solve for $y[n]$ and you get the rule to run forward:

It needs the stored values $y[-1],\ y[-2],\dots$, which are zero at initial rest. The equation is recursive when some $a_k$ with $k\ge1$ is not zero.

Worked example

1. Iteration. Take $y[n]=0.5\,{y[n-1]}+x[n]$ with $x=\delta[n]$ and initial rest. Running forward gives 1, 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625, 0.0078125 for $n=0$ to $7$. This is $y[n]=0.5^n$, so $y[5]=0.03125$.

2. Moving average and leaky integrator. Let $x=3,\ 6,\ 9$ at $n=0,\ 1,\ 2$, and 0 after. The 3-point moving average gives $y=1,\ 3,\ 6,\ 5,\ 3,\ 0$ for $n=0$ to $5$. For example $y[2]=(9+6+3)/3=6$. Its last nonzero output is at $n=4$, two samples after the last input, and it is 0 from $n=5$ on: no loop. The leaky integrator $y[n]=0.5\,{y[n-1]}+0.5\,x[n]$ with $x=\delta[n]$ gives $0.5,\ 0.25,\ 0.125,\ 0.0625,\dots$ and $0.00390625$ at $n=7$, never exactly 0.

3. Zero-input plus zero-state. Take $y[n]=0.8\,{y[n-1]}+u[n]$ with $y[-1]=2$. The zero-input part is $y_\text{zi}[n]=2\cdot{0.8}^{n+1}$, which gives 1.6, 1.28, 1.024, 0.8192 for $n=0$ to $3$. The zero-state part is $y_\text{zs}[n]=5\,(1-{0.8}^{n+1})$, which gives 1, 1.8, 2.44, 2.952. The totals are 2.6, 3.08, 3.464, 3.7712. Direct iteration agrees:

   $$\begin{aligned} 0.8(2)+1&=2.6\\ 0.8(2.6)+1&=3.08\\ 0.8(3.08)+1&=3.464\\ 0.8(3.464)+1&=3.7712 \end{aligned}$$

4. Equation to impulse response. Take $y[n]=x[n]+2\,{x[n-1]}+0.5\,{y[n-1]}$ with $x=\delta[n]$. Then $h[0]=1$. For $n=1$ the three terms are $x[1]=0$, $2\,x[0]=2$ and $0.5\,h[0]=0.5$, so $h[1]=2.5$. After that only the loop acts, so each value is half the one before: $h$ is 1, 2.5, 1.25, 0.625, 0.3125. In the card’s convention, $b_0=1$, $b_1=2$, $a_0=1$ and $a_1=-0.5$. In SciPy this is lfilter([1, 2], [1, -0.5], x).

5. Diagram to equation. Fig. 1 shows a diagram. It has three wires into the adder, so the equation has three terms. The straight wire from the input gives $x[n]$. The wire from the input through one delay box and a gain of $-1$ gives $-{x[n-1]}$. The wire from the output through one delay box and a gain of $0.9$ gives $0.9\,{y[n-1]}$. Together,

   $$\begin{aligned} y[n]&=x[n]-{x[n-1]}\\ &\quad+0.9\,{y[n-1]}, \end{aligned}$$

   and it is recursive, because it has a loop. In the card’s convention $b_0=1$, $b_1=-1$ and $a_1=-0.9$. From rest, a unit sample gives $h[0]=1$ and $h[1]=0-1+0.9\cdot1=-0.1$. After that each value is $0.9$ times the one before: $h$ is 1, −0.1, −0.09, −0.081.

Where you’ll meet this

Most digital filters in audio effects, sensors and radios run a difference equation, one sample at a time, in a loop like the ones on this page. The leaky integrator is one of the cheapest smoothers: one stored number and two multiplications per sample. The code that runs these filters carries the stored values from one block of samples to the next and calls them the filter’s state. The stored values produce the zero-input part of this page, and SciPy’s lfiltic turns them into the zi argument that lfilter takes.

The next page, Differential equations and analog systems (6.2), does the same job for a circuit, where time is not counted in samples. There the stored part and the pushed-in part appear again, along with a second way to cut the answer in two. Then First- and second-order systems (6.3) names the number that sets how fast the loop here shrinks.

Reference card