作者： Mustahsin Zarif

之前我們在《信號處理簡介》一文中已經(jīng)見過了兩類濾波器：有限脈沖響應 (FIR) 濾波器和無限脈沖響應 (IIR) 濾波器。我們看到了移動平均濾波器如何同時以 FIR 和 IIR 形式表示，但如果它們相互比較又各具哪些優(yōu)勢呢？

回顧我之前博客中的例子，我們可以將 FIR 濾波器展開為如下形式：

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

在這里，我們看到我們需要：

1. 5 次乘法和
2. 4 次求和運算。

乘法運算的計算成本特別高。因此，如果我們再次查看 IIR 形式，我們會發(fā)現(xiàn)它只需要：

1. 3 次乘法和
2. 2 次求和運算

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

這大大降低了計算成本！這對于單片機等嵌入式設備來說非常好，因為它們在每個離散時間步驟上執(zhí)行計算所消耗的資源更少。

例如，當我對采用 FIR 和 IIR 形式的 11 點移動平均濾波器使用 Python 函數(shù) ‘time.time’ 時，所有參數(shù)（窗口大小、采樣率、樣本大小等）相同，我分別得到以下運行時間結果：51 ms、27 ms。

離散時間IIR 濾波器示例

現(xiàn)在我們已經(jīng)了解了為什么 IIR 濾波器在單片機上表現(xiàn)更好了，讓我們看一個使用 Arduino UNO 和 DFRobot MPU6050 慣性測量單元 (IMU) 的示例項目（圖 1）。我們將對 IMU 數(shù)據(jù)應用指數(shù)移動平均 (EMA) 濾波器，以查看原始數(shù)據(jù)和平滑數(shù)據(jù)之間的差異。

圖 1：MPU6050 與 Arduino Uno 之間的連接框圖。（圖片來源：Mustahsin Zarif）

圖 2：MPU6050 與 Arduino Uno 之間的連接。（圖片來源：Mustahsin Zarif）

指數(shù)移動平均濾波器具有遞歸形式：

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

它之所以是遞歸，是因為我們測量的任何當前輸出也取決于先前的輸出；即系統(tǒng)具有記憶。

常數(shù) alpha () 決定了我們想要賦予當前輸入相對于先前輸出多大的權重。為了清楚起見，讓我們展開方程得到：

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

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

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

我們看到，alpha 越大，當前輸入對當前輸出的影響就越大。這是好事，因為如果系統(tǒng)在不斷演進，過去的值不能代表當前的系統(tǒng)。另一方面，如果系統(tǒng)突然發(fā)生瞬間異常變化，情況就會變得很糟糕；在這種情況下，我們希望我們的輸出能夠遵循之前輸出所遵循的趨勢。

事不宜遲，現(xiàn)在讓我們看看 EMA 濾波器的代碼是如何用于 MPU6050 的。

EMA 濾波器代碼：

副本#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 >

當我們運行此代碼并檢查串口繪圖儀時，我們可以看到 x、y 和 z 軸方向加速度的成對粗糙和平滑線條，其中窗口大小為 11 和 alpha 值為 0.2（圖 3 至 5）。

圖 3：x 方向的原始加速度值和濾波后的加速度值。（圖片來源：Mustahsin Zarif）

圖 4：y 方向的原始加速度值和濾波后的加速度值。（圖片來源：Mustahsin Zarif）

圖 5：z 方向的原始加速度值和濾波后的加速度值。（圖片來源：Mustahsin Zarif）

讓代碼的智能化更進一步

我們現(xiàn)在知道，與 FIR 濾波器相比，IIR 濾波器更適合用作控制器，因為所需的求和和乘法計算明顯較少。然而，當我們實現(xiàn)這段代碼時，執(zhí)行的計算并不只有求和和乘法：每當有新的時間樣本進入時，我們都必須移動樣本，而這個過程在后臺需要計算能力。因此，我們可以借助循環(huán)緩沖區(qū)，而不是在每個采樣時間間隔移動所有樣本。

我們的做法是：用一個指針來記住傳入的數(shù)據(jù)樣本的索引。然后，每次指針指向緩沖區(qū)中的最后一個元素時，它接下來都會指向緩沖區(qū)的第一個元素，新數(shù)據(jù)將替換之前存儲在這里的數(shù)據(jù)，因為這是現(xiàn)在我們不再需要的最舊數(shù)據(jù)（圖 6）。因此，這種方法允許我們跟蹤緩沖區(qū)中最舊的樣本并替換該樣本，而不必每次都移動樣本以將新數(shù)據(jù)放入數(shù)組的最后一個元素中。

圖 6：循環(huán)緩沖區(qū)示例圖。（圖片來源：Mustasin Zafir）

這是使用循環(huán)緩沖區(qū)的 EMA 濾波器實現(xiàn)的代碼。您能嘗試對陀螺儀而不是對加速計運行這段代碼嗎？也可以嘗試使用不同的系數(shù)！

使用循環(huán)緩沖區(qū)代碼的 EMA 濾波器：

副本#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 >

結語

在這篇博客中，我們討論了 IIR 和 FIR 濾波器之間的區(qū)別，重點討論了它們的計算效率。通過從 FIR 到 IIR 所需運算次數(shù)減少這一小例子，我們可以想象當應用規(guī)?；瘯r IIR 濾波器的效率會有多高，這對于硬件能力有限的實時應用非常重要。

我們還研究了一個使用 Arduino Uno 和 MPU6050 IMU 的示例項目，其中我們部署了一個指數(shù)移動平均濾波器來降低傳感器數(shù)據(jù)中的噪聲，同時仍然捕捉底層信號行為。最后，為了提高效率，我們提供了一個更智能的示例代碼，即采用循環(huán)緩沖區(qū)而不是在每個時間間隔移動數(shù)據(jù)。

在下一篇博客中，我們將利用 [Red Pitaya] 的 FPGA 功能來實現(xiàn)一個 4 抽頭 FIR 濾波器數(shù)字電路！