“You don’t understand anything until you learn it more than one way.” – Marvin Minsky

In my short Web Audio book, I covered the BiquadFilterNode, but didn’t
have any sense for how it worked. As I sat down to read Human and Machine
Hearing, it became clear that I needed to catch up on some digital
filtering fundamentals.

What follows is an introduction to digital filters via explorable
explanation I built to help myself better understand some DSP concepts.
The approach I took was to try to present the concept as visually and aurally as
possible, maximizing opportunities to build intuition. I learned a lot in the
process. Read on for a introduction, jump ahead to the Filter Playground,
or check out this video:

Cycles are everywhere

The world is full of examples of cyclical phenomena. Large cycles include
planetary motion, seasons, tides, and ocean waves. Societies are governed by
cycles: empires rise and fall, economies boom and bust, and fashion keeps
repeating itself. On an individual scale, the human lifecycle, eating and
sleeping, heartbeats and breathing, and locomotion are all periodic. And so are
sound and light, the very nature of the world we percieve.

Wouldn’t it be nice to understand and manipulate cyclical phenomena? A couple
hundred years ago, Fourier showed that any repeating periodic signal,
regardless of its complexity, can be represented as a sum of sine functions,
paving the way for much deeper signal understanding. More recently, Shannon
showed that you can take analog signals and represent them as numbers,
which brings us to digital signal processing. For modifying digital signals, we
need to look to digital filters.

The surprising link between periodic signals and difference equations

Bear with me on a brief mathematical tangent. Let’s start with a simple
difference equation that gives you an output sequence of numbers y[n] given
an input sequence x[n]:

Given `x = [1, 2, 3], a = b = 0.5`, we can calculate y[n] with simple
arithmetic, assuming that `y[-k] = 0`:

And so, given `x[n] = [1, 2, 3]` and the above difference equation, we get
`y[n] = [0.5, 1.25, 2.125]`. Computers can do this sort of arithmetic really
really quickly. But so what? And what does this have to do with our goal of
manipulating periodic signals in general? Let’s explore a bit more generally.

We can take any function x(t) and sample it numerically to get a sequence x[n].
For example, if we take `x(t) = t^2` and sample every integer from 1 to 10, we
would get `x[n] = [x(1), x(2), …, x(10)]`, in blue. Next, calculate y[n] and
color it green:

Let’s take several x(t), sample them into x[n], and see what the above
difference equation gives us for y[n]. When we plot x[n] in blue and y[n] in
green, we get the following graphs:

The first row shows the result of our difference equation on some mathematical
functions: `t^2`, `2^t`, `tan(0.1 t)`, and very little relationship
between x and y. The second row shows some sinusoidal functions with varying
period and phase in blue, sampled and run through the same difference equation
and shown in green. The pattern starts to become pretty clear: sine functions
are special in the same way. The generalized result is surprising and
awesome: if you take any sine function, and feed it into any difference
equation, you end up with another sine function with the same frequency, but a
different amplitude (gory details here).

We’ve looked at just one specific difference equation: `y[n] = 0.5 x[n] + 0.5
y[n-1]`. Let’s look at difference equations in general, and see how they
affect the relationship between x[n] and y[n].

Transfer functions describe the behavior of filters

We just saw that for an input sine function with a frequency and amplitude, a
difference equation will produce an output sine function with another amplitude.
In the `sin(0.1 t)` graph above, we can see that the amplitude of the green
graph is 0.5. If we look across all frequencies `x(t) = sin(omega t)` and see
how difference equations change the amplitude for all `omega`. This will give
us the frequency response for this specific difference equation. But what
happens in general?

If we bring our difference equation into a canonical form, we can apply the
Z-transform and get a corresponding transfer function. The math to
derive this is complicated, but the bottom line is that for any difference
equation in the general form:

The corresponding transfer function is of this form:

Let’s go back to our example, `y[n] = ax[n] + by[n − 1]`. We can
see that this fits the general form, with `a_0 = 1`, `b_0 = b`, `a_1 =
a`, with all of the other `a_i = b_i = 0`. Plugging in these values, the
transfer function for this example is:

This transfer function can tell us a lot about the behavior of the difference
equation, and ultimately its frequency response.

Visual intuition around transfer functions

But what the heck does this transfer function tell us?! Let’s try to build some
visual intuition. Firstly you see that our H(z) is a rational polynomial.
This means the numerator and denominator both have roots. Numerator roots are
called zeros and denominator zeros are called poles. From high school math,
recall that we can plot roots on a number line. Except in this case, we’re also
interested in imaginary roots, and we can plot them on the complex plane. O’s
represent zeros, X’s represent poles. For our H(z), the numerator is `0.5 z`,
which has a root (zero) at zero, and the denominator is `z – 0.5` with a root
(pole) at 0.5, and so plotting poles and zeros, we get:

