Signal Processing: Exponentially Moving Average (EMA) Filter

Previously, in Introduction to Signal Processing, we've seen the two classes of filters: Finite Impulse Response (FIR) and Infinite Impulse Response (IIR). We saw how the moving average filter can be expressed in both a FIR and IIR form, but what are the benefits of one over the other?

Looking back at the example from my previous blog, the FIR filter expanded has the form:

y[5] = (x[5]+x[4]+x[3]+x[2]+x[1]+x[0])/5,

Here, we see that we require:

1. 5 multiplication and
2. 4 summation operations.

Multiplication operations are particularly computationally expensive. Therefore, if we look at the IIR form again, we see that it requires only:

1. 3 multiplication and
2. 2 summation operations

y[6]=(x[6]+y[5]-x[1])/5

This reduces the computation cost significantly! This is good for embedded devices such as microcontrollers since they spend less resources at every discrete time step to perform calculations.

As an example, when I use the Python function ‘time.time’ for the 11 point moving average filter in FIR and IIR form, with all the parameters (window size, sample rate, sample size, etc.) the same, I get the following runtime results respectively: 51 ms, 27 ms.

Discrete time IIR filter example

Now that we have an intuition as to why IIR filters perform better on microcontrollers, let's see an example project using an Arduino UNO and a DFRobot MPU6050 inertial measurement unit (IMU) (Figure 1). We will apply the exponentially moving average (EMA) filter to the IMU data to see differences between raw and smoothed data.

Figure 1: Block diagram connection between the MPU6050 and the Arduino Uno. (Image source: Mustahsin Zarif)

Figure 2: Connection between the MPU6050 and the Arduino Uno. (Image source: Mustahsin Zarif)

The exponentially moving average filter is of the recursive form:

y[n] = α*x[n] + (1- α)*y[n-1]

It is recursive because any current output that we are measuring depends on the previous outputs as well; i.e., the system has memory.

The constant alpha (𝞪) determines how much weight we want to give to the current input as opposed to the previous outputs. For clarity, lets expand the equation to get:

y[n] = α*x[n] + (1- α )*(α*x[n−1]+(1−α)*y[n−2])

y[n] = α*x[n] + (1- α )*x[n−1]+α*(1−α)2*x[n−2])+ ...

y[n] = k=0nα*(1−α)k*x[n−k]

We see that the greater the alpha, the more the current input affects the current output. This is good since if the system is evolving, values far in the past are not as representative of the current system. On the other hand, this would be bad if, for example, there is a sudden, momentary change in the system from normal; in this case, we would like our output to follow the trend that the previous outputs were following.

Now without further ado, let's see how the code for an EMA filter would work for the MPU6050.

EMA Filter Code:

Copy#include <wire.h>
#include <mpu6050.h>

MPU6050 mpu;

#define BUFFER_SIZE 11  // Window size

float accelXBuffer[BUFFER_SIZE];
float accelYBuffer[BUFFER_SIZE];
float accelZBuffer[BUFFER_SIZE];
int bufferCount = 0;  

void setup() {
  Serial.begin(115200);
  Wire.begin();

  mpu.initialize();

  if (!mpu.testConnection()) {
    Serial.println("MPU6050 connection failed!");
    while (1);
  }

  int16_t ax, ay, az;
  for (int i = 0; i < BUFFER_SIZE; i++) {
    mpu.getMotion6(&ax, &ay, &az, NULL, NULL, NULL);
    accelXBuffer[i] = ax / 16384.0;
    accelYBuffer[i] = ay / 16384.0;
    accelZBuffer[i] = az / 16384.0;
  }
  bufferCount = BUFFER_SIZE;
}

void loop() {
  int16_t accelX, accelY, accelZ;

  mpu.getMotion6(&accelX, &accelY, &accelZ, NULL, NULL, NULL);

  float accelX_float = accelX / 16384.0;
  float accelY_float = accelY / 16384.0;
  float accelZ_float = accelZ / 16384.0;

  if (bufferCount < BUFFER_SIZE) {
    accelXBuffer[bufferCount] = accelX_float;
    accelYBuffer[bufferCount] = accelY_float;
    accelZBuffer[bufferCount] = accelZ_float;
    bufferCount++;
  } else {
    for (int i = 1; i < BUFFER_SIZE; i++) {
      accelXBuffer[i - 1] = accelXBuffer[i];
      accelYBuffer[i - 1] = accelYBuffer[i];
      accelZBuffer[i - 1] = accelZBuffer[i];
    }
    accelXBuffer[BUFFER_SIZE - 1] = accelX_float;
    accelYBuffer[BUFFER_SIZE - 1] = accelY_float;
    accelZBuffer[BUFFER_SIZE - 1] = accelZ_float;
  }

//calculate EMA using acceleration values stored in buffer
  float emaAccelX = accelXBuffer[0];
  float emaAccelY = accelYBuffer[0];
  float emaAccelZ = accelZBuffer[0];
  float alpha = 0.2;

  for (int i = 1; i < bufferCount; i++) {
    emaAccelX = alpha * accelXBuffer[i] + (1 - alpha) * emaAccelX;
    emaAccelY = alpha * accelYBuffer[i] + (1 - alpha) * emaAccelY;
    emaAccelZ = alpha * accelZBuffer[i] + (1 - alpha) * emaAccelZ;
  }

  Serial.print(accelX_float); Serial.print(",");
  Serial.print(accelY_float); Serial.print(",");
  Serial.print(accelZ_float); Serial.print(",");
  Serial.print(emaAccelX); Serial.print(",");
  Serial.print(emaAccelY); Serial.print(",");
  Serial.println(emaAccelZ);

  delay(100);
}
</mpu6050.h></wire.h>

