IIR filter with circular buffer

Using Arduino Programming Questions
1

Hello,
I am trying to implement a butterworth lowpass filter with a circular buffer. The code below does not really filter as intended and I am not really able to understand what is wrong. Any suggestion?

#define PERIOD_MICROSECS 1000

static uint32_t lastRead = 0;
     

int analog_pin = 0;   
int32_t analog_input0 = 0;       
double analog_input0_lp_filtered = 0;



// The coefficients are calculated offline. I set the sample frequency to 1000 Hz, the cutoff to 50 Hz.
// The filter has order = 4.
double c[] = {0.000416599204407, 0.001666396817626, 0.002499595226440, 0.001666396817626, 0.000416599204407};
double d[] = {1.000000000000, -3.180638548875, 3.861194348994, -2.112155355111, 0.438265142262};

int m = 10;
double x[10];
double y[10];

int n = 0;



double lp_butterworth(double input_sample, double *x, double *y, int m, int n)
{
    x[n] = input_sample;
    y[n] = d[0] * x[n] + d[1] * x[(n-1+m)%m] + d[2] * x[(n-2+m)%m] + d[3] * x[(n-3+m)%m] + d[4] * x[(n-4+m)%m]
                     - c[1] * y[(n-1+m)%m] - c[2] * y[(n-2+m)%m] - c[3] * y[(n-3+m)%m] - c[4] * y[(n-4+m)%m];
    
    return y[n];
}            





// ************************* //

void setup() {

  Serial.begin(115200);
}




void loop() {
  
  if (micros() - lastRead >= PERIOD_MICROSECS) {
        lastRead += PERIOD_MICROSECS;

        analog_input0 = analogRead(analog_pin);  
        analog_input0_lp_filtered = lp_butterworth(analog_input0, x, y, m, n);
        n = (n + 1) % m; 
        
        //Check the original and filtered signals with the serial plotter
        Serial.print(analog_input0); 
        Serial.print(" ");
        Serial.println(analog_input0_lp_filtered); 
  } 
}
2

Interesting, where did you find this formula?

3

The code below does not really filter as intended

The code does something. You forgot to tell us what it actually does.
You intend for the code to do something, but you forgot to tell us what you intended it to do.

Repeatedly calculating (n-x+m)%m is rather expensive, when n and m are constant for any given call to lp_butterworth.

4

Hi, I found that formula here:

I simply want to implement a lowpass butterworth filter.
There must be something wrong in the code. I am getting the coefficients from Signal Processing .

5

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!

It seems that this code was already problematic

6

You don't indicate what processor you're using, but I suspect that code will be problematic on any AVR processor where "double" is only 32-bit (same as float). The coefficients range over 5 decades, and floats have a range of only about 6 decades. So, you may well be processing mostly noise, and rounding errors.

I suggest you try the same code on a PC or ARM-based Arduino, using "real" 64-bit doubles, and see if it behaves any better.

Regards,
Ray L.

7

Thanks you are all right. Actually the error was that I confused the c and d coefficients. Yes, ARM 32-bit is the way to go rather than AVR

8

FYI, there is the ARM DSP library which provides everything you need for your digital filtering (I think one of the convolution functions should do the trick):

https://www.keil.com/pack/doc/CMSIS/DSP/html/group__groupFilters.html

9

Thanks a lot