Fixed-point Exponential Moving Average filter with given cutoff frequency

Projects Programming
1

Hello,
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 filters - Exponential moving average cut-off frequency - Signal Processing Stack Exchange 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: Ten Little Algorithms, Part 2: The Single-Pole Low-Pass Filter - Jason Sachs but I am getting an offset that I can't explain,
Here my 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 attached 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?

error.png
error.png1684×1326 94.5 KB

2

If you're using Arduino Uno, the size of integer is 16 bit, with a range -32768 to 32767.

Your alpha_scaled constant is 1097, 11 bit. If (x-y) size is >= 5 bit, the signed integer multiplication overflows.

For example, if (x-y) = 32 then alpha_scaled * (x-y) = 1097 * 32 = 35104 > 32767 --> overflow.

So, if you're on Arduino Uno, try using long integers, 32 bit.

3

Many thanks! Of course this was an error. However, even substituting all int with int32_t the problem persists. Here the updated code:

#define PERIOD_MICROSECS 1000

static unsigned long lastRead = 0;

int32_t analog_pin = 0;   // sensor connected to analog pin 0
int32_t analog_input0 = 0;  // variable to store the read value; 
int32_t 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
int32_t 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
int32_t alpha_scaled = 1097;
// N ~= -log_2(alpha) = -log_2(0.267730531659313) = 1.901146423365515 ~ 2  
int32_t N = 2;




int32_t low_pass_EMA_fixed_point(int32_t x, int32_t y, int32_t alpha_scaled, int32_t Q, int32_t 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);  
  }
}

Any idea about how to solve the issue?

4

Well

y += (alpha_scaled * (x-y)) >> (Q-N);

in my opinion, that instruction can't work.

In input to function low_pass_EMA_fixed_point(), y is an integer.

But

(alpha_scaled * (x-y))

is a fixed point number with Q binary digits after decimal point, and

(alpha_scaled * (x-y)) >> (Q-N)

is a fixed point number with N binary digits after decimal point.

So

y += (alpha_scaled * (x-y)) >> (Q-N);

is a misaligned sum of an integer with a fixed point number.

The right way to do this, is to align y before sum

int32_t low_pass_EMA_fixed_point(int32_t x, int32_t y, int32_t alpha_scaled, int32_t Q, int32_t N)
{  

 int32_t scaled_dy = (alpha_scaled * (x-y)) >> (Q-N); // scaled_dy is a fixed point number with N binary digit after decimal point

 scaled_y = y << N; // scaled_y is now a fixed point number with N binary digit (all 0) after decimal point

 scaled_y += scaled_dy; // aligned sum of fixed point numbers

 return(scaled_y >> N); // return an integer, converting the fixed point with N binary digit after decimal point
}

But this way I don't understand the N thing meaning: the function low_pass_EMA_fixed_point() always skips the N extra binary digit without incrementing filter precision, I think something is missed.

Maybe you must preserve scaled_y value, so the intended function is:

int32_t scaled_y = 0; // static fixed point variable for running sum with N binary digit after decimal point

int32_t low_pass_EMA_fixed_point(int32_t x, int32_t y, int32_t alpha_scaled, int32_t Q, int32_t N)
{  
 scaled_y += (alpha_scaled * (x-y)) >> (Q-N); // scaled_y is a fixed point number with N binary digit after decimal point

 return(scaled_y >> N); // return an integer, converting the fixed point with N binary digit after decimal point
}