“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`:

`y[0] = 0.5 x[0] + 0.5 y[-1] = 0.5 * 1 + 0.5 * 0 = 0.5`

`y[1] = 0.5 x[1] + 0.5 y[0] = 0.5 * 2 + 0.5 * 0.5 = 1.25`

`y[2] = 0.5 x[2] + 0.5 y[1] = 0.5 * 3 + 0.5 * 1.25 = 2.125`

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.

via : Boris Smus