import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
    matplotlib.RcParams._get = dict.get

9.6 Recursive filters

Every filter we have seen so far computes its output purely from the input. But there is no reason a difference equation cannot also refer to past outputs. Filters that do are called recursive filters, and they open up a large and powerful new class of behaviors.

Signal-flow diagrams

Recursive filters are often best understood visually, as a signal-flow diagram. This is yet another perspective on filters, complementing difference equations, convolution, and impulse responses. The diagrams are built from three elements: wires that carry a signal, a summing junction (drawn as a circled plus) that adds signals together, and a delay block labeled \(z^{-1}\), which delays its input by exactly one sample, mapping \(x[n]\) to \(x[n-1]\).

The notation \(z^{-1}\) comes from the z-transform, a generalization of the DFT that is the standard tool for analyzing recursive filters. We won’t cover the z-transform in this course, but we will still adopt the conventional \(z^{-N}\) notation in signal flow diagrams for a delay of \(N\) sample. Note that delaying by one sample is itself just convolution with the impulse response \(\color{red}{h} = [0, 1]\).

Fig. 55 Left: a feedforward filter, \(y[n] = x[n] + x[n-1]\), whose output depends only on the input (a delayed copy of the input is added in). Right: a feedback filter, \(y[n] = x[n] + y[n-1]\), whose output depends on itself (a delayed copy of the output is added in). The feedback loop is what makes a filter recursive.

The left diagram is an ordinary filter: the input flows forward through a delay and a sum to the output. The right diagram adds a feedback loop, “tapping” the output, delaying it, and feeding it back into the sum. This feedback is what makes the filter recursive. Contrasting the two side by side, the only difference is whether the delayed copy fed into the sum comes from the input (\(x[n-1]\), feedforward) or from the output (\(y[n-1]\), feedback).

Generalized difference equation

We can generalize the difference equation to include both past inputs and past outputs.

Definition 27 (General recursive difference equation)

A recursive filter is defined by

\[\begin{split} \begin{aligned} \purple{y[n]} = {}& b_0\,\blue{x[n]} + b_1\,\blue{x[n-1]} + \cdots + b_M\,\blue{x[n-M]} && \text{(feedforward)} \\ & \phantom{b_0\,\blue{x[n]}}{} + a_1\,\purple{y[n-1]} + \cdots + a_L\,\purple{y[n-L]} && \text{(feedback)} \end{aligned} \end{split}\]

The \(M{+}1\) feedforward coefficients \(b_i\) act on past inputs (like convolution before), while the \(L\) feedback coefficients \(a_j\) act on past outputs. Note there is no \(a_0\) term, since \(y[n]\) cannot depend on itself, only on past outputs.

The largest input and output delays are \(M\) and \(L\). The order of the filter is the larger of the two, \(\max(M, L)\): the number of samples of “memory” it must keep.

Two facts about this general form are worth committing to memory. First, every filter of this form is LTI, feedback and all. Second, its order \(\max(M, L)\) is the largest delay it uses. For example, \(y[n] = x[n] + x[n-1] + \tfrac{1}{3}y[n-2]\) has \(M = 1\) and \(L = 2\), so it is a second-order filter.

Recursive filters are commonplace in computer music because they can achieve higher-quality frequency responses with very few coefficients (and thus very little computation) compared to the equivalent non-recursive filter. More on “higher-quality” frequency responses in the next section!

Finite and infinite impulse responses

Feedback has a striking consequence for the impulse response. Consider the simplest recursive filter,

\[y[n] = x[n] + y[n-1].\]

What is its response to the unit impulse \(\delta = [1, 0, 0, \ldots]\)? We can read the output off the difference equation one sample at a time, recalling that \(y[n] = 0\) for \(n < 0\):

\[\begin{split} \begin{aligned} y[0] &= x[0] + y[-1] &= 1 + 0 = 1, \\ y[1] &= x[1] + y[0] &= 0 + 1 = 1, \\ y[2] &= x[2] + y[1] &= 0 + 1 = 1, \\ y[3] &= x[3] + y[2] &= 0 + 1 = 1, \\ &\;\;\vdots \end{aligned} \end{split}\]

