Biquad Filter C++ Source Code

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]));
}
https://www.earlevel.com/main/2012/11/26/biquad-c-source-code/