# Fixed-point Exponential Moving Average filter with given cutoff frequency

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) {

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?

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:

• Thanks. Unfortunately your solution does not work. Actually regarding the resolution I changed all int to int32_t to have a bigger resolution. Still the problem persists... the filter does not work. Any idea? – L_T May 21 at 10:15
• Hi sorry, I assumed that your int was 32 bits, also I have some issues with the signs. Also, with this code you'll have to scale one of the signals that you're sending, either x<<N or y>>N since the output will be in Q12. I'll update the code – Hkxs May 21 at 15:06
• Thanks. Using your code the resulting signal does not really seem to be affected by any filtering....not really sure why. – L_T May 21 at 18:01
• I just tested the code in gdbonline.com with an impulse response, I missed the parenthesis around (1<<Q) and I'm getting the FFT, it seems to be filtering the high frequencies – Hkxs May 21 at 18:45