I don’t want to get into the business of teaching people how to code—there are a huge number of free resources available on the internet to that do that. But I’ll give a small taste for those trying to get started.

First, I want to be clear that biquads, especially the direct forms, get a lot of coverage simply because they are a straight-forward and basic building block—not because they are always the best choice in filters. The direct forms are terrible choices if you want a sweepable synthesizer filter, for instance. But they are often an adequate and computationally efficient choice for many types of audio and control signal filtering. And they are simple and easy to understand. And since beginners need a place to start, they serve as a good example for coding a DSP example.

C++ is a popular language for DSP coding, and coding in general. It’s object-oriented, and that makes for a good model in developing our audio processing building blocks. Here’s how we approach a biquad filter object.

### Filter object requirements

Our filter object should carry it’s own settings and states. That way, we can have multiple instances of the filter, each remaining independent, with its own filter type, frequency, and internal state. We define the state and settings variables (.h) file, giving them a name and reserving space for them.

For this filter, we’ll have a filter type (lowpass, highpass, bandpass, notch, peak, lowshelf, and highshelf), frequency, Q, and peak gain (used only for the peak and shelf types). I’ll use normalized frequency, where 1.0 is the sample rate—divide the frequency you want to set, in Hertz or cycles per second, by your sample rate to get normalized frequency (2438.0/44100 for 2438 Hz at a sample rate of 44100 Hz). We’ll set peak gain in dB, where negative numbers are for cut, and positive numbers for boost.

Biquads need to convert the frequency, Q, and peak gain parameters to the a0, a1, a2, b1, and b2 factors. Since that takes a fair amount of computation, it makes sense to calculate the factors only when the parameters change, which is almost always much lower than the sample rate at which the filter processes audio. So, we’ll want to reserve space to hold those factors. And because we want to be able to set the parameters individually, and because they rely on each other, we’ll keep copies of the current parameter settings as well.

For this example, we’ll use the transposed direct form II architecture (which has good floating point characteristics, and requires only two unit delays), so we’ll need to reserve two memory locations for the delay elements.

To specify the filter type, we assign a number for each type, using enumerated types so that we can attach a more meaningful name to the type. Here’s a look at that definition:

enum { bq_type_lowpass = 0, bq_type_highpass, ... bq_type_highshelf };

And here are the variables we’re allocating for each filter instance:

int type; double a0, a1, a2, b1, b2; double Fc, Q, peakGain; double z1, z2;

### Filter functions

Next, we need to think about the software functions. We’ll need functions to set the type, frequency, Q, and peak gain individually. For convenience, we’ll also include a function that accepts all of the parameters at once, particularly handy for setting fixed filters. Also, we always need a constructor, called to create a filter instance, and a destructor, called (sometimes implicitly) to remove a filter instance and free its memory. For convenience we’ll include a bare-bones constructor, and a second one that includes setting the filter—especially handy for filters that need to be set just once.

And we’ll need to function that processes and input sample and returned the filtered result—this is called for each sample, and we want it to be as efficient as possible. If our processing routine is small, we’d like to make it an inline function, which saves some overhead of the function call by letting the compiler place the code contained in the function inline. To do that, we define the function as inline, and place the code in the header file. Since all all of these functions need to be accessible by other code, they need to be defined as public; since all of the variables are internal to the filter, we define them as protected—external code has no access. (We make the code a little more tidy by splitting off the computation biquad factors into a separate function; since it’s only accessed internally, we make it protected.)

We’ll use the double type (64-bit floating point) for our computations, and float (32-bit floating point) as the input and output sample size. (We could define 64-bit samples, or use templates to accommodate different sample types, but that’s overdoing it for this example.) Here’s our completed header file (with a common technique to avoid nested includes from being a problem for the compiler):

// // Biquad.h // // Created by Nigel Redmon on 11/24/12 // EarLevel Engineering: earlevel.com // Copyright 2012 Nigel Redmon // // For a complete explanation of the Biquad code: // http://www.earlevel.com/main/2012/11/26/biquad-c-source-code/ // // License: // // This source code is provided as is, without warranty. // You may copy and distribute verbatim copies of this document. // You may modify and use this source code to create binary code // for your own purposes, free or commercial. // #ifndef Biquad_h #define Biquad_h enum { bq_type_lowpass = 0, bq_type_highpass, bq_type_bandpass, bq_type_notch, bq_type_peak, bq_type_lowshelf, bq_type_highshelf }; class Biquad { public: Biquad(); Biquad(int type, double Fc, double Q, double peakGainDB); ~Biquad(); void setType(int type); void setQ(double Q); void setFc(double Fc); void setPeakGain(double peakGainDB); void setBiquad(int type, double Fc, double Q, double peakGainDB); float process(float in); protected: void calcBiquad(void); int type; double a0, a1, a2, b1, b2; double Fc, Q, peakGain; double z1, z2; }; inline float Biquad::process(float in) { double out = in * a0 + z1; z1 = in * a1 + z2 - b1 * out; z2 = in * a2 - b2 * out; return out; } #endif // Biquad_h

Now, for the implementation of the rest of our functions in the .cpp file:

