## Evaluating filter frequency response

A question that pops up for many DSP-ers working with IIR and FIR filters, I think, is how to look at a filter’s frequency and phase response. For many, maybe they’ve calculated filter coefficients with something like the biquad calculator on this site, or maybe they’ve used a MATLAB, Octave, Python (with the scipy library) and functions like freqz to compute and plot responses. But what if you want to code your own, perhaps to plot within a plugin written in c++?

You can find methods of calculating biquads, for instance, but here we’ll discuss a general solution. Fortunately, the general solution is easier to understand than starting with an equation that may have been optimized for a specific task, such as plotting biquad response.

### Plotting an impulse response

One way we could approach it is to plot the impulse response of the filter. That works for any linear, time-invariant process, and a fixed filter qualifies. One problem is that we don’t know how long the impulse response might be, for an arbitrary filter. IIR (Infinite Impulse Response) filters can have a very long impulse response, as the name implies. We can feed a 1.0 sample followed by 0.0 samples to obtain the impulse response of the filter. While we don’t know how long it will be, we could take a long impulse response, perhaps windowing it, use an FFT to convert it to the frequency domain, and get a pretty good picture. But it’s not perfect.

For an FIR (Finite Impulse Response) filter, though, the results are precise. And the impulse response is equal to the coefficients themselves. So:

For the FIR, we simply run the coefficients through an FFT, and take the absolute value of the complex result to get the magnitude response.

(The FFT requires a power-of-2 length, so we’d need to append zeros to fill, or use a DFT. But we probably want to append zeros anyway, to get more frequency points out for our graph.)

### Plotting the filter precisely

Let’s look for a more precise way to plot an arbitrary filter’s response, which might be IIR. Fortunately, if we have the filter coefficients, we have everything we need, because we have the filter’s transfer function, from which we can calculate a response for any frequency.

The transfer function of an IIR filter is given by

### $$H(z)=\frac{a_{0}z^{0}+a_{1}z^{-1}+a_{2}z^{-2}…}{b_{0}z^{0}+b_{1}z^{-1}+b_{2}z^{-2}…}$$

z0 is 1, of course, as is any value raised to the power of 0. And for normalized biquads, b0 is always 1, but I’ll leave it here for generality—you’ll see why soon.

To translate that to an analog response, we substitute e for z, where ω is 2π*freq, with freq being the normalized frequency, or frequency/samplerate:

### $$H(e^{j\omega})=\frac{a_{0}e^{0j\omega}+a_{1}e^{-1j\omega}+a_{2}e^{-2j\omega}…}{b_{0}e^{0j\omega}+b_{1}e^{-1j\omega}+b_{2}e^{-2j\omega}…}$$

Again, e0jω is simply 1.0, but left so you can see the pattern. Here it is restated using summations of an arbitrary number of poles and zeros:

### $$H(e^{j\omega})=\frac{\sum_{n=0}^{N}a_{n}e^{-nj\omega}}{\sum_{m=0}^{M}b_{m}e^{-mj\omega}}$$

For any angular frequency, ω, we can solve H(e). A normalized frequency of 0.5 is half the sample rate, so we probably want to step it from 0 to 0.5—ω from 0 to π—for however many points we want to evaluate and plot.

### Coding it

From that last equation, we can see that a single FOR loop will handle the top or the bottom coefficient sets. Here, we’ll code that into a function that can evaluate either zeros (a terms) or poles (b terms). We’ll refer to this as our direct evaluation function, since it evaluates the coefficients directly (as opposed to evaluating an impulse response).

You’ve probably noticed the j, meaning an imaginary part of a complex number—the output will be complex. That’s OK, the output of an FFT is complex too, and we know how to get magnitude and phase from it already.

Some languages support complex arithmetic, and have no problem evaluating “e**(-2*j*0.5)”—either directly, or with an “exp” (exponential) function. It’s pretty easy in Python, for instance. (Something like, coef[idx] * math.e**(-idx * w * 1j), as the variable idx steps through the coefficients array.)

