NOTE: A .pdf version of this with much more readable equations is attached.

This tutorial is going to be a crash course in going from a difference-equation representation of a filter to a filter using Pd's "raw" filters. The raw filters are [rzero~], [rpole~], [czero~], and [cpole~] (there's also [rzero-rev~] and company, but we won't be using them here). It won't go into detail on everything and will likely be overwhelming if you're new to filter design or the math involved, but hopefully it will get you started, help you understand some basic concepts needed for filter design, and/or at least clear up some things.

NOTATION

I'll first start with a simple first-order difference equation to explain the notation I'll be using and to give the basics of what you should be seeing in the equation:

y[n] = b_0*x[n] + b_1*x[n-1] + a_1*y[n-1]

x[n] (read: "x of n") refers to input, and y[n] refers to output. n is the current sample, n-1 is the previous, and so forth. b_m (read: "b sub m", I'm using the underscore for this since subscripts don't work on this forum) are the coefficients of x[n-m], and a_m are the coefficients of y[n-m]. Using b_m for x[n-m] and a_m for y[n-m] is the standard convention that you will typically see in academic papers.

So another way to read this equation is "the current output is b_0 times the current input plus b_1 times the previous input plus a_1 times the previous output". As such, a_m are commonly referred to as the feedback coefficients (since they multiply the output feeding back into the filter), and b_m are referred to as the feedforward coefficients.

One way to look at x[n-1] is to think of it as the input delayed by one sample. Likewise, y[n-2] would be the output delayed by two samples. So, as you can see, the m in b_m and a_m matches up with the number of delayed samples (i.e. b_m*x[n-m]), and m is the delay in samples. The order of the filter is the maximum delay m. So if you had a difference equation with x[n], x[n-1], and y[n-2], it would be a second-order filter, because the maximum delay is 2 samples.

Also, we're going to be discussing complex numbers here (numbers with a real and imaginary part). While you may have used i = sqrt(-1) in math classes, I will be using j instead. j is used in engineering instead of i, and since DSP is more of an engineering thing, you will come across j more often when researching this stuff.

In addition, the ^ will be used for exponents. x^2 is "x to the power of 2", or [expr pow(x, 2)].

RAW FILTERS

Okay, now that we have notation out of the way, let's quickly look at the difference equations for the raw filters.

[rzero~] and [czero~]: y[n] = x[n] - b_1*x[n-1]

[rpole~] and [cpole~]: y[n] = x[n] + a_1*y[n-1]

NOTE THAT THE SIGNS ARE DIFFERENT. If you don't remember this, it will bite you in the ass later. They are different in the difference equation, but they line up in the z-transform, which we will discuss in the next section. (Don't know why they're negative, though.)

The difference between [rzero~] and [czero~] is that [rzero~] uses only real numbers for b_1, and the input is also real. [czero~], on the other hand, uses complex numbers, meaning it has a real and imaginary part. The real and imaginary parts are represented as separate inputs. So for a complex number x + j*y, x would go in the left inlet and y would go in the right (you don't have to calculate j...I mean, it isn't even real). The same is true with [rpole~] and [cpole~]. While most filters for audio are real filters, we will see later why the complex filters are important.

Z-TRANSFORM

The z-transform converts a difference equation in the time-domain to a transfer function in the z-domain. The usefulness of the z-domain is beyond the scope of this tutorial. But, you don't need to know everything about it for it to be useful here. Just accept that it works, just like you have to accept that imaginary numbers work and that infinity exists at least conceptually even though our brains are too small to comprehend, and you'll be fine. You'll be happy to know that doing the z-transform is much easier than wrapping your head around imaginary numbers, anyway.

The basic representation of a generic filter in the z-domain is this:

B(z)
H(z) = --------
A(z)

Again, B(z) is read "B of z". The difference between square brackets [] and parenthesis () is that brackets represent signals with discrete points (like a digital sampling rate) and parenthesis represent continuous signals (like analogue). It's not a huge deal here, it's just another convention.

As you might have guessed, B(z) is named because it represents the part of the difference equation with b_m coefficients, while A(z) represents the part with the a_m coefficients. H(z) is just another naming convention for filters. h[n] is the impulse response in the time-domain, so H(z) is the z-transform of the impulse response.

So, how do we get from the time-domain to the z-domain? Well, it's actually quite simple, but it's probably easiest to explain with an example. So let's use the biquad filter as an example. The difference equation of a conventional digital biquad is this:

y[n] = b_0*x[n] + b_1*x[n-1] +b_2*x[n-2] - a_1*y[n-1] - a_2*y[n-2]