// // Biquad.cpp // // Created by Nigel Redmon on 11/24/12 // EarLevel Engineering: earlevel.com // Copyright 2012 Nigel Redmon // // For a complete explanation of the Biquad code: // http://www.earlevel.com/main/2012/11/26/biquad-c-source-code/ // // License: // // This source code is provided as is, without warranty. // You may copy and distribute verbatim copies of this document. // You may modify and use this source code to create binary code // for your own purposes, free or commercial. // #include <math.h> #include "Biquad.h" Biquad::Biquad() { type = bq_type_lowpass; a0 = 1.0; a1 = a2 = b1 = b2 = 0.0; Fc = 0.50; Q = 0.707; peakGain = 0.0; z1 = z2 = 0.0; } Biquad::Biquad(int type, double Fc, double Q, double peakGainDB) { setBiquad(type, Fc, Q, peakGainDB); z1 = z2 = 0.0; } Biquad::~Biquad() { } void Biquad::setType(int type) { this->type = type; calcBiquad(); } void Biquad::setQ(double Q) { this->Q = Q; calcBiquad(); } void Biquad::setFc(double Fc) { this->Fc = Fc; calcBiquad(); } void Biquad::setPeakGain(double peakGainDB) { this->peakGain = peakGainDB; calcBiquad(); } void Biquad::setBiquad(int type, double Fc, double Q, double peakGainDB) { this->type = type; this->Q = Q; this->Fc = Fc; setPeakGain(peakGainDB); } void Biquad::calcBiquad(void) { double norm; double V = pow(10, fabs(peakGain) / 20.0); double K = tan(M_PI * Fc); switch (this->type) { case bq_type_lowpass: norm = 1 / (1 + K / Q + K * K); a0 = K * K * norm; a1 = 2 * a0; a2 = a0; b1 = 2 * (K * K - 1) * norm; b2 = (1 - K / Q + K * K) * norm; break; case bq_type_highpass: norm = 1 / (1 + K / Q + K * K); a0 = 1 * norm; a1 = -2 * a0; a2 = a0; b1 = 2 * (K * K - 1) * norm; b2 = (1 - K / Q + K * K) * norm; break; case bq_type_bandpass: norm = 1 / (1 + K / Q + K * K); a0 = K / Q * norm; a1 = 0; a2 = -a0; b1 = 2 * (K * K - 1) * norm; b2 = (1 - K / Q + K * K) * norm; break; case bq_type_notch: norm = 1 / (1 + K / Q + K * K); a0 = (1 + K * K) * norm; a1 = 2 * (K * K - 1) * norm; a2 = a0; b1 = a1; b2 = (1 - K / Q + K * K) * norm; break; case bq_type_peak: if (peakGain >= 0) { // boost norm = 1 / (1 + 1/Q * K + K * K); a0 = (1 + V/Q * K + K * K) * norm; a1 = 2 * (K * K - 1) * norm; a2 = (1 - V/Q * K + K * K) * norm; b1 = a1; b2 = (1 - 1/Q * K + K * K) * norm; } else { // cut norm = 1 / (1 + V/Q * K + K * K); a0 = (1 + 1/Q * K + K * K) * norm; a1 = 2 * (K * K - 1) * norm; a2 = (1 - 1/Q * K + K * K) * norm; b1 = a1; b2 = (1 - V/Q * K + K * K) * norm; } break; case bq_type_lowshelf: if (peakGain >= 0) { // boost norm = 1 / (1 + sqrt(2) * K + K * K); a0 = (1 + sqrt(2*V) * K + V * K * K) * norm; a1 = 2 * (V * K * K - 1) * norm; a2 = (1 - sqrt(2*V) * K + V * K * K) * norm; b1 = 2 * (K * K - 1) * norm; b2 = (1 - sqrt(2) * K + K * K) * norm; } else { // cut norm = 1 / (1 + sqrt(2*V) * K + V * K * K); a0 = (1 + sqrt(2) * K + K * K) * norm; a1 = 2 * (K * K - 1) * norm; a2 = (1 - sqrt(2) * K + K * K) * norm; b1 = 2 * (V * K * K - 1) * norm; b2 = (1 - sqrt(2*V) * K + V * K * K) * norm; } break; case bq_type_highshelf: if (peakGain >= 0) { // boost norm = 1 / (1 + sqrt(2) * K + K * K); a0 = (V + sqrt(2*V) * K + K * K) * norm; a1 = 2 * (K * K - V) * norm; a2 = (V - sqrt(2*V) * K + K * K) * norm; b1 = 2 * (K * K - 1) * norm; b2 = (1 - sqrt(2) * K + K * K) * norm; } else { // cut norm = 1 / (V + sqrt(2*V) * K + K * K); a0 = (1 + sqrt(2) * K + K * K) * norm; a1 = 2 * (K * K - 1) * norm; a2 = (1 - sqrt(2) * K + K * K) * norm; b1 = 2 * (K * K - V) * norm; b2 = (V - sqrt(2*V) * K + K * K) * norm; } break; } return; }

Download the Biquad source code.

### Usage examples

Job done. Now how to we use it? When we need a biquad, we just instantiate the Biquad object, request its initial setting, then call the process function for every sample we want to filter.