For languages that don’t, we can use Euler’s formula, ejx = cos(x) + j * sin(x); that is, the real part is the cosine of the argument, and the imaginary part is the sine of it.

(Remember, j is the same as i—electrical engineers already used i to symbolize current, so they diverged from physicists and used j. Computer programming often use j, maybe because i is a commonly used index variable.)

So, we create our function, run it on the numerator coefficients for a given frequency, run it again on the denominator coefficients, and divide the two. The result will be complex—taking the absolute value gives us the magnitude response at that frequency.

### Revisiting the FIR

Since we already had a precise method of looking at FIR response via the FFT/DFT, let’s compare the two methods to see how similar they are.

To use our new method for the case of an FIR, we note that the denominator is simply 1, so there is no denominator to evaluate, no need for division. So:

For the FIR, we simply run the coefficients through our evaluation function, and take the absolute value of the complex result to get the magnitude response.

Does that sound familiar? It’s the same process we outlined using the FFT.

### And back to IIR

OK, we just showed that our new evaluation function and the FFT are equivalent. (There is a difference—our evaluation function can check the response at an arbitrary frequency, whereas the FFT frequency spacing is defined by the FFT size, but we’ll set that aside for the moment. For a given frequency, the two produce identical results.)

Now, if the direct evaluation function and the FFT give the same results, for the same frequency point, and the numerator and denominator are evaluated by the same function, by extension we could also get a precise evaluation by substituting an FFT process for both the numerator and denominator, and dividing the two as before. Note that we’re no longer talking about the FFT of the impulse response, but the coefficients themselves. That means we no longer have the problem of getting the response of an impulse that can ring out for an unknown time—we have a known number of coefficients to run through the FFT.

### Which is better?

In general, the answer is our direct evaluation method. Why? We can decide exactly where we want to evaluate each point. That means that we can just as easily plot with log frequency as we can linear.

But, there may be times that the FFT is more suitable—it is extremely efficient for power-of-2 lengths. (And don’t forget that we can use a real FFT—the upper half of the general FFT results would mirror the lower half and not be needed.)

### An implementation

We probably want to evaluate ω from 0 to π, corresponding to a range of half the sample rate. So, we’d call the evaluation function with the numerator coefficients and with the denominator coefficients, for every ω that we want to know (spacing can be linear or log), and divide the two. For frequency response, we’d take the absolute value (equivalently, the square root of the sum of the squared real and imaginary parts) of each complex result to obtain magnitude, and arc tangent of the imaginary part divided by the real part (specifically, we use the atan2 function, which takes into account quadrants). Note that this is the same conversion we use for FFT results, as you can see in my article, A gentle introduction to the FFT.

$$magnitude:=\left |H \right |=abs(H)=\sqrt{H.real^2+H.imag^2}$$ $$phase := atan2(H.imag,H.real)$$

For now, I’ll leave you with some Python code, as it’s cleaner and leaner than a C or C++ implementation. It will make it easier to transfer to any language you might want (Python can be quite compact and elegant—I’m going for easy to understand and translate with this code). Here’s the direct evaluation routine corresponding to the summation part of the equation (you’ll also need to “import numpy” to have e available—also available in the math library, but we’ll use numpy later, so we’ll stick with numpy alone):

import numpy as np

# direct evaluation of coefficients at a given angular frequency
def coefsEval(coefs, w):
res = 0
idx = 0
for x in coefs:
res += x * np.e**(-idx * 1j * w)
idx += 1
return res

Again, we call this with the coefficients for each frequency of interest. Once for the numerator coefficients (the a coefficients on this website, corresponding to zeros), once for the denominator coefficients (b, for the poles—and don’t forget that if there is no b0, the case for a normalized filter, insert a 1.0 in its place). Divide the first result by the second. Use use abs (or equivalent) for magnitude and atan2 for phase on the result. Repeat for every frequency of interest.

