作者： Mustahsin Zarif

之前我們?cè)凇?a target="_blank">信號(hào)處理簡(jiǎn)介》一文中已經(jīng)見過(guò)了兩類濾波器：有限脈沖響應(yīng) (FIR) 濾波器和無(wú)限脈沖響應(yīng) (IIR) 濾波器。我們看到了移動(dòng)平均濾波器如何同時(shí)以 FIR 和 IIR 形式表示，但如果它們相互比較又各具哪些優(yōu)勢(shì)呢？

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

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

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

1. 5 次乘法和
2. 4 次求和運(yùn)算。

乘法運(yùn)算的計(jì)算成本特別高。因此，如果我們?cè)俅尾榭?IIR 形式，我們會(huì)發(fā)現(xiàn)它只需要：

1. 3 次乘法和
2. 2 次求和運(yùn)算

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

這大大降低了計(jì)算成本！這對(duì)于單片機(jī)等嵌入式設(shè)備來(lái)說(shuō)非常好，因?yàn)樗鼈冊(cè)诿總€(gè)離散時(shí)間步驟上執(zhí)行計(jì)算所消耗的資源更少。

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

離散時(shí)間IIR 濾波器示例

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

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

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

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

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

它之所以是遞歸，是因?yàn)槲覀儨y(cè)量的任何當(dāng)前輸出也取決于先前的輸出；即系統(tǒng)具有記憶。

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

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 越大，當(dāng)前輸入對(duì)當(dāng)前輸出的影響就越大。這是好事，因?yàn)槿绻到y(tǒng)在不斷演進(jìn)，過(guò)去的值不能代表當(dāng)前的系統(tǒng)。另一方面，如果系統(tǒng)突然發(fā)生瞬間異常變化，情況就會(huì)變得很糟糕；在這種情況下，我們希望我們的輸出能夠遵循之前輸出所遵循的趨勢(shì)。

事不宜遲，現(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 >

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

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

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

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

讓代碼的智能化更進(jìn)一步

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

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

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

這是使用循環(huán)緩沖區(qū)的 EMA 濾波器實(shí)現(xiàn)的代碼。您能嘗試對(duì)陀螺儀而不是對(duì)加速計(jì)運(yù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 >

結(jié)語(yǔ)

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

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

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