Biquad *lpFilter = new Biquad(); // create a Biquad, lpFilter; lpFilter->setBiquad(bq_type_lowpass, Fc / sampleRate, 0.707, 0); // filter a buffer of input samples, in-place for (int idx = 0; idx < bufSize; idx++) { in[idx] = lpFilter->process(in[idx]); }

Here’s an example of cascading two peak filters to yield a 3 dB dip at 200 Hz and a 6 dB bump at 6 kHz:

Biquad *filter1 = new Biquad(bq_type_peak, 200.0 / sampleRate, 1.0, -3.0); Biquad *filter2 = new Biquad(bq_type_peak, 6000.0 / sampleRate, 5.0, 6.0); // filter a buffer of input samples, in-place for (int idx = 0; idx < bufSize; idx++) { in[idx] = filter2->process(filter1->process(in[idx])); }

Fantastic post, and fantastic blog. This is an incredibly useful resource for DSP beginners like myself to study.

How well do simple biquad filters like this respond to modulation of cutoff or q values?

Excellent question, Tom, and the answer is “not well”. Plus there is the complication of having to change multiple coefficients. Contrast that, for instance, with the Chamberlin state-variable filter, with independent frequency and q, or similarly with the Moog-style digital filters (four series one-pole filters plus overall feedback). I don’t suggest the Chamberlin SVF as-is—it’s only adequate with oversampling, and there are much-improved versions. I hope to get to the topic of a usable synth filter in the coming months. Here’s a teaser—the wavetable oscillator (sawtooth), through a quick and dirty improved SVF (no oversampling), modulated by an envelope generator that I wrote over the holidays (article in progress).

Saw through SVF, modulated by envelope

Hey Nigel,

Have you got any plans to get this article out ? I’d be really interested in reading through your approach. I’ve just about got my head around VA filter design and implementations of moog ladder filters etc.

Would love to see an ear level take on usable synth filters.

Interesting timing, Josh. I wrote a quick filter article a few days ago (I’ll edit it this weekend and post). Mainly, I point out VA filter design papers from other people, so it won’t give you any new information, though. Mainly, I wanted to be sure that readers understand that there are much better alternatives for modulatable synthesizer filters than posted here so far (liquids have their place, but a synth filter isn’t one of them, and there are much improved state variable implementations over the Chamberlin).

Thanks for the source code, very informative and very helpful!

Hi Nigel

I read the above code and your codes also runs relativly smoothly.

I have a more basic Q – since I havn’t done much filter design before…

I would like to generate Low/band/high pass filter using your code

could you please elaborate or refer me to information about what is the meaning of following input parameter in the constructor

Fc ? (I am guessing this is cutoff, but do I sepcify Fpass )

Q ? how I can I determine it based on the shape of filter which I do want to have

peakGainDB ?

Thanks a lot

Yehuda

Hi Yehuda,

You’re right—I should have included some notes on how to use the code. The best way to see which parameters are used for a given filter shape it to try it out on the biquad calculator. The one difference is that Fc, in the code, is “normalized” cutoff or center frequency. That is, a value of 0.25 corresponding to one-fourth the sample rate. So, to set a frequency of 1000 Hz when running at 44.1 kHz sample rate, use 1000.0 / 44100.0.

A you’ll see in the calculator, setPeakGain affects only the low and high shelf and the peak filter. The setQ function affects all but the shelving filters. You can set these for any filter, but the values aren’t used by the filters that don’t need them.

Nigel

Thanks 😉

Thanks for posting! This article has been extremely helpful.

I have one question though:

I need a low pass with a much steeper curve for what I am doing, and based on the Biquad Calculator v2, I actually switched to using a negative high shelf instead. It’s still not steep enough though.

So if I wanted to create a fourth-order or higher filter, would this be the correct way to go about doing it?

inline float Biquad::process(float in) {

double out = in;

int order = 4;

for(int i = 0; i < order; ++i)

out = out * a0 + z1;

z1 = in * a1 + z2 – b1 * out;

z2 = in * a2 – b2 * out;

return out;

}

Thanks again

Hi Nicholas,

That would not work, because each second order section needs its own set of delays. Instead of modifying the function for higher orders, simply create two (or more) of the liquids, and pass the output of one into the next. Something like:

Biquad *filter1 = new Biquad();

Biquad *filter2 = new Biquad();

… set them…

// Process

out = filter2->process(filter1->process(input));

Hi Nigel,

Thanks for this implementation, am using it in my Synthedit BiQuad filter module.

Am sharing the source code, included yours.

Cheers!

Great! Thanks for the heads-up, and the link to your website.

How would you handle invalid Fc or Q? Should the filter implementation handle these and return some sane coefficients (instead of NaN for example), or just leave it to whomever is using the filter.

First, the user interface should handle it to some degree. That is, if you’re running this from real (in embedded hardware) or virtual (a plug-in or app) controls, the control ranges should be valid. That doesn’t limit everything, though, since you could sum an LFO or envelope generator with a setting (say, filter cutoff), and get an out-of-range parameter. You could have limits checking and return error codes, but that would be silly, requiring a lot of error testing and decisions in callers. The parameter accessors, such as setFc, should simply clip the parameters appropriately–limiting Fc to between zero and 0.5 (half the sample rate), for instance. (if (Fc > 0.5) Fc = 0.5;…)

Hi Nigel,