Here’s a python function that evaluates numerator and denominator coefficients at an arbitrary number of points from 0 to π radians, with linear spacing, returning array of magnitude (in dB) and phase (in radian, between +/- π):

# filter response, evaluated at numPoints from 0-pi, inclusive
def filterEval(zeros, poles, numPoints):
magdB = np.empty(0)
phase = np.empty(0)
for jdx in range(0, numPoints):
w = jdx * math.pi / (numPoints - 1)
resZeros = coefsEval(zeros, w)
resPoles = coefsEval(poles, w)

# output magnitude in dB, phase in radians
Hw = resZeros / resPoles
mag = abs(Hw)
if mag == 0:
mag = 0.0000000001  # limit to -200 dB for log
magdB = np.append(magdB, 20 * np.log10(mag))
phase = np.append(phase, math.atan2(Hw.imag, Hw.real))
return (magdB, phase)

Here’s an example of evaluating biquad coefficients at 64 evenly spaced frequencies from 0 Hz to half the sample rate (these coefficients are right out of the biquad calculator on this website—don’t forget to include b0 = 1.0):

zeros = [ 0.2513643668578741, 0.5027287337157482, 0.2513643668578741 ]
poles = [ 1.0, -0.17123074520885395, 0.1766882126403502 ]

(magdB, phase) = filterEval(zeros, poles, 64)

print("\nMagnitude:\n")
for x in magdB:
print(x)

print("\nPhase:\n")
for x in phase:
print(x)

Next up, a javascript widget to plot magnitude and phase of arbitrary filter coefficients.

### Extra credit

The direct evaluation function performs a Fourier analysis at a frequency of interest. For better understanding, reconcile it with the discrete Fourier transform described in A gentle introduction to the FFT. In that article, I describe probing the signal with cosine and sine waves to obtain the response at a given frequency. Look again at Euler’s formula, which shows that e is cosine (real part) and sine (imaginary part), which the article alludes to this under the section “Getting complex”. You should understand that the direct evaluation function presented here could be used to produce a DFT (given complete evaluation of the signals at appropriately spaced frequencies). The main difference is that for this analysis, we need not do a complete and reversible transform—we need only analyze frequency response values that we want to graph.

Posted in Biquads, FFT, Filters, FIR Filters, IIR Filters | 18 Comments

Sometimes we’d like to cascade biquads to get a higher filter order. This calculator gives the Q values for each filter to achieve Butterworth response for lowpass and highpass filters.

Order:
Q values:

You can calculate coefficients for all biquad (and one-pole) filters with the biquad calculator.

Sometimes we’d like a steeper cutoff than a biquad—a second order filter—gives us. We could design a higher order filter directly, but the direct forms suffer from numerical problems due to limited computational precision. So, we typically combine one- and two-pole (biquad) filters to get the order we need. The lower order filters are less sensitive to precision errors. And we maintain the same number of math operations and delay elements as the equivalent higher order filter, so think of cascading as simply rearranging the math.

The main problem with cascading is that if you take two Buterworth filters in cascade, the result is no longer Butterworth. Consider a Butterworth—maximally flat passband—lowpass filter. At the defined corner frequency, the magnitude response is -3 dB. If you cascade two of these filter, the response is now -6 dB. We can’t simply move the frequency up to compensate, since the slope into the corner is also not as sharp. Increasing the Q of both filters to sharpen the corner would degrade the passband’s flatness. We need a combination of Q values to get the correct Butterworth response.

### How to calculate Q values

The problem of figuring out what the Q should be for each stage of a biquad cascade becomes very simple if we look at the pole positions of the Butterworth filter we want to achieve in the s-plane. In the s-plane, the poles of a Butterworth filter are distributed evenly, at a constant radius from the origin and with a constant angular spacing. Since the radius corresponds to frequency, and the pole angle corresponds to Q, we know that all of the component filters should be set to the same frequency, and their Q is simple to calculate from the pole angles. For a given pole angle, θ, Q is 1 / (2cos(θ)).