(Sidenote: [biquad~] actually reverses the signs for the feedback coefficients)

Each of the samples in the equation gets replaced by an exponent of z where the exponent is the delay length. Let's start by just figuring out B(z). In this case, we get:

B(z) = b_0 + b_1*z^-1 + b_2*z^-2

Pretty straight-forward. x[n-m] becomes z^-m. (z^0 = 1, which is why b_0 stands alone here.)

Now for the feedback part: A(z). With this part, y[n-m] becomes -(z^-m). Notice the signs get reversed. Also, a_0 is actually the number that multiplies y[n], which is 1, and it's sign doesn't get reversed because it's on the left side of the difference equation. This gives us:

A(z) = 1 + a_1*z^-1 + a_2*z^-2

So our final z-transform of the biquad filter is:

b_0 + b_1*z^-1 + b_2*z^-2
H(z) = ----------------------------------
1 + a_1*z^-1 + a_2*z^-2

So what we have is a second-order polynomial B(z) over a second-order polynomial A(z).

CONVERT TO FIRST-ORDER FILTERS IN SERIES

The biquad is a second-order filter. The raw filters in Pd are first-order filters. We need to represent this second-order filter using nothing but first-order filters. So, we need to break it down into first-order polynomials. To do this, we need to find the poles and zeros.

Well, let's back up a second. How can we even break this apart at all? As it happens, multiplication in the z-domain is the same thing as filters running in series in the time-domain. As a very simple example, look at it this way:

B(z) B(z) 1
H(z) = ------ = ------ * ------
A(z) 1 A(z)

In this instance, B(z) is one filter, and 1/A(z) is another. The signal goes through B(z) first, then 1/A(z). But it has exactly the same output as B(z)/A(z).

Since we're working with polynomials here, we need to factor them out into first-order polynomials to get first-order filters in series. To do this, we need to find the roots. The roots of B(z) are the zeros, and the roots of A(z) are the poles. It is generally less confusing to find them if we rewrite them as functions of z instead of z^-1. We can do this by multiplying H(z) by z^2/z^2 = 1:

b_0*z^2 + b_1*z + b_2
H(z) = -----------------------------
z^2 + a_1*z + a_2

Next, we need to factor out b_0 from B(z). It's not necessary for finding the roots, but it will be needed later for implementation purposes. This is simply a matter of dividing it by b_0, which we will rename g for gain. We'll rename b_m/b_0 to β_m.

z^2+ β_1*z +β_2
H(z) = g * ------------------------
z^2 + a_1*z + a_2

What we now have is a quadratic equation in both the numerator and the denominator. Finding the roots of those equations (i.e., the value of z that will make them equal zero) can simply be done using that quadratic formula we all thought was useless in high school:

For y = a*x^2 + b*x +c, the roots are

-b ± sqrt(b^2 - 4*a*c)
y = -----------------------------
2*a

This translate for the zeros to:

-β_1 ± sqrt(β_1^2 - 4*β_2)
q = ------------------------------------
2

and the poles to:

-a_1 ± sqrt(a_1^2 - 4*a_2)
p = ------------------------------------
2

Note that there are two answers each because of the ±. Also note that if b^2 - 4*a*c < 0, you are going to end up with complex numbers as your zero or pole, and this is going to happen often. Once we find the zeros/poles, H(z) becomes:

(z - q_1)(z - q_2) (1 - q_1*z^-1)(1 - q_2*z^-1)
H(z) = g * --------------------- = g * ------------------------------------
(z - p_1)(z - p_2) (1 - p_1*z^-1)(1 - p_2*z^-1)

The last step just undoes the z^2/z^2 multiplication we did earlier to bring us back to negative powers of z. But the answers are the same. If you substitute q_1 for z, you will get 0 in both versions.

AND NOW FOR THE PD IMPLEMENTATION

As stated earlier, multiplication in the z-domain is the same thing as running filters in series. We now have four first-order polynomials and a gain. This is what we need to use the raw filters. Breaking it down into multiplications, we get:

g * (1 - q_1*z^-1) * (1 - q_2*z^-1) * (1/(1 - p_1*z^-1)) * (1/(1 - p_2*z^-1))

Assuming the poles and zeros are complex, this is the same thing as:

[inlet~]
|
| g
| /
[*~ ]
| q_1
| //
[czero~]
| /
| / q_2
| / //
[czero~]
| /
| / p_1
| / //
[cpole~]
| /
| / p_2
| / //
[cpole~]
|
|
|
[outlet~]