great post, very nice code and useful Biquad calculator.

I have the question. If I would like to create a cascade of filters, is there any way of vizualization the resulting curve?

Thanks

Niro

A cascade to implement higher orders, such as fourth, sixth? Basically, you need to have each biquad implement two poles and two zeros each…I might write about that at some point, but I have more interesting topics to get to first…

How would i go about doing an equalizer with this , im new to c and finding it difficult to implement it,

what i am trying to do is make a 10 band graphic equalizer to reduce or enhance certain frequencies inside a sample.

You would create multiple peaking filters (you could also put a low shelf and a high shelf at the extremes if you want). Your control logic would set the peak gain of each according to your controls. Then at process time you would just pass the input sample to the first biquad, and it’s output to the second, that output to the third, and so on.

// create and initialize

Biquad *bq1 = new Biquad();

Biquad *bq2 = new Biquad();

bq1->setBiquad(bq_type_peak, 1000.0, 0.707, 0.0);

bq2->setBiquad(bq_type_peak, 2000.0, 0.707, 0.0);

…

// process

output = bq1->process(input);

output = bq2->process(ouput);

Thanks Nigel.

tried it and i am getting no where with it.

I am working at 48000 with single samples ( floats ) that are being placed in a buffer then sent to the sound card ( the highest sample is (-6.0 0.0 +6.0) i think the such low samples are the problem).

and its completely stumping me why it wont apply.

for some reason the 3 band low mid and high ( from an other site ) works well but i cant apply anything else that ive tried , its a old sid ( c64 music player ) im working on btw.

any tips would be great on how to get it working.

Hi Dave—Probably not enough info for me to make helpful comments, but I’m not sure what you mean about the sample values (-6.0…). If your sound driver accepts floats, normally the full amplitude range is -1.0…1.0. If it accepts integers, then of course they would be large numbers, depending on whether 16 bit or 24.

But you should be able to look at the sample values you’re getting out of the 3-band that works, and the code you’re trying to get to work. You can look at raw numbers, but it’s probably easier to write them to a file and look in a audio editor, such as Audacity. You could also put in a signal (a sine wave perhaps) and see what it looks like coming out of just one biquad, initially set with no gain change, and go from there.

Hi Nigel , that is a good idea and thankyou for your help , It’s appreciated.

This source code was incredibly usefull for me to implement four LPF filter.

I was able to implement the algorithm with a fixed point cpu and it works fine on a really low cpu : dingo a320 400Mhz MIPS without floating point and 32Mo Ram.

The source code is available for anyone on my github account : https://github.com/yoyz/audio

I’d like to really thank you because i was seeking this for a long time.

Johann

Hello, Nigel, I would like to thank for sharing this. Amazing what I found.

I not using this for audio signal. My data is the capture from ADC of my controller.

each value could be from 0 to 15000 max. The capture rate is 200 Hz and I want to remove 3 Hz from the captured signal. I added the code it compiles and runs but the

‘out’ is unfortunately a zero all the time.

Do I need to go to Biquad calculator first and find out coefficients and then put those into

class where init (Biquad::Biquad()) ?

My coefficients are:

a0 = 0.002080565890575604

a1 = 0.004161131781151208

a2 = 0.002080565890575604

b1 = -1.8668911626483358

b2 = 0.8752134262106381

some my snippet:

Biquad *lpFilter = new Biquad(); // create a Biquad, lpFilter;

int Fc=3; // 3Hz to remove

int sampleRate=200;

// void Biquad::setBiquad(int type, double Fc, double Q, double peakGainDB)

lpFilter->setBiquad(bq_type_lowpass, Fc / sampleRate, 0.707, 1);

size_t flt=12;

char captureFilename[40];

sprintf(captureFilename,”%s”,”filtered.txt”);

remove(captureFilename);

debugfile1=fopen(captureFilename,”w+”);

char ta_float[16];

for(;;)

{

sprintf(ta_float,”%d”,truck_array[1][flt]);

fprintf(debugfile1, “%5d,%5.2f,\r\n”,

truck_array[1][flt],

lpFilter->process(atof(ta_float)));

flt++;

if(flt>forrecno1) break;

}

fflush(debugfile1);

fclose(debugfile1);

Hi Alex—You’ve got the general idea; the call to lpFilter->setBiquad calculates the coefficients, so you don’t need to use the biquad calculator. I think that you problem might simply be that your declare Fc and SampleRate to be type int, so when you pass Fc/sampleRate to set the filter, 3/200 evaluates to 0. You need to cast one to float.

Nigel

Thanks Nigel,

I successfully implemented a 10-band graphic equalizer (1 octave spanning 32 to 16000Hz) using the hints and snippets on your blog (using 10 peak filters).. The question I have is about choosing the right Q value. At the moment I set Q to 0.707 for every band but here raising the gain of one band raises the gain of the neighbor bands also significantly.. So is there a way to choose Q depending on the number of bands (2/3 octave, 1/3 octave etc.)?

Thanks a lot! Mario

Interaction and Q in graphic equalizers is a fairly deep subject, Mario. But there is a lot of good material available. Do a web search on something like “digital graphic equalizer choice of Q”. I’ve seen some good papers from Rane on the subject, for instance.