Calculating pole positions is easy: For a filter of order n, poles are spaced by an angle of π/n. For an odd order, we’ll have a one-pole filter on the real (horizontal) axis, and the remaining pole pairs spaced at the calculated angle. For even orders, the poles will be mirrored about the real axis, so the first pole pairs will start at plus and minus half the calculated angle. The biquad poles are conjugate pairs, corresponding to a single biquad, so we need only pay attention to the positive half for Q values.

### Examples

For a 2-pole filter, a single biquad, the poles are π/2 radians apart, mirrored on both sides of the horizontal axis. So, our single Q value is based on the angle π/4; 1/(2cos(π/4)) equals a Q value of 0.7071.

For a 3-pole filter, the pole spacing angle is π/3 radians. We start with a one-pole filter on the real (σ) axis, so the biquad’s pole angle is π/3; 1/(2cos(π/3)) equals a Q of 1.0.

For a 4-pole filter, we have two biquads, with poles spaced π/4 radians apart, mirrored about the real axis. That means the first biquad’s pole angle is π/8, and the second is 3π/8, yielding Q values of 0.541196 and 1.306563.

## Filters for synths—the 4-pole

The last post noted that the two most popular synthesizer filters are the 2-pole state variable, and the 4-pole “Moog style”. And we started with the state variable—simple, popular, and delivering multiple filter outputs (lowpass, bandpass…) simultaneously. Here, we’ll follow up with comments on the filter associated with Moog (and especially the Minimoog). In general, we’ll refer to this as a 4-pole synth filter.

While this filter is usually thought of as a lowpass filter, the other popular filter types can be derived easily. Many people first saw this in the Oberheim Xpander (and Matrix-12) synths of the ’80s, the idea came from Bernie Hutchins’ Electronotes in the ’70s. So don’t feel that you must go the direction of the state variable is you want multiple filter types, including 2-pole response.

Lowpass response is the overwhelming choice for typical synth use. Note that a 4-pole lowpass is not necessarily better then a 2-pole (as in the state variable)—they are just choices. You might want a 4-pole for the darker Minimoog bass sounds, and a 2-pole for the brassy OB8-style sounds.

### Basic construction

The 4-pole is implemented by a string of four one-pole lowpass filter in series. We need corner peaking and resonance control for a synth filter, and we get that by feeding back the output to the input. While trivial in the analog domain, this feedback is the tricky part in the digital recreations. The reason is that it’s not a continuous system, and the obvious way to handle it is to put a delay in that part, so the output of the current sample period is available as input for the next. But this creates some bad side effects, particularly for tuning. In the past, people dealt with this by accounting for those errors.

### Approaching analog

But it’s not just tuning errors—if it were, that would be simple to fix. The Minimoog popularity, in part, is that it is designed to easily overdrive the filter, to get a “fat” tone. This is another thing that is simple in the analog domain, but doing the same in the digital domain produces noticeably digital artifacts. And if your goal is to make something that sound analog, this is a source of spectacular “fail”.

So instead of this simple delay approach in the 4-pole feedback path, modern ideas use more complex techniques to avoid the large-scale errors in an effort to get closer to how the analog counterpart works. And part of the effort is in dealing with an overdriven feedback path. The result reduced digital artifacts, makes the filter’s behavior more closely resemble its analog counterpart when overdriven, and also gives a smoother, more predictable and more musical sound at high resonance.

Note that these techniques are often called “zero feedback delay” (ZDF) filters. That is meant to highlight the fact a trivial delay is not used. I’m not a huge fan of that, since it’s not meaningful to someone who doesn’t know of the delay it refers to, and of course there are always sources of internal delay in an such filter design. But I mention ZDF so that if you’ve heard it before, be assured that we are talking about those sort of techniques here.

A great resource for this topic is Vadim Zavalishin’s The Art of VA Filter Design (“VA” for “Virtual Analog”).

Posted in Digital Audio, Filters, IIR Filters, Synthesizers | 5 Comments