With a real filters, of course, you would just use the real versions and get rid of the imaginary patch cords. But with typical higher order audio filters, the zeros and poles will come in complex conjugate pairs, so you can expect to use the complex filters. Besides, a real number x = x + j*0, anyway, so the complex ones will work in both situations.

With filter orders greater than 2, the quadratic formula obviously doesn't hold. However, it is quite common for high-order filters to be implemented as a series of biquads because they are easier to work with that way (Butterworth and Chebychev filters work nicely as biquad series). So the formula may still be useful. In addition, programs like Octave or Matlab have functions that can find the roots of arbitrary polynomials for you for those difficult cases.

PHEW!

That might be a lot to swallow in such a short space. For something more thorough, I recommend the online book The Scientist and Engineer's Guide to Digital Signal Processing by Steven W. Smith, which actually uses very accessible language despite the title. It's also more than just filters. I should point out, however, that the author does go against convention and switches b_m and a_m, so watch out. But other than that, it's pretty great. Julius O. Smith III also has a good book on digital filters online. It's a little less accessible, but there's a lot there.

Hopefully that didn't suck?

http://www.pdpatchrepo.info/hurleur/filter_tut_1.pdf

Thanks a lot for the detailed write-up, filter-guru Maelstorm. Much better than the regular concise description 'find the points where H(z) is zero, you're done'.

Still, I am confused. Can't a complex first order feedback filter be considered a real second order feedback filter? If you follow the signal flow in a cpole~ from real input to output, and leave out the imaginary input, the one sample delay for imaginary input functions as a two sample delay for the real input.

The illustration is from a page on complex resonators I did some years ago:

http://www.katjaas.nl/complexintegrator/complexresonator.html

If my question makes no sense, I'll remove it as to not spoil your tutorial.

Katja

This does not suck, Maelstorm.

Thanks.

Thanks guys!

katjav, I believe the short answer is yes. The long answer is going to require more explanation than the tutorial above, so, fuck it, I'll just go ahead and make a bit of a PART II here, and I'll use your filter as an example at the end.

NOTE: A .pdf version with much more readable equations is attached below.

To explain this filter as a purely real filter (not an implementation of a complex filter) would require a longer explanation of the z-transform, as the simple approach I gave above doesn't really work in this case. The reason being is that the filter needs to be described using two difference equations (because of that one-sample feedback loop in the imaginary section that is within the overall two-sample loop). For the sake of clarity, I'm going to make these substitutions/names:

c_1 = r*cos(&#969;)
c_2 = r*sin(&#969;)
y_1[n] = real output
y_2[n] = imaginary output

The equations become:

y_1[n] = x[n] + c_1*y_1[n-1] - c_2*y_2[n-1]
y_2[n] = c_2*y_1[n] + c_1*y_2[n-1]

We can see one sample delays here, but they're in two equations. If we find the transfer function that gives us y_1[n], we can see exactly the order of that specific section, and even find a way to implement it without the need for y_2[n]. However, the z-transform as explained earlier is just a shortcut for when you have a simple single-difference-equation filter. It doesn't cut it here. You have to do it the long way.

THE Z-TRANSFORM-extended

When we are in the z-domain, z) is the z-transform of the input, and Y(z) is the z-transform of the output (z itself can be an arbitrary complex number that you can plug in to find useful information about the system, but we don't care here. It's just a variable). Since we are assuming an arbitrary input, we don't know what either of these are, nor to we really care about that either. What we are trying to find is what causes z) to become Y(z). That is H(z), the transfer function of the filter. In the z-domain, the basic formula looks like this:

Y(z) = H(z)*X(z)

In other words, the z-transform of the input times the z-transform of the filter function gives us the z-transform of the output. So, of course, a little algebra gives the basic formula of the transfer function:

Y(z)
H(z) = -----
z)

Remember that. This is the goal when implementing the z-transform.

In the time domain, some thing like x[n-k] is a mathematical way of saying "the input delayed by k samples". In the z-domain, the same thing is said as (z^-k)*X(z). z^-k applies a delay of k samples, and z) is the z-transformed signal that the delay is applied to. This is really all the z-transform itself is. Once we do that, it's just simple algebra to get to H(z).

To illustrate, let's do a simple example:

y[n] = b_0*x[n] + b_1*x[n-1] + a_1*y[n-1]

First, we do the transformation:

Y(z) = b_0*X(z) + b_1*X(z)*z^-1 + a_1*Y(z)*z^-1