When we run this code and check the serial plotter, we can see rough and smooth lines in pairs for accelerations in the x, y, and z axes, using a window size of 11 and an alpha value of 0.2 (Figure 3 to 5).

Figure 3: Raw and filtered acceleration values in the x-direction. (Image source Mustahsin Zarif)

Figure 4: Raw and filtered acceleration values in the y-direction. (Image source Mustahsin Zarif)

Figure 5: Raw and filtered acceleration values in the z-direction. (Image source Mustahsin Zarif)

Making the code one step smarter

We now have an idea of how IIR filters are better for controllers compared to FIR filters because of the significantly fewer summation and multiplication calculations required. However, when we implement this code, summation and multiplication aren’t the only calculations being performed: we have to shift the samples every time a new time sample comes in, and this process, under the hood, requires compute power. Therefore, instead of shifting all of the samples at every sampling time interval, we can employ the help of circular buffers.

Here is what we do: we have a pointer that remembers the index of the data sample that came in. Then, every time the pointer points to the last element in the buffer, it points to the first element of the buffer next, and the new data replaces the data that was stored here before, since this is now the oldest data that we do not need anymore (Figure 6). Consequently, this method allows us to keep track of the oldest sample in the buffer and replace that one without having to shift samples every time to put the new data in the last element of the array.

Figure 6: Example illustration of a circular buffer. (Image source: Mustasin Zafir)

This is what the code looks like for an EMA filter implementation using circular buffers. Can you try to run this for a gyroscope instead of an accelerometer? Play around with the coefficients too!

EMA Filter Using a Circular Buffer Code:

Copy#include <wire.h>

#include <mpu6050.h>

MPU6050 mpu;

#define BUFFER_SIZE 11  // Window size

float accelXBuffer[BUFFER_SIZE];

float accelYBuffer[BUFFER_SIZE];

float accelZBuffer[BUFFER_SIZE];

int bufferIndex = 0;  

void setup() {

  Serial.begin(115200);

  Wire.begin();
 

  mpu.initialize();


  if (!mpu.testConnection()) {

    Serial.println("MPU6050 connection failed!");

    while (1);

  }

  int16_t ax, ay, az;

  for (int i = 0; i < BUFFER_SIZE; i++) {

    mpu.getMotion6(&ax, &ay, &az, NULL, NULL, NULL);

    accelXBuffer[i] = ax / 16384.0;

    accelYBuffer[i] = ay / 16384.0;

    accelZBuffer[i] = az / 16384.0;

  }

}

void loop() {

  int16_t accelX, accelY, accelZ;

  mpu.getMotion6(&accelX, &accelY, &accelZ, NULL, NULL, NULL);

  float accelX_float = accelX / 16384.0;

  float accelY_float = accelY / 16384.0;

  float accelZ_float = accelZ / 16384.0;

  accelXBuffer[bufferIndex] = accelX_float;

  accelYBuffer[bufferIndex] = accelY_float;

  accelZBuffer[bufferIndex] = accelZ_float;

  bufferIndex = (bufferIndex + 1) % BUFFER_SIZE; //circular buffer implementation 

  float emaAccelX = accelXBuffer[bufferIndex];

  float emaAccelY = accelYBuffer[bufferIndex];

  float emaAccelZ = accelZBuffer[bufferIndex];

  float alpha = 0.2;

  for (int i = 1; i < BUFFER_SIZE; i++) {

    int index = (bufferIndex + i) % BUFFER_SIZE;

    emaAccelX = alpha  accelXBuffer[index] + (1 - alpha)  emaAccelX;

    emaAccelY = alpha  accelYBuffer[index] + (1 - alpha)  emaAccelY;

    emaAccelZ = alpha  accelZBuffer[index] + (1 - alpha)  emaAccelZ;

  }

  Serial.print(accelX_float); Serial.print(",");

  Serial.print(emaAccelX); Serial.print(",");

  Serial.print(accelY_float); Serial.print(",");

  Serial.print(emaAccelY); Serial.print(",");

  Serial.print(accelZ_float); Serial.print(",");

  Serial.println(emaAccelZ);

  delay(100);

}
</mpu6050.h></wire.h>

Summary

In this blog, we discussed the difference between IIR and FIR filters with an emphasis on their computational efficiencies. By taking a small example of the reduction in number of operations required from FIR to IIR, we can imagine how efficient IIR filters will be when applications are scaled, which is important for real-time applications on limited hardware power.

We also took a look at an example project using an Arduino Uno and MPU6050 IMU, where we deployed an exponentially moving average filter to reduce noise in sensor data while still capturing the underlying signal behavior. Finally, in the interest of efficiency, we also saw an example of smarter code by employing circular buffers instead of shifting data at every time interval.

In the next blog, we will use Red Pitaya’s FPGA functionality to implement a 4 tap FIR filter digital circuit!