Implement 2nd order low pass filter in C++, how to compute coefficients?

Question

I am struggling to find the proper algorithms for generating the coefficients for low pass filters. I wrote the following taking the butterworthLowPass code from another SO question:

class Filter {
    float fs;
    float a1, a2, b0, b1, b2;
    float v1, v2;

   public:
    Filter(float fs) : fs(fs) {}

    float compute(float x)
    {
        float v0 = x - a1 * v1 - a2 * v2;
        float y = b0 * v0 + b1 * v1 + b2 * v2;
        v2 = v1;
        v1 = v0;
        return y;
    }

    void butterworthLowPass(float fc)
    {
        const float ita = 1.0 / tan(M_PI * fc / fs);
        const float q = sqrt(2.0);
        b0 = 1.0 / (1.0 + q * ita + ita * ita);
        b1 = 2 * b0;
        b2 = b0;
        a1 = 2.0 * (ita * ita - 1.0) * b0;
        a2 = -(1.0 - q * ita + ita * ita) * b0;
    }

    void chebyshevLowPass(float f, float ripple)
    {
        // TODO: implement
    }

    void ellipticLowPass(float f, float ripple, float stopband)
    {
        // TODO: implement
    }

    void displayCoefficients()
    {
        fprintf(stderr, "a = [ %f, %f, %f ]\n", 1.f, a1, a2);
        fprintf(stderr, "b = [ %f, %f, %f ]\n", b0, b1, b2);
    }
};

Here the main:

int main()
{
    const float fs = 10000;
    Filter biquad(fs);
    biquad.butterworthLowPass(1000);
    biquad.displayCoefficients();
}

When taking the computed coefficients and using Matlab to display the filter response, I get a -15dB gain which is not what I expect.

a = [ 1.000000, 1.142980, -0.412802 ];
b = [ 0.067455, 0.134911, 0.067455 ];
[mag, phase, wout] = bode(tf(b,a,1/(fs/2)));
subplot(2,1,1);
semilogx(wout(:,1)/(2*pi), 20*log10(squeeze(mag)), '-b'); zoom on; grid on;
axis tight
ylim([-40 0])
title('magnitude'); xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)');
subplot(2,1,2);
semilogx(wout(:,1)/(2*pi), squeeze(phase), '-r'); zoom on; grid on;
axis tight

Matlab gives me something very different:

>> fs = 10000;
>> fc = 1000;    
>> [b, a] = butter(2, 1000/10000)
b =
    0.0201    0.0402    0.0201
a =
    1.0000   -1.5610    0.6414

Further experiments

I tried to clone ruohoruotsi/Butterworth-Filter-Design.git and test with:

#include "Butterworth.h"

using namespace std;

int main() {
  float fs = 20000;
  float fc = 500;
  int order = 2;

  vector<Biquad> coeffs; // second-order sections (sos)
  Butterworth butterworth;

  bool designedCorrectly = butterworth.loPass(fs, fc, 0, order, coeffs, 1.0f);

  std::cout << "Designed correctly? " << designedCorrectly << std::endl;
  std::cout << "a = [" << 1.f << ", " << coeffs[0].a1 << ", " << coeffs[0].a2
            << "]" << std::endl;
  std::cout << "b = [" << coeffs[0].b0 << ", " << coeffs[0].b1 << ", "
            << coeffs[0].b2 << "]" << std::endl;
}

It gives me:

a = [1, 1.77863, -0.800803]
b = [1, 2, 1]

Which again doesn't really fit.

Answer 1

I dont have access to Matlab, but on python

butter(2,.2) = [0.06745527, 0.13491055, 0.06745527], [ 1., -1.1429805, 0.4128016]

and

butter(2,.1) = [0.02008337, 0.04016673, 0.02008337],[ 1. , -1.56101808, 0.64135154]

So your values from Matlab are off by a factor of "2" on FF (I once made same mistake and posted a comment in Calculate Coefficients of 2nd Order Butterworth Low Pass Filter)

The "Butterworth.h" code has normalized the "b" values by multiplying by (1.0 + q*ita + ita*ita). Your "a" values have minus signs compared to what Matlab, Scipy, and the original SO post use, I didnt check how that contributes too.

Answer 2

For Butterworth you could use this:

void butter(float fc, float fs, bool lowPass = true) {
    const float k = tanf(M_PI * fc/(fs*2));
    const float k2 = k * k;
    const float q = 1 / sqrtf(2);
    const float norm = 1 / (1 + k / q + k2);
    if (lowPass) {
        b0 = k2 * norm;
        b1 = 2 * b0;
        b2 = b0;
    } else {
        b0 = 1 * norm;
        b1 = -2 * b0;
        b2 = b0;
    }
    a1 = 2 * (k2 - 1) * norm;
    a2 = (1 - k / q + k2) * norm;
}

In Matlab you get the exact same coefficients with butter(2, fc/fs)

Answer 3

If the transfer function after a z-transform is defined as

H(z) = ( b0 + b1/z + b2/z^2 ) / ( a0 + a1/z + a2/z^2 )