Again, we're trying to find H(z) = Y(z)/X(z) here, so we need to rearrange it to get Y(z)/X(z) on one side.

Y(z) - a_1*Y(z)*z^-1 = b_0*X(z) + b_1*X(z)*z^-1

Y(z) * (1 - a_1*z^-1) = z) * (b_0 + b_1*z^-1)

Y(z) b_0 + b_1*z^-1
----- = ------------------- = H(z)
z) 1 - a_1*z^-1

And there you have it. As you can see, that lines up with the simple explanation given earlier, and also explains why the y[n] coefficients change sign while the x[n] ones don't.

Hopefully you can also see why it is so useful for manipulating compared to difference equations. In the time domain, the delay and the signal are linked. y[n] and y[n-1] are separate things. However, in the z-domain, Y(z) and Y(z)*z^-1 allow you to extract Y(z). You can pull the signal away from the delay!

APPLYING TO MULTIPLE DIFFERENCE EQUATIONS

So let's see how to use our newfound knowledge of the z-transform to find the transfer function of a filter that requires multiple difference equations using the one katjav provided. I'll copy the equations again here so you don't have to scroll up:

y_1[n] = x[n] + c_1*y_1[n-1] - c_2*y_2[n-1]
y_2[n] = c_2*y_1[n] + c_1*y_2[n-1]

We want to know the transfer function that gives us the output y_1[n]. First things first, let's transform and roll out! (You knew a Transformers reference was coming at some point, let's be real with ourselves.)

Y_1(z) = z) +c_1*Y_1(z)*z^-1 - c_2*Y_2(z)*z^-1
Y_2(z) = c_2*Y_1(z) + c_1*Y_2(z)*z^-1

We want to find the H(z) for the output that gives us y_1[n], so we need to get to Y_1(z)/X(z). But first, we need to isolate Y_2(z) in the second equation:

Y_2(z) = c_2*Y_1(z) + c_1*Y_2(z)*z^-1

Y_2(z) - c_1*Y_2(z)*z^-1 = c_2*Y_1(z)

Y_2(z) * (1 - c_1*z^-1) = c_2*Y_1(z)

c_2
Y_2(z) = ---------------- * Y_1(z)
1 - c_1*z^-1

Now that we know Y_2(z), we can substitute it into the first equation and find Y_1(z)/X(z). It gets a little messy here, so I'll try not to skip any not-so-obvious steps, for clarity. Hopefully you can follow. Let me know if I screw up somewhere or if something isn't clear:

(c_2*z^-1) * c_2
Y_1(z) = z) + c_1*Y_1(z)*z^-1 - --------------------- * Y_1(z)
1 - c_1*z^-1

(c_2^2)*z^-1
Y_1(z) - c_1*Y_1(z)*z^-1 + ------------------ * Y_1(z) = z)
1 - c_1*z^-1

(c_2^2)*z^-1
Y_1(z) * (1 - c_1*z^-1 + ----------------- ) = z)
1 - c_1*z^-1

Y_1(z) 1
-------- = ---------------------------------------
z) (c_2^2)*z^-1
( 1 - c_1*z^-1 + ----------------- )
1 - c_1*z^-1

Y_1(z) 1
-------- = --------------------------------------------------------------------------
z) 1 - c_1*z^-1 (1 - c_1*z^-1)*c_1*z^-1 (c_2^2)*z^-1
( ---------------- - ----------------------------- + ---------------- )
1 - c_1*z^-1 1 - c_1*z^-1 1 - c_1*z^-1

Y_1(z) 1 - c_1*z^-1
-------- = ---------------------------------------------------------------------
z) 1 - c_1*z^-1 - c_1*z^-1 + (c_1^2)*z^-2 + (c_2^2)*z^-1

Y_1(z) 1 - c_1*z^-1
-------- = ---------------------------------------------------- = H_1(z)
z) 1 + (c_2^2 - 2*c_1)*z^-1 + (c_2^2)*z^-2

And there you have it. Assuming I did that right, that's a second-order filter with two poles and one zero. At least, it is in the z-domain. That's why I said "I believe the short answer is yes" in the beginning and then made you read through all of that. ;-p

http://www.pdpatchrepo.info/hurleur/filter_tut_2.pdf

Puzzling on your equations, still overwhelmed by it. One zero and two poles, that seems plausible from the flow graph as well: one minus and two plus'es in the flow of the real part. The real output gives a band filter. With a similar set of equations, the H(z) could be found for the imaginary output, then leaving out y_1[n], it gives a lo-pass filter. To find the poles, would still require the quadratic equation trick, no?