## Filters for synths–starting out

We haven’t developed a synth filter here yet…

### Filters we’ve presented

Biquads. While they are useful for many simple cases of filtering, they are not a good choice for analog synthesizer emulation. Most notably, they are poorly suited to time-varying parameters such as smooth filter sweeps.

Hal Chamberlin’s digital state variable filter, which has the advantage of independent control of filter frequency and resonance (Q), as well as simultaneous output of lowpass, bandpass, highpass, and notch filtering. While this served in many synthesizers over the years, using oversampling to fix its most glaring problem (limited frequency control range before losing stability), more modern design approaches offer far better performance.

### What analog synths use

It’s helpful to discuss what most people are used to before getting to digital implementations.

For analog synthesizer emulation, two basic analog designs stand out, having withstood the test of time the best. The 24 dB/octave (“four pole”) Moog style (“transistor ladder”) lowpass filter, and the 12 dB/octave (“two pole”) state variable with its multiple outputs.

There are other designs, including the “diode ladder” design by Roland, the Steiner-Parker filter (basically Sallen-Key)—these are some of the ways voltage control was solved in analog filters, resulting in characteristic sounds. We’ll limit discussion to the two standouts.

### Where to start

The state-variable filter is the best place to start—simple and versatile.

And my choice for your first stop, for a simple filter that’s scarcely more complex than a Chamberlin svf is Andrew Simper’s trapezoidally integrated svf. There is enough detail and code in his pdf for anyone to implement a useable, flexible, and good sounding synthesizer filter—also useable for other filtering applications such as equalization.

http://www.cytomic.com/files/dsp/SvfLinearTrapOptimised2.pdf

Dig into that for a start! (And further reading from Andy)

### Bonus audio sample

Here is the sound of two of our sawtooth wavetable oscillators, detuned, through the SVF coded directly from Andy’s paper, with frequency swept from 20 kHz down through the audio range by our ADSR; note that such a sweep would blow up the Chamberlin implementation without oversampling:

SVF lowpass sweep from 20 kHz

## Audio Signal Processing for Music Applications

I’d like to recommend this excellent—and free—online course:

Audio Signal Processing for Music Applications
by Prof Xavier Serra, Prof Julius O Smith, III

The brief: In this course you will learn about audio signal processing methodologies that are specific for music and of use in real applications. You will learn to analyse, synthesize and transform sounds using the Python programming language.

This is the second session of the course, which will start on September 21, 2015. Enrollment is open here:

https://www.coursera.org/course/audio

I took the first session last fall. My primary motivation was to try a MOOC (Massive Open Online Course), and I wanted to pick a topic I already had knowledge in, to better evaluate the experience. For me, this topic was one in which I had an idea of the basic techniques used, but not the details and no practical experience, so it was a good fit. Plus, the course requires the use of Python, which I had interest in and wanted an excuse to learn.

Professor Serra delivers the course videos, and does an excellent job—clear and well-paced. The course helps to give a better understanding of musical components of sound, and the techniques used to alter individual aspects (especially pitch and duration, independently).

The ten-week course is not easy, and requires at least a few hours per week—a relatively small percentage of the large number of enrollees completed the course with a passing grade (yes, I passed). It requires a relatively modest amount of programming, as the sms-tools package handles most of the work, but will be tough for non-programmers. However, you can watch the videos even if you do not wish to do the weekly assignments and quizzes. The assignments will give you a much better understanding, of course.

Posted in Uncategorized | 1 Comment

## Dither—The Naked Truth video

This video presents the “naked truth” on dither and truncation error, by stripping away the original signal of a musical clip and listening at different bit levels. I boost the error to a normalized audio volume for easy comparison of sound quality between different sample sizes, so your listening environment is not critical, but headphones will be a plus.

Posted in Digital Audio, Dither, Video | 2 Comments

## Dither widget

 Sine Run Rounded Connect Dither Added Dither + Rounded

