0
$\begingroup$

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);  
  }
}

Serial plotter screenshot

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?

$\endgroup$

migrated from dsp.meta.stackexchange.com May 20 at 19:06

This question came from our discussion, support, and feature requests site for practitioners of the art and science of signal, image and video processing.

1
$\begingroup$

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: enter image description here

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – L_T May 21 at 10:15
  • $\begingroup$ 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 $\endgroup$ – Hkxs May 21 at 15:06
  • $\begingroup$ Thanks. Using your code the resulting signal does not really seem to be affected by any filtering....not really sure why. $\endgroup$ – L_T May 21 at 18:01
  • $\begingroup$ 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 $\endgroup$ – Hkxs May 21 at 18:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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