Thank you for this! I’m trying to find my way from implementing DSP in Reaktor to doing it in C++ and I think this will help me. I have a question though. How do you implement unit delays in C++? That’s the only thing stopping me from implementing the filters I created in Reaktor in C++. In Reaktor, there’s a module exactly for unit delays, called “z-1” (or something like that).

Hi Jordan. Unit delays are simply a variable that holds the last state. See the Biquad::process function above in Biquad.h. See how the code uses z1 in the output calculation, then calculates z1 for the next time the function is called. So z1 is a unit delay. Also note that z2 is calculated, but not used until next time for the z1 calculation, and the resulting value is used inn the output calculation the time after that—so the z2 calculation is delayed two samples. Longer delays are usually done with a circular buffer (it’s faster to update pointers that to move data from variable to variable).

Maybe this is more obvious; assuming you have a class variable z1 (double), this returns the previous input sample:

Respected Nigel,

I have a question how to pass real audio data to your filter?

For example, if I have a raw audio data from a WAV file, what is a basic concept of data passing to your filter?

You’d first read the wav file into a floating point array. You can either roll you own code for that (I do), or use something like libsndfile (I haven’t used it, but it is popular). Then it’s a matter of calling the filter’s process routine for each value in the array, and writing that to an array (it could be the same one), then writing that out as a file.

I have tried to implement your code and it works very well. With your example, when I set a gain of 6dB Bump at 200Hz, the sound is crackling and appears some artifacts. Do you know why? ¿What could I do to fix it?

lpFilterRight200->setBiquad(bq_type_peak, 200.0 / 44100.0, 1.0, 6.0);

Clipping your audio output, possibly?

Hallo, Nigel.

An impressive code and i managed to implement it in my first vst plugin. Thank you for sharing! I don’t have much of programming experience, except for some small stuff in c++ and also ActionScript way back, so i am still quite the newbie.

The only problem is that, i can hear a lot of artifacts being added to the signal, already at +1db. I read some of the other comments, where the same problem is mentioned. Obviously nothing is clipping and i have the Filter creation seperated from the processing part. I am modulating the variables there also, which, as i read in your comment, is not really suggested.

Filter_5->setPeakGain(mFilter_5 * 100);

*output = Filter_5->process(*output);

and i defined the Variable like this:

GetParam(mFilter_5)->InitDouble(“Hi Mids”, 0, -12, 12, 0.1);

I had to multiply it with 100, for it to work at all, which is what i also don’t really understand, because the value should be given in in dB, not in dB * 100, right?

Or should i try getting my feet wet with oversampling?

Your message (and the follow-up attempt) were kicked to the spam folder—I’m not sure why offhand.

OK, it’s clear you’re using wdl with IPlug (wdl-ol or other). First, I assure you that the filter code works fine, and does not glitch or distort, so this is really an issue of debugging your first plugin. Your best bet is to get in an interactive discussion on a board like kvraudio, or the wdl discussion in the cockos forums. You’ll need to show a lot more of your code.

In a nutshell, be sure that you understand that the filter is instantiated in your initialization, that control changes (the setPeakGain call, for instance) are handled via your OnParamChange handler, and the process call is in your ProcessDoubleReplacing handler.

If you really want to modulate the filter on each sample: First, you probably should look at a different filter (one that has frequency, Q, and gain as separate controls, basically). Second, for any type of filter, you’d probably want to minimize the math (see how the calcBiquad call has some relatively expensive math functions, especially tan), if you intend to do a significant amount of DSP.

Thank you for the tips. What would be the simplest form of Filter calculations, as in which types? And what is the advantage of biquads over other forms, then, if they are so CPU heavy? Is it just that they are way more precise, or is it the flexibility of setting up a specific filter with diverse roll offs and resonances?

They are not CPU heavy. The problem is that you need to change all of the coefficients to change one parameter (Q, gain, or frequency). Biquads aren’t bad for general EQ work, but if you need to modulate the parameters quickly, other choices such as the state variable are much better suited. Don’t use a biquad for a synth filter. See my recent articles.

P.S. as i noticed just now – the amount of artefacts seems to decrease drastically, as i increase the buffer size. 64 was the lowest i was using and the oscilloscope i was using showed artefacts all over the sinus wave, where as at 2048 samples it’s just 1 spikes appearing, dissapearing and reappearing again.

So does this mean, that at the lower Buffer size the filter cannot calculate as fast, as at the higher buffer sizes? And how would i fight against this problem, without having to introduce a lot of latency?

Your filter is instantiated once—not in the processing function? And filter settings are remaining constant (or not changing often)? It seems like the filter is functioning fine during processing, but the states are being disrupted in between.

Oskar, I think you figured it out already. But I want to help others searching for the same Issue: Crackling noise that decreases with higher buffer sizes. My (dumb) mistake was I was processing a stereo signal with only one biquad filter that ran over left and right buffer sequentially. The more the stereo signal differs (at the beginning/end of the buffer from left to right) there is a “step” in z1/z2 that makes that crackling noise and obviously decreases with longer buffer sizes.

Simple solution: use a separate biquad filter object for left and right buffer.

Took me a while to figure that out but is just so logic, I hit my head on the table 😉 . Thank you Nigel for the great work.

Thanks for that great example code!

