Butterworth Lowpass Filter Coefficients in C++

When I searched the web it wasn't easy to find a simple function to calculate Butterworth lowpass filter coefficients, so I looked up the theory and did things on my own.

If you just want coefficients for a single filter specification, you can use the code generated by this nice online code generator. It will produce C-code that you can use for a specified samplerate and cutoff-frequency. It will also do Chebyshev and Bessel filters as well as lowpass, highpass, bandpass and bandstop filters for you.

However, if your program needs adjustable cutoff or samplerate, you can not use this script.

The first part of this post will take a glance at the required analytical steps and the second part will simply give you the code.

One of the standard techniques is to use the analytical frequency response of the filter and then use bilinear transformation to obtain IIR filter coefficients.

The required steps require some math for which I found this lecture very helpful:

1) Take analytical transfer function: H(s)
2) Do bilinear transform: H( (z-1)/(z+1) )
3) "Warp" cutoff frequency to find cutoff in the "bilinear domain"
4) Express function as
5) Use these coefficients for time-domain filtering with the linear difference equation:

Notice that in order to simplify the derivation of the filter, you can safely set the sample time T to 2.


Doing this for a 2 pole Butterworth filter gives the following code in C++:

void getLPCoefficientsButterworth2Pole(const int samplerate, const double cutoff, double* const ax, double* const by)
{
    double PI      = 3.1415926535897932385;
    double sqrt2 = 1.4142135623730950488;

    double QcRaw  = (2 * PI * cutoff) / samplerate; // Find cutoff frequency in [0..PI]
    double QcWarp = tan(QcRaw); // Warp cutoff frequency

    double gain = 1 / (1+sqrt2/QcWarp + 2/(QcWarp*QcWarp));
    by[2] = (1 - sqrt2/QcWarp + 2/(QcWarp*QcWarp)) * gain;
    by[1] = (2 - 2 * 2/(QcWarp*QcWarp)) * gain;
    by[0] = 1;
    ax[0] = 1 * gain;
    ax[1] = 2 * gain;
    ax[2] = 1 * gain;
}

which then can be used in a filter like this:

double xv[3];
double yv[3];

void filter(double* samples, int count)
{
   double ax[3];
   double by[3];

   getLPCoefficientsButterworth2Pole(44100, 5000, ax, by);

   for (int i=0;i<count;i++)
   {
       xv[2] = xv[1]; xv[1] = xv[0];
       xv[0] = samples[i];
       yv[2] = yv[1]; yv[1] = yv[0];

       yv[0] =   (ax[0] * xv[0] + ax[1] * xv[1] + ax[2] * xv[2]
                    - by[1] * yv[0]
                    - by[2] * yv[1]);

       samples[i] = yv[0];
   }
}

And that is it. To reduce summation error you may want to use Kahan summation!