The input contributes only its single initial \(1\), but the feedback keeps copying the previous output forward forever. The impulse response is \([1, 1, 1, 1, \ldots]\), infinitely long. This particular filter accumulates a running sum of its input.

This distinguishes two families of filters:

1. A filter with only feedforward coefficients has a finite impulse response (FIR) equal to the coefficients themselves. Its impulse response has as many nonzero samples as it has coefficients, and then stops.

2. A recursive filter (with feedback) generally has an infinite impulse response (IIR). The feedback keeps the response going forever.

With feedback comes a new danger: an IIR filter can be unstable. Compare two filters. The filter \(y[n] = x[n] + 0.9\,y[n-1]\) is stable: each pass through the loop shrinks the signal by a factor of \(0.9\), so its impulse response \([1, 0.9, 0.81, \ldots]\) decays toward zero. But \(y[n] = x[n] + 1.1\,y[n-1]\) is unstable: each pass amplifies the signal by \(1.1\), so its impulse response \([1, 1.1, 1.21, \ldots]\) grows without bound and quickly explodes into a deafening blowup. Designing stable recursive filters is a central concern of filter design.

For implementation, an IIR filter must generally be run as a difference equation, computing each output from previous outputs, rather than as a direct convolution (its impulse response is infinite, so we cannot convolve with all of it). That said, the impulse response of a stable IIR filter decays, so in practice we can approximate it by a finite one: run the impulse response until it has decayed below some threshold (say \(60\) dB down, \(|y[n]| \le 0.001\)), truncate it there, and convolve with the result. The example below computes the impulse response of the stable recursive filter \(y[n] = x[n] + 0.9\,y[n-2]\) and finds where it crosses the \(-60\) dB truncation threshold:

# hide
import numpy as np
import matplotlib.pyplot as plt
import pyquist as pq

# The stable recursive filter y[n] = x[n] + 0.9*y[n-2]. Its impulse response is
# a train of "echoes" 2 samples apart, each 0.9x the previous one, and it decays
# forever. To implement a stable IIR filter as a finite convolution, we truncate
# its impulse response once it falls below a threshold, typically -60 dB.
a, delay, N = 0.9, 2, 250
x = np.zeros(N)
x[0] = 1.0                                   # unit impulse
h = np.zeros(N)                              # the impulse response we build up
for n in range(N):
    h[n] = x[n] + a * (h[n - delay] if n >= delay else 0.0)

peak_n = np.arange(0, N, delay)              # the response peaks once per echo

plt.figure(figsize=(10, 4))
ml, sl, bl = plt.stem(np.arange(N), h)       # vertical bars + dots, one per sample
plt.setp(ml, color="C4", markersize=3)
plt.setp(sl, color="C4")
plt.setp(bl, visible=False)

# Mark a few decay thresholds. Each is db_to_amplitude(db) in linear amplitude.
# Because the decay is exponential, equal drops in dB happen at equal spacings.
for db, color in [(-20, "#f4a259"), (-40, "#e76f51"), (-60, "#c1121f")]:
    thresh = pq.helper.db_to_amplitude(db)
    cross = peak_n[np.argmax(np.abs(h[peak_n]) < thresh)]   # first echo below it
    plt.axhline(thresh, color=color, linestyle="--", linewidth=1.2)
    plt.axvline(cross, color=color, linestyle=":", linewidth=1.2,
                label=f"{db} dB (amp {thresh:g}) at n = {cross}")

plt.xlabel("Sample index n")
plt.ylabel("Impulse response amplitude")
plt.ylim(-0.06, 1.05)                        # offset below 0 so the lines clear the axis
plt.legend()
plt.show()

print("Equal drops in dB occur at equally spaced samples (exponential decay). "
      "Truncating at the -60 dB crossing keeps the audible tail.")

Equal drops in dB occur at equally spaced samples (exponential decay). Truncating at the -60 dB crossing keeps the audible tail.