Fixed-point Exponential Moving Average filter with given cutoff frequency

Question

I want to implement in fixed point arithmetic an exponential moving average filter, with a specific cutoff frequency for a given sampling rate. The formula for the filter is $$ y_n = \alpha x_n + (1 - \alpha) y_{n-1} $$ Looking at Exponential moving average cut-off frequency one can derive the value for alpha for a given sampling rate and a wanted cutoff. For instance, for a sampling rate of 1000 Hz and a cutoff of 50 Hz, the value of alpha is 0.267730531659313:

Fs = 1000
f3db = 50 %This is the cutoff frequency that I want 
format long;
omega3db = f3db * pi/(Fs/2)
alpha = cos(omega3db) - 1 + sqrt(cos(omega3db).^2 - 4*cos(omega3db) + 3)  =  0.267730531659313

Now I tried to implement the exponential moving average filter in fixed point arithmetic as described here: https://www.embeddedrelated.com/showarticle/779.php but I am getting an offset that I can't explain.

Here my Arduino code:

#define PERIOD_MICROSECS 1000

static unsigned long lastRead = 0;

int analog_pin = 0;   // sensor connected to analog pin 0
int analog_input0 = 0;  // variable to store the read value; 
int analog_input0_filtered_fixed_point = 0;


/*
 * Formula to compute the alpha parameter of the EMA given a desired cutoff frequency. 
 * Taken from: https://dsp.stackexchange.com/questions/40462/exponential-moving-average-cut-off-frequency
 * 
 * Fs = 1000
 * f3db = 50 %This is the cutoff frequency that I want 
 * format long;
 * omega3db = f3db * pi/(Fs/2)
 * alpha = cos(omega3db) - 1 + sqrt(cos(omega3db).^2 - 4*cos(omega3db) + 3)  =  0.267730531659313
 * 
*/

//I calculate in advance the coefficients for the filter
int Q = 12; // I can choose Q freely 
//I want a cutoff of 50Hz at a sampling frequency of 1000 Hz. So I calculate the corresponding alpha with formula above
//For cutoff 50 Hz, alpha = 0.267730531659313; then alpha_scaled = round(2^Q*alpha) = round(2^12* 0.267730531659313) = round(4096* 0.267730531659313) = 1097
int alpha_scaled = 1097;
// N ~= -log_2(alpha) = -log_2(0.267730531659313) = 1.901146423365515 ~ 2  
int N = 2;




int low_pass_EMA_fixed_point(int x, int y, int alpha_scaled, int Q, int N){
  y += (alpha_scaled * (x-y)) >> (Q-N);
  return y >> N;
}



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

void setup() {
  Serial.begin(115200);
  analog_input0_filtered_fixed_point = analogRead(analog_pin); //set EMA y for t=1 
}

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

        analog_input0 = analogRead(analog_pin);  
        analog_input0_filtered_fixed_point = low_pass_EMA_fixed_point(analog_input0, analog_input0_filtered_fixed_point, alpha_scaled, Q, N); 

        //Check the original and filtered signals with the serial plotter
        Serial.print(analog_input0); 
        Serial.print(" ");
        Serial.println(analog_input0_filtered_fixed_point);  
  }
}

See above a plot with the serial plotter.

Any idea about how to solve the offset problem? What am I doing wrong? Any suggestion about how to solve the general problem of a fixed point exponential moving average filter with a specific cutoff frequency?

Answer 1

You should scale your analog input to match your desired resolution, it's better to have bigger resolution (specially when dealing with IIR filters), I don't know the reason of returning y>>N, probably the scaling problem is there, I would implemented it like this:

int low_pass_EMA_fixed_point(int x, int y, int alpha_scaled, int Q, int N)
{
  int32_t y_temp;
  y_temp = (alpha_scaled*(x<<N))+((1<<Q)-alpha_scaled)*y;
  return (int)(y_temp >> Q);
}

Also a good practice is test your filter with an impulse to check if it was implemented correctly before using the actual data.

EDIT: tested this program with onlinegdb.com (the printf are used to copy the result into octave):

#include <stdio.h>
#include <stdlib.h>

#define NUM 1024

int low_pass_EMA_fixed_point(int x, int y, int alpha_scaled, int Q, int N);

int main()
{
    int alpha_scaled = 1097; //Q12
    int Q = 12;
    int N = 2;
    static int x[NUM];
    int y = 0;

    x[0] = 1024; //input in Q10
    printf("y=[");
    for(int i=0;i<NUM;i++)
    {
        y = low_pass_EMA_fixed_point(x[i], y, alpha_scaled, Q, N);
        if(i == NUM-1)
        {
            printf("%d];",y); //print the last output
        }
        else
        {
            printf("%d,",y);
        }
    }
    return 0;
}

int low_pass_EMA_fixed_point(int x, int y, int alpha_scaled, int Q, int N)
{
  int32_t y_temp;
  y_temp = (alpha_scaled*(x<<N))+(((1<<Q)-alpha_scaled)*y);

  return (int)(y_temp>>Q);
}

And I'm getting this frequency response in octave:

Answer 2

Try analog_input0 << Q, just like what https://www.embeddedrelated.com/showarticle/779.php did:

x1q = np.round(x1*(1<<Q)).astype(int)