then I think that the formulation for the coefficients a1 and a2 is producing values of the wrong sign. (At least, that's what I found when reworking it, but I might be wrong).

The following code gives the opposite signs for a1 and a2.

#include <iostream>
#include <iomanip>
#include <complex>
#include <vector>
#include <cmath>
using namespace std;

const double PI = 4.0 * atan( 1.0 );
const complex<double> iPI( 0.0, PI );
const double SQRT2 = sqrt( 2.0 );

// const double factor = 0.5;    // use for Matlab comparison
const double factor = 1.0;       // as original post

//----------------------------------------------------------------------

void getCoefficients( double fd, double fs, vector<double> &a, vector<double> &b )
{
   double alpha = 1.0 / tan( factor * PI * fd / fs );
   b[0] = 1.0 / ( 1.0 + SQRT2 * alpha + alpha * alpha );
   b[1] = 2 * b[0];
   b[2] = b[0];
   a[0] = 1.0;
   a[1] = 2.0 * ( 1.0 - alpha * alpha ) * b[0];
   a[2] = ( 1.0 - SQRT2 * alpha + alpha * alpha ) * b[0];
}

//----------------------------------------------------------------------

void bode( const vector<double> &a, const vector<double> &b, double f, double fs, double &mag, double &phase )
{
   complex<double> invz = ( 1.0 - iPI * factor * f / fs ) / ( 1.0 + iPI * factor * f / fs );
   complex<double> H = ( b[0] + b[1] * invz + b[2] * invz * invz ) / ( a[0] + a[1] * invz + a[2] * invz * invz );
   mag = abs( H );
   phase = atan2( H.imag(), H.real() );
}

//----------------------------------------------------------------------

int main()
{
   vector<double> a(3), b(3);
   double fs = 10000;
   double fd = 1000;

   // Get Butterworth coefficients in H(z) = ( b0 + b1/z + b2/z^2 ) / ( a0 + a1/z + a2/z^2 )
   getCoefficients( fd, fs, a, b );
   cout << "a:\t";
   for ( double e : a ) cout << e << '\t';
   cout << "\nb:\t";
   for ( double e : b ) cout << e << '\t';
   cout << "\n\n";

   double logmin = 1.0, logmax = 4.0;
   int nw = 31;
   double dlog = ( logmax - logmin ) / ( nw - 1 );
   #define FMT << scientific << setw(15) << setprecision(4) <<

   cout FMT "f(Hz)" FMT "mag" FMT "db" FMT "phase(deg)" << '\n';
   for ( int i = 0; i < nw; i++ )
   {
      double f, mag, phase;
      f = pow( 10.0, logmin + i * dlog );
      bode( a, b, f, fs, mag, phase );

      double phasedeg = phase * 180.0 / PI;
      double db = 20.0 * log( mag ) / log( 10.0 );
      cout FMT f FMT mag FMT db FMT phasedeg << '\n';
   }
}

Output:

a:  1   -1.14298    0.412802    
b:  0.0674553   0.134911    0.0674553   

          f(Hz)            mag             db     phase(deg)
     1.0000e+01     1.0000e+00    -3.7956e-08    -7.8347e-01
     1.2589e+01     1.0000e+00    -9.5341e-08    -9.8635e-01
     1.5849e+01     1.0000e+00    -2.3949e-07    -1.2418e+00
     1.9953e+01     1.0000e+00    -6.0156e-07    -1.5634e+00
     2.5119e+01     1.0000e+00    -1.5111e-06    -1.9683e+00
     3.1623e+01     1.0000e+00    -3.7956e-06    -2.4783e+00
     3.9811e+01     1.0000e+00    -9.5341e-06    -3.1205e+00
     5.0119e+01     1.0000e+00    -2.3949e-05    -3.9296e+00
     6.3096e+01     9.9999e-01    -6.0156e-05    -4.9494e+00
     7.9433e+01     9.9998e-01    -1.5110e-04    -6.2354e+00
     1.0000e+02     9.9996e-01    -3.7954e-04    -7.8588e+00
     1.2589e+02     9.9989e-01    -9.5331e-04    -9.9113e+00
     1.5849e+02     9.9972e-01    -2.3942e-03    -1.2513e+01
     1.9953e+02     9.9931e-01    -6.0114e-03    -1.5821e+01
     2.5119e+02     9.9826e-01    -1.5084e-02    -2.0052e+01
     3.1623e+02     9.9566e-01    -3.7791e-02    -2.5501e+01
     3.9811e+02     9.8920e-01    -9.4310e-02    -3.2581e+01
     5.0119e+02     9.7352e-01    -2.3312e-01    -4.1849e+01
     6.3096e+02     9.3720e-01    -5.6339e-01    -5.3957e+01
     7.9433e+02     8.6132e-01    -1.2967e+00    -6.9313e+01
     1.0000e+03     7.3050e-01    -2.7276e+00    -8.7273e+01
     1.2589e+03     5.5943e-01    -5.0451e+00    -1.0563e+02
     1.5849e+03     3.9180e-01    -8.1387e+00    -1.2189e+02
     1.9953e+03     2.5949e-01    -1.1718e+01    -1.3493e+02
     2.5119e+03     1.6715e-01    -1.5538e+01    -1.4496e+02
     3.1623e+03     1.0636e-01    -1.9464e+01    -1.5262e+02
     3.9811e+03     6.7339e-02    -2.3435e+01    -1.5850e+02
     5.0119e+03     4.2546e-02    -2.7423e+01    -1.6305e+02
     6.3096e+03     2.6859e-02    -3.1418e+01    -1.6660e+02
     7.9433e+03     1.6951e-02    -3.5416e+01    -1.6939e+02
     1.0000e+04     1.0696e-02    -3.9415e+01    -1.7159e+02

As for Matlab ... well, if you change the factor at the top of that code to the commented-out value (0.5) then you would get the Matlab values:

a:      1       -1.56102        0.641352        
b:      0.0200834       0.0401667       0.0200834