I needed a subsonic filter for an application I’m developing, and with the code you posted here I was able to implement it without having to spend a lot of time reading up on IIR filter implementations & design.

I have a question though: Could you recommend a forum or mailing list for discussing stuff related to digital audio processing? We would like to integrate a compressor into our application, which is something that I definitely will require help with. And I have no idea where the right place to ask such a question would be.

The kvraudio DSP forum is a good place to ask questions, especially for filters—do a search there to find some excellent discussions on synth filters. For more purely technical DSP discussions, perhaps less on specific implementations, try the music-dap mailing list (search “music dsp columbia”).

Thank you so much! I guess I will try my luck at the kvraudio DSP forum.

Thanks for your code. It would be very helpful.

I used your code is applied for my ECG filtering. Should it work?

I tried to apply the LPF, HPF and Notch filter as following:

bufout_float[idx] = lpFilter->process(hpFilter->process(notchFilter->process(bufin_float[idx])));

But the output seems not correct. May I know where wrong is? Thanks for your time.

One more thing, Can you help to explain the Q parameter. Does it affect to the filtered result? If i have the sampling rate is 125, then what should be Q value? Thanks

Yes, that chain should work fine. I couldn’t tell you what might be wrong.

For an understanding of Q setting, play with the biquad calculator. Q is independent of sample rate.

Dear Nigel,

Thank you for the Biquad code. I only need one more piece of information before I can put this code to use. Could you provide the math for bq_type_allpass?

Regards,

Ben H

How do I implement theDirect Form II Canonic biquad filter structure, if I have one zero and 2 poles with a common gain. My transfer func in freq domain is follows:

s + 400

—————

s^2 + 12 s + 20 .

After bilinear transform I have my z-domain tf as follows:

0.005 + 0.004706 z^-1

—————————-

1 – 1.941 z^-1 + 0.9418 z^-2

How shall I implement the same. Does this be achieved with one section?

Thanks in advance.