This is the widget I used for the animated dither chart in my Audio Dither Explained video. “Run” animates the dither process, otherwise it changes only when you change the amplitude setting. “Connect” connects the dots; this is simply a visual aid, making it easily to follow the sequence of samples.

## Audio Dither Explained video

This video discusses the how and why of dither.

Posted in Digital Audio, Dither, Video | 1 Comment

I was admonished (in not a nice way—not terribly rude, but quite to the point that I don’t know what I’m doing) by an anonymous visitor, who concluded that I don’t know much about C++ and maybe should have written it in a different language (why, when the purpose is to give people C++ code?). I would reply to him if he had left a real email address, but maybe I should make a few things clear to all.

I supply bare bones code. I try to leave out anything not necessary, and avoid even fairly obvious optimizations (and some less obvious that would make the code more difficult to understand. For instance, the wavetables could be one sample larger and let you avoid testing for wrap-around and the resulting branch). See About source code examples.

There is a common misconception that if your class is suitable for subclassing, then the destructor must be declared virtual. This is not true. If it were, the C++ architects would have simply built it in.

So, to some (including the anonymous commenter), it may seem confusing that I don’t declare destructors as virtual, which they might conclude to imply no subclassing, while declaring variables as protected (not private), implying support for subclassing. There is a method to this apparent madness: First, I don’t expect that the object will be subclassed by you, and I want to keep it efficient. So, I don’t declare destructors virtual, which would create some overhead we don’t need. However, if you do choose to subclass, then it’s important that you have access to certain variables, hence the protected declarations where appropriate.

However, that fact is that base class destructors need not be virtual when subclassed—that’s a requirement only if you want to delete them via the base class. By not declaring them virtual, we get the lack of vtable in the typically case (of no subclasses), while still allowing subclassing.

Also, I don’t mess with namespaces, and avoid templates even though they would be useful—again, trying to focus on the code that implements the DSP algorithm. I try to fit as much of the pertinent code on the screen for you as possible, so while I like comments in my code, I keep them brief.

And I don’t go through and declare every non-changing variable as const. This is not open source code. I’m not trying to protect you or some group you are working on code with from yourselves. I’m giving you seeds, not a framework, and it’s not in a Git repository. If I keep it simple enough, I figure you’ll be able to pull the bits out that matter for you into your own code. I don’t have time for both open source and tutorials—you can find a ton of open source oscillators, filters, and modulators on the internet, and I think you’ll find that the authors do not teach you how and why the code works. The latter is few and far between, and that’s why I’m here.

I’m open to suggestions, but realize that I do think about my choices, even if it’s not apparent. So, please, let’s make it a discussion. And the anonymous fellow is invited to point us to his contributions to the public. 😉

Posted in Source Code | 7 Comments

## Pole-Zero placement v2

 Pair Pole mag Pole angle Pair Zero mag Zero angle Sample rate (Hz) linearlog Plot

### A new pole-zero calculator

An JavaScript remake of the old Java-based pole-zero placement applet—visit that page for tips on pole-zero locations for standard biquads. The main additions are input fields for precision pole-zero placement, and an option to display the response with a log frequency scale.

The basic idea is that poles blow, zeros suck. Think of poles as controlling a frequency-dependent feedback or resonance—the impulse response of a pole inside the unit circle decays, while one outside is like runaway feedback (think of a mic feeding back into a loudspeaker). A pole on the unit circle gives a sustained oscillation (but watch out for numerical errors—keep your poles inside the unit circle, typically). Zeros absorb a particular frequency; when on the unit circle, they absorb the corresponding frequency completely.

So, poles push the frequency response up around their corresponding frequency, and zeros pull down around theirs. Keep in mind that the frequency response graph is normalized, just as the filter coefficients are. So, while a pole pushes up the response, it appears as though all other frequencies are being pushed down instead. Of course, normalization is important in practical application, but be aware of it when visualizing how poles and zeros interact.

Posted in Biquads, Digital Audio, Filters, IIR Filters, Widgets | 13 Comments