Here’s how you can think about zeros and poles: if z is really close to a zero,
|H(z)| will be close to zero. If it’s really close to a pole, |H(z)| will blow
up, approaching infinity. Now, imagine you took a rubber sheet and draped it
over the complex plane. Now put stones where the O’s are, and telescoping tent
poles where the X’s are, and then extended the tent poles to be really tall. You
would end up with a circus tent. This is the surface that we get by plotting
complex values of z (along the x-y plane) and using |H(z)| as the height
corresponding to each z.

We can recover the pole-zero diagram by taking a birds-eye view at the same
circus tent, and plotting the poles and zeros on the xy-plane:

The neat thing about this circus tent is that it tells you the frequency
response of the filter. To do that, we look at the unit circle as it sits on top
of the circus tent diagram:

Which we can unwrap into a frequency response plot by looking at `|H(e^(i
omega))|` with `omega in [0, pi]`:

Now let’s see to what our example filter does to a white noise sample. Here is
the frequency response of the sample played through our filter according to a
Web Audio AnalyserNode (you can also hear the filter in the Filter
Playground):

Since noise has equal power across the frequency range, it is a good end-to-end
test. We expect the AnalyserNode‘s frequency response to line up closely with
the bode plot in the previous figure, and they do.

The filter we created above is a low pass filter, meaning that it allows low
frequencies to pass, but attenuates high frequencies. Let’s look at other kinds
of filters.

More complex filters

Our first filter had one pole and one zero. What happens if we make a filter
with two poles and two zeros? Here’s an example:

`H(z) = (z^2-1)/(z^2-1.975z+0.99) `

This example has two poles and two zeros, which is also called a Biquad filter.
The Web Audio BiquadFilterNode has 8 different filter types, all
of which are implemented with two poles and two zeros.

The filter in question is called bandpass filter, because it allows a band of
frequencies through, and attenuates the rest. Go ahead and open this in the
filter playground, and you’ll see a variety of views of this
bandpass filter.

The point of a playground isn’t just to look at other people playing, it’s to
play with them! So I invite you to try it out. You can generate a biquad filter
by selecting parameters with the filter wizard, or input any H(z) manually,
or move around poles and zeros visually. Using the pole-zero view, you can add
poles and zeros with buttons, or remove them by dragging them far enough out of
the unit circle. Check out my YouTube video for more examples of playing
around on the playground.

Implementation notes and thanks

I’d like to thank a handful of people for their help on this side project.
Firstly, to Raymond Toy for continuing to update the Web Audio API spec with
useful goodies. This side project wouldn’t be possible without the recently
added IIRFilterNode. Raymond has a few filter-related projects on the web,
including Digital Filter Design, which lets you create more complex
digital filters using a cascade of second order filters.

I built the 3D complex function plots with Mathbox, a very powerful
WebGL based visualization toolbox. If you haven’t seen it yet, check out How to
Fold a Julia Fractal, which is awesome in its own right, but also
illustrates the power of mathbox. It’s also a great introduction to complex
numbers, which I glossed over in this here post. Huge credit to Steven Wittens
for both mathbox and the inspiring blog post.

Finally, my thanks to Dick Lyon for writing an interesting and challenging
book and responding to my email (squee!), which ultimately inspired this
project.

A mathematical appendage

I tried to make the post understandable as possible by reducing analytical math
and leaning heavily on interactive illustrations. Inevitably I have waved my
hands and collected massive mathematical debt along the way, most notably
everything to do with deeply understanding the Z-transform, but
also:

• The filter playground requires poles to be inside the unit circle. If a
  pole is outside of the unit circle, the filter will become
  unstable, meaning that it will tend to blow up a signal that is
  fed into it.
• High order IIR filters tend to become numerically unstable. The
  solution is to break up high order rational functions into products of lower
  order rational functions, but the filter playground doesn’t currently do this.
• In most real life applications, complex filters are implemented as a cascade
  of first or second order filters chained together. This Digital Filter Design
  example illustrates it well.
• You might have noticed that points on the pole zero plot of the filter
  playground are never found floating alone on the complex plane. All of the
  coefficients of the numerator and denominator polynomials are real since they
  were taken from the original difference equation y[n], so by the complex
  conjugate root theorem, their roots (poles and zeros) must be either
  purely real, or appear in complex conjugate pairs.

Over and out

Thank you for getting this far, especially because I suspect this post may fall
into a sort of uncanny valley: too technical for a casual reader, and too
trivial for a DSP expert. At the very least, building the filter playground
helped me wrap my head around digital filters. Ultimately I hope the filter
playground can serve as a useful teaching tool for DSP novices.

I’d love to hear whether reading the post, watching the video and playing with
the filter playground helped you better understand digital filters. Please let
me know via twitter or by email.