OK, from the tf, a0 = 0.005, a1 = 0.004706 (and a2 = 0.0); b0 = 1.0 (which we won’t use in a normalized structure, because we omit the multiplier), b1 = -1.941, b2 = 0.9418. In my biquad source code, there is no method to set coefficients (a0, a1…) directly, but you could add one. The Biquad::process inline function in the .h file is a transposed direct form II, so once you’ve set the coefficient, use that for the sample processing. The transformed version of DFII is better than the non-transformed version—better floating point accuracy. (In the filter frequency response grapher, put “0.005,0.004706” in the left field and “1,-1.941,0.9418” in the right, click the graph to see the frequency response.

Is there any way to calculate the time delay introduced by this filter?

There is no time delay, per se; if you look at the block diagrams, the input connects directly to the output without an intervening delay element. What you have is group delay, which is frequency-dependent (it’s the derivative of phase with respect to frequency). If you put lowpass biquad coefficients in the filter frequency response grapher and look at the log chart, you’ll get a good idea of how the phase shift is zero for much of the passband, until getting near the cutoff and increasing into the stop band.

Thank you for explaining biquad filters. It really helps me out a lot. I only have one question regarding the follow situation: I am building a 32-bands graphic equalizer and I would like to calculate the samples with 32 cascading biquad filters. Now what I’ve got out of it is 1.INF#000000 after couple of filters. Am I doing something wrong or is there anything I can do?

Thanks in advance.

It seems you filter is “blowing up”, which happens when a pole is outside the unit circle. Using my code without modification? I assume you’re using cascaded peaking filters.

Yes, that is indeed the case

First, thank you for your blog, it has been a great help for me to understand DSP.

I have converted your C++ codet to java, and am having some problems with the results. I have set up a bandpass filter with a sample rate of 48kHz and a center frequency of 20kHz. I want a 500Hz -3db bandwidth so I have set Q to 40 (20kHz/500Hz). The result I get has a much tighter bandwidth (about 95Hz). I have tried this on other tools (Purepath and http://engineerjs.com) and it seems to work as I expected. If I plug a 95Hz bandwidth into these other tools I get roughly the same coefficients as your Biquad C++ source code gives for Q = 40 ).

I am not experienced enough to understand how you derived you formulas, so I don not know if there is an error in them or if I am doing something wrong. Any help would be appreciated.

I followed Robert Bristow-Johnson’s lead on the definition of Q. Search for “rbj cookbook”. I found a couple of online calculators that produced the same results as mine. Unfortunately, I suspect the first one might use my code, but I found another that uses rbj’s, and gives the same results: BiQuadDesigner. So, I have to say it seems to be working as designed!

Thank you for the quick reply.

I found Robert Bristow-Johnson’s cook book and manually calculated the coefficients. After normalizing the the coefficients, they match what your code produced.

I also gather from the cook book that “the EE kind of definition” is used for Q in a band pass filter… (Q = fc/delta_f I assume). Am I confusing Q with something else?

Thank you for your help. It seems the problem was with my definition of Q / frequency warping. I got an answer at https://dsp.stackexchange.com/questions/47522/help-with-audio-eq-cookbook-bpf-filters-and-q

Ah—I noticed your spec of 20 kHz, but neglected to mention the frequency warping. Thanks for following up with the link.

How to design 2 biquads filters

Question?

Great stuff – out of about 10 different BiQuad source codes out there – this is the only one I’ve found that works well – thank you!

Here’s the code to get the Frequency response curve from the Coeff’s if you want to visually graph it:

double ZTransMag(double freq)

{

double w0 = 2 * M_PI * freq;

double cosW0 = cos(w0)*mag;

double cos2W0 = cos(2.0 * w0)*mag;

double num = (a0*a0) + (a1*a1) + (a2*a2) + (2.0 * cosW0 * ((a0 * a1) + (a1 * a2))) + (2.0 * cos2W0 * a0 * a2);

double den = 1.0 + (b1*b1) + (b2*b2) + (2.0 * cosW0 * (b1 + (b1 * b2))) + (2.0 * cos2W0 * b2);

return (abs(num / den));

}

What are the coefficients you used for this ?

Can you add an allpass mode to the filter?

I’ve found these coefficients:

b0 = 1 – K/Q + K^2

b1 = 2*(K^2 – 1)

b2 = 1

a0 = 1 + K/Q + K^2

a1 = 2*(K^2 – 1)

a2 = 1 – K/Q + K^2

But they don’t match your variable naming convention.

I am looking for an All Pass filter too.

I too have found other implementations that use sin/cos and not sure how to translate it to match the same processing here.

If you want second order, I just derived it…the s-domain transfer function is (s^2 – s/Q + 1) / (s^2 + s/Q + 1). If you add another case bq_type_allpass (don’t forget to add to the enum list too):

Thanks for sharing this, Nigel!

One small thing – in the header file, the definition of “void setBiquad” misses the “DB” at the end of “peakGain”.

Well, setPeakGain is in dB, I’m lazy so I don’t want to type setPeakGainDB when I’m always going to use dB…I name the variable peakGainDB as a reminder, but it’s not something you’d type…

That’s fine – I wasn’t referring to “setPeakGainDB()”, but the function “setBiquad()”. If you look closely, the last parameter of this function has a different name in the header, vs. the source file.

Ah. Yes, it should be in the header file, as it’s informative, thanks. (fixed)

How does one go about making this sound creamier?

I got it to work in a plugin I’m developing, and it behaves as it should, but it sounds significantly lower quality than other EQ plugins. Would oversampling change much?

“Creamier” often refers to the filter’s behavior when overdriven a bit, adding thickness to the sound. Oversampling is one way to get better overdriven behavior. But the filter topology can be very important if you’re going that route, especially when adding adding non-linearities in feedback. Check the kvraudio forum and other places on the web for discussion of “zero delay feedback” filters. Though these are primarily focussed on synth filters, which are more intrusive on the sound by nature than EQ filters.

Dear Nigel,

Is that scenario possible with the implementation of the Biquad lowpass filter with 6dB/octave? Thanks!

Hi Alex. So, obviously 6 dB is first order, it can be implemented as a biquad but some of the coefficients will be zero (three coefficients for a one-pole lowpass, or three for a one-pole, one-zero lowpass). There is separate code on the site for a one-pole lowpass. I don’t have source code written up for a 1p1z lowpass, but the new Biquad calculator v3 widget will calculate coefficients that you can use in the biquad C++ filter process (or, you can make a first-order process that deletes the unneeded parts). If you need the equations, you can use the web inspector of your browser to look at the widget’s source code, and get the part you need—trivial to convert to C++:

Thank you so much for sharing these Nigel!

I was wondering if you had any ideas why high pass filter isn’t giving me the results I normally expect. I am trying to filter out frequencies below 3000 Hz and the effect is barely audible or visible in a spectrogram.

By contrast, Low pass filter works perfectly and almost completely removes high frequencies, as expected.

I have tried stacking three highpass filters in a row, but the overall signal level is reduced with hardly perceivable change in low frequencies.

Is there a way to strengthen how this Biquad code cuts the lows?

Thank you very much.

Hi Blake. I don’t see anything wrong in the code on this page, and it matches the code for the biquad calculator, which plots a correct highpass. I Suppose I’d have to know more about exactly what you’re doing. Maybe a good first step is to see that your code and the biquad calculator produce the same coefficients. The biquad calculator plots directly from the coefficients, so if it produced incorrect coefficients, the plot would also be wrong.

Hi Nigel,

following this topic https://forum.juce.com/t/iir-peak-juce-vs-nigel-implementation/54517/1, would it possible to implement a variations on the code above to have q “engineering” instead of “proportional”?

Not so able for the match part hehe.

Simple changing some variable? Or the whole coefficients need to be refactored?

Thanks

For the approach of designing from analog prototypes (s-domain), you’d need to start with the prototype that implements Q the way you want, then perform the bilinear z transform. I’ll have to leave that up to you.

Hi Nigel.

First, thanks a lot for good explanation you provide in your articles.

I have a question re biquad coefficients calculation. There’s a software that refers to this article https://webaudio.github.io/Audio-EQ-Cookbook/Audio-EQ-Cookbook.txt but implements a bit different formulae for band-reject filter, e.g.

omega = 2.0 * M_PI * fc / fs;

sn = sin(omega);

cs = cos(omega);

alpha = sn * M_SQRT1_2 * exp(-M_LN10 / 20 * r);

a0r = 1.0 / (1.0 + alpha);

d.b0 = a0r;

d.b1 = a0r * (-2.0 * cs);

d.b2 = a0r;

d.a1 = a0r * (2.0 * cs);

d.a2 = a0r * (alpha – 1.0);

They use some ‘r’, that’s called ‘resonance’ and is neither Q nor bandwidth.

When it’s zero, BW and center frequency attenuation are maximal, BW is about fc*sqrt(2), that means that Q is about 1/sqrt(2). BW and attenuation decrease non-linearly with this ‘r’ increasing.

Maybe you can give me a clue about what this ‘resonance’ is, what its physical meaning is, and how it’s mathematically related to Q and BW? Can we assume that

alpha = sn * M_SQRT1_2 * exp(-M_LN10 / 20 * r) = sn/(2*Q) ?

Software authors are silent, unfortunately.

Thanks.

I’m working from my phone in a far away place, but note there needn’t be “meaning” if this is aimed at musical use, other than it gives the kind of peaking the authors want.

Good point.

The most difficult issue is that there’s no convention for coefficients calculation, as you mentioned here https://www.earlevel.com/main/2017/01/04/about-coefficient-conventions/. Will read further, maybe I’ll found a way to calculate Q knowing this ‘r’.

I don’t understand the absence of b0? How can I put it into matlab code that calculates y(n) = B0 * x(n) + B1 * x(n_1) + B2 * x(n_2) – A1 * y(n_1) – A2 * y(n_2);?

Be aware that matlab uses the convention of B on the input side of the type 1 biquad, feeding forward, and A coefficients feeding back the output. So use my B1 for A1 in your formula, etc. Sorry for the confusion, but unfortunately there is no universal convention. I chose to use the order than most of my text books used at the time.

Hi!

The high pass option does not work for me, I set 20000hz for a 44100hz file and get very bass sounding sound. seems like it does some inverse stuff. I really made sure the right type was chosen. Quite useless for me. I use cbuilder. Very strange.

Thanks anyway.

What you describe seems consistent with have the

aandbcoefficients swapped. Are you using my code to both calculate and execute the filter? I can double check when I have more time, but this same code is used in commerical plugins, including multiple highpass filters, and gives the expected response. For a start, maybe you could list out the coefficients it generates for the setting.more :

result :

— hipass : Fc = 0.4989, Q = 0.7070, norm = 0.0000, V = 1.0000, K = 280.7481, a0 = 0.0000, a1 = -0.0000, a2 = 0.0000, b1 = 1.9899, b2 = 0.9900

code i use:

double norm;

double V = pow(10.0, fabs(peakGain) / 20.0);

double K = tan(M_PI * Fc);

case bq_type_highpass:

// error here ?? wrong pass ?? seems to do the opposite

norm = 1.0 / (1.0 + K / Q + K * K);

a0 = 1.0 * norm;

a1 = -2.0 * a0;

a2 = a0;

b1 = 2.0 * (K * K – 1.0) * norm;

b2 = (1.0 – K / Q + K * K) * norm;

printf(”

— hipass : Fc = %2.4f, Q = %2.4f, norm = %2.4f, V = %2.4f, K = %2.4f, a0 = %2.4f, a1 = %2.4f, a2 = %2.4f, b1 = %2.4f, b2 = %2.4f

“, Fc, Q, norm ,V, K, a0 ,a1 ,a2, b1 , b2);

break;

…

The norm value looks suspicious.

The result with 22000.0 is quite good, just the opposite of what you would expect.

If I go down to 20000.0 the audio sounds almost like unprocessed.

Could this be a precision problem ? Should I also run the two channels separately? One class for each channel?

…

Here’s what I don’t understand: Why are you using a highpass filter set close to half the sample rate? That would only pass high frequencies. set at 20k, the response would be down -12 dB by 10k, -24 dB at 5k…

Hi again !

I tried to split the two channels in my file (stereo) into samples thrown into two different classes of biquad. That seemed to do the trick. I have no clue why.

Since I read continuous buffers from a file I only though the mixing of left / right samples would effect the quality slightly. Very strange .

Maybe you have the explanation.

Thanks!

Ah—yes, the channels must each have their own instance of a filter. Realize that IIR filters have feedback, so if you mix the samples, the output for the right will be affected by the memory of the left, and so on. And it’s not the same as summing the two channels, since they alternate into the filter. Glad you got that figured out.

Well, if I set the freq to over 22050, the program crashes .

Some K = tan(M_PI * Fc) problem, precision problem ?

And I am only interested in what you can actually hear anyway.

And why would the response be down in Db?

Is there a way to extract what actually is there?

And how do I input arguments for a bq_type_bandpass solution?

Thanks again.

You can’t put the frequency over 22050 because that’s over the Nyquist frequency (half the sample rate), it’s undefined.

I think you’re missing the fact that a highpass filter blocks what’s below. So if you set it to 20000, you’re blocking everything in the audio band. The rolloff is 12 dB/oct, since it’s a second other filter, so that’s why I said the response is down 12 dB at 10k.

Put your settings in Biquad calculator v3 to see the response.

Hi!

Well, I did not miss that fact, I can actually hear the sound when setting 20000hz, but low as you explained the reason for.

Thanks for your help and tips !!!