Implementation of IIR Lowpass Filter in C++

Question

I have the filter coefficients for a Butterworth lowpass filter, coming from the Matlab function [b, a] = butter(3, 0.4, 'low') and I am implemented the same computation Matlab is using according to their documentation for filter(b, a, X). The results for filtering a constant signal of 5.0, for example, are the same, but only for the first 10 values!?

I suppose that my circular buffer is wrong, but I cannot find any problems. The values are written correctly in x using the filter method, the arrays are initialized with zeros, the circular buffer pointer n has the right values... Do you have any ideas?

// Interface
class LowpassFilter {
private:
    double x[10]; // input vector
    double y[10]; // output vector
    int n;    // pointer to the current array index

public:
    LowpassFilter(); 
    double filter(double sample);
};


// Implementation

// filter coefficients a and b
const double a[] = {1.0, -0.577240524806303, 0.421787048689562, -0.056297236491843};
const double b[] = {0.098531160923927, 0.295593482771781, 0.295593482771781, 0.098531160923927};
static int c = 0;

LowpassFilter::LowpassFilter() : x{0.0}, y{0.0}, n(0) { } // Constructor

double LowpassFilter::filter(double sample)
{
    x[n] = sample;
    y[n] = b[0] * x[n] + b[1] * x[(n-1)%10] + b[2] * x[(n-2)%10] + b[3] * x[(n-3)%10]
                       - a[1] * y[(n-1)%10] - a[2] * y[(n-2)%10] - a[3] * y[(n-3)%10];

    std::cout << c++ << ": " << y[n] << std::endl; // for checking the result with the Matlab results

    double result = y[n];
    n = (n + 1) % 10; // new pointer index 
    return result;
}

What happens with the modulus when say n == 0 and (n-2)%10? [Check this](http://stackoverflow.com/questions/7594508/modulo-operator-with-negative-values). — emsr, Mar 20 '14 at 16:15
@MichaelDorner: The modulo *does* have negative values, and that's the cause of the problem. Replace e.g. `(n-1)%10` with `(n+10-1)%10` to ensure you always get the expected positive value. — Mike Seymour, Mar 20 '14 at 16:42

score 3 · Accepted Answer · edited May 23 '17 at 10:30

Thanks to Mike Seymour and emsr the problem were the negative indexes in the computation of y[n]. To solve the problem only one line has to be adopted:

y[n] = b[0] * x[n] + b[1] * x[(n-1+m)%m] + b[2] * x[(n-2+m)%m] + b[3] * x[(n-3+m)%m]
                   - a[1] * y[(n-1+m)%m] - a[2] * y[(n-2+m)%m] - a[3] * y[(n-3+m)%m];

to make sure that the index is positive. Now it works fine. Thanks alot!

Implementation of IIR Lowpass Filter in C++

1 Answers1