I still have to understand how H(z)=Y(z)/X(z) from part II relates to H(z)=B(z)/A(z)from part I. Both seem familiar. B(z) is the feed-forward part, from the input samples x, and 1/A(z) is the feedback or IIR part, from the output samples y. Y(z) is the output transform and z) is the input transform. If the output is the product of the input and the filter, the input could be retrieved from the output divided by the filter. Vice versa, the filter could be found from the output divided by the input. So far, it makes sense. But now this: Y(z)/X(z)=B(z)/A(z)? But B(z) comes from the input samples, and A(z) from the output! Seems like nominator and denominator swapped. Why do I suddenly fail to understand things I've always taken for granted?

Is it erh... because of the negative powers of z? z^-1 = 1/z. The negative powers go clockwise on the complex plane, but the signal forward direction is anti-clockwise. That's where the inversion comes from, no? I'm always so confused when it comes to delay.

Maelstorm, where did you learn about z-transform? When I was at school during one year (Sonology Institute Den Haag) our dsp teacher used to say 'z-transform is not interesting for you guys, it is for engineers' (Sonology is for artists). But for me, the complex plane is a place of imagination, and it's secrets play a silent music in my head.

Well, do not spend more time on my questions and confusion, it's endless.

Katja

To be honest, when it comes to filters that are defined by more than one difference equation, I'm a bit in "self-taught" territory, which is why I sort of covered myself with the "I believe" language. I'm pretty sure it's right, but I haven't had it explained to me by someone I'm certain is an authority on the subject. I'm just slightly hesitant because, in simple one-difference-equation filters, a pole corresponds to a feedback delay and a zero corresponds to a feedforward delay, and that's an area I'm much more comfortable with. But I was just playing around today with plotting a state variable filter's response using the same approach, and it looks right. So I guess it makes sense.

H(z) is simply "generally defined" as Y(z)/X(z) in a universal way. H(z) itself is the z-transform of the impulse response. In the time-domain, x[n] convolved with the impulse response h[n] will give you y[n]. Convolution in the time domain translates to multiplication in the z domain. This corresponds to what I mentioned earlier when I said that multiplication in the z domain is the same thing as filters-in-series in the time domain. ANY filter can be represented in the time domain as its impulse response. Take an arbitrary input and convolve it with an IR, and you get the output. Convolve it with another IR, and you get the next output.

H(z) = B(z)/A(z) is simply a mathematical representation of the impulse response in the z domain using the polynomials B(z) and A(z).

In the time domain, if you take a finite number of output samples and de-convolve them with the same number of input samples, you'll have an impulse response of the same length. Likewise, in the z domain, if you take the z-transform of a finite number of output samples and divide them by the z-transform of the same number of input samples, you'll have the z-transform of an impulse response of the same length. In this way, you can generalize an arbitrary filter using the impulse response. But you'll only get an FIR filter this way (unless you use some trickery to estimate the finite H(z)). For an IIR filter, an impulse response is impractical because it is of infinite length. However, it can still be defined as an equation, and that's what B(z)/A(z) does (in the z domain).

Anyway, I may have given more info than needed there to explain, but to see how Y(z)/X(z) = B(z)/A(z), take a look at this step again from part II:

Y(z) * (1 - a_1*z^-1) = z) * (b_0 + b_1*z^-1)

Y(z) b_0 + b_1*z^-1
----- = ------------------- = H(z)
z) 1 - a_1*z^-1

Again, H(z) is generally defined as Y(z)/X(z) for an arbitrary filter. B(z)/A(z) is the function that explicitly computes H(z).

@katjav said:

Maelstorm, where did you learn about z-transform?

I learned about it a little over a year ago while attempting a Master's degree at McGill University. I came from a musician background, but the Music Tech department there is much more engineering-oriented. A large portion of what I know about filters came from my first semester there in an undergraduate-level course. Biggest mind-fuck of a class I've ever taken.

@Maelstorm said:

I'm pretty sure it's right, but I haven't had it explained to me by someone I'm certain is an authority on the subject.

Well, watching a teacher writing formulae on the blackboard is only a starting point for learning. In my view, solid knowledge doesn't come from reading a book or listening to a teacher, but from building things and verifying that they work. And that is what you have done, with your filter tools for Pd.

By the way, I once did a MaxMsp patch with a non-realtime model of a complex integrator (like [cpole~] is) excited by an impuls, acting as a resonator. I've today ported this patch to Pd, see attached.

Katja

http://www.pdpatchrepo.info/hurleur/cpole-model.pd