0
$\begingroup$

I am trying to implement a butterworth lowpass filter with a circular buffer. The code below does not 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); 
  } 
}

Below you can see the result of the code in the serial plotter using as analog input an accelerometer at steady state. You can notice a huge offset and not really a filtering effect.

Serial plotter result

$\endgroup$
  • 1
    $\begingroup$ Can you describe in what way it doesn't work like it should? Also, I'm always very skeptical when people say "I used a digital Butterworth", because that is often just because that's the first thing they found in their textbook, but almost never what they actually need. So, could you explain why you're using a Butterworth IIR here? What's the purpose/greater applicational scope of all this? $\endgroup$ – Marcus Müller May 22 '19 at 12:20
  • $\begingroup$ Simply the result is not really a signal result of filtering. There must be something wrong in the code. I am getting the coefficients from exstrom.com/journal/sigproc . I have update the post with a screenshot of the serial plotter. I am willing to filter various analog sensors (including accelerometers and infrared sensors) and I want to be able to control the cutoff frequency. The butterworth is known to be have maximally flat response, which is what I want. $\endgroup$ – L_T May 22 '19 at 12:56
  • $\begingroup$ you'll need to be more precise than "Not a result of filtering". Please test with a sensible test signal (for example $1,0,0,\ldots, 0$), and see how the result differs from your expectation! $\endgroup$ – Marcus Müller May 22 '19 at 13:09
  • $\begingroup$ There is no error in my testing methodology: the signal I am using is more than reasonable, and I know what I should expect from a butterworth filtering. Do you have any suggestion about the actual code? Simply the result is that there is no filtering. $\endgroup$ – L_T May 22 '19 at 13:12
1
$\begingroup$

First of all :

have you initialized all the previous data to 0 ? The x and y buffer should be initialized to 0. Otherwise, you could get some unexplained transients at the beginning. Maybe this is the offset that you see at the beginning. Maybe the compiler will set the initial data in the data at 0, but you should not rely on the compiler, it should be explicitely reset to 0.

Secondly : I think you have mixed up the c and d coefficients. I think the d coefficients should be used with the output (y) and the c coefficients for the input (x). I say that because d[0] = 1 which makes me think that this the a0 coefficient. Usually the a0 coefficient is 1 in an IIR. It can be different than 1, but usually we try to normalize the coefficients by setting a0 to 1. https://en.wikipedia.org/wiki/Infinite_impulse_response

Third : Have you tried the coefficients in Octave/Matlab/Python ? Use freqz or Bode in Octave/Matlab. Do you get the expected frequency response?

Finally : This is not related to the matter at hand, but I notice that you use the modulo function. The modulo function is slow unless you're using powers of 2. You should consider having a circular buffer with a size of 4 instead of 10. You don't need more since you only go 4 samples in the past.

$\endgroup$
  • $\begingroup$ Thanks, the issue was that I inverted c and d coefficients... $\endgroup$ – L_T May 22 '19 at 14:58
  • $\begingroup$ @L_T without trying to be mean (honestly!), that's something that you would have started to suspect when you tested with a $1,0,\ldots$ signal, because that isolates problems pretty well :) anyway, I'm happy that the issue was resolved! $\endgroup$ – Marcus Müller May 22 '19 at 16:22

Not the answer you're looking for? Browse other questions tagged or ask your own question.