信号処理：指数移動平均（EMA）フィルタ

前回までの信号処理入門では、有限インパルス応答（FIR）と無限インパルス応答（IIR）の2種類のフィルタを見てきました。移動平均フィルタがFIR型とIIR型の両方で表現できることを確認しましたが、それぞれにはどのような利点があるのでしょうか。

前回のブログの例を振り返り、FIRフィルタを展開すると、次のような形式になります。

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

ここで、次のことが必要であることがわかります。

1. 5回の乗算および
2. 4回の加算

乗算は特に計算負荷の大きい操作です。したがって、IIRの形式をもう一度見てみると、必要な演算は次の通りです。

1. 3回の乗算および
2. 2回の加算

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

これにより、計算コストが大幅に削減されます。これは、マイクロコントローラのような組み込みデバイスにとっては、好都合です。なぜなら各離散時間ステップで計算を実行するために消費するリソースが少なくて済むからです。

例として、Pythonの関数「time.time」を使用して、FIR形式とIIR形式の11点移動平均フィルタを実行し、すべてのパラメータ（ウィンドウサイズ、サンプルレート、サンプルサイズなど）を同じに設定した場合、それぞれ51ms、27msの実行時間が得られました。

離散時間IIRフィルタの例

IIRフィルタがマイクロコントローラでより優れた性能を発揮する理由について理解できたところで、Arduino UNOとDFRobot MPU6050慣性計測ユニット（IMU）を使用したプロジェクトの例を見てみましょう（図1）。IMUデータに指数移動平均（EMA）フィルタを適用し、生データと平滑化データの違いを確認します。

図1：MPU6050とArduino Unoの接続ブロック図。(画像提供：Mustahsin Zarif氏）

図2：MPU6050とArduino Unoの接続。（画像提供：Mustahsin Zarif氏）

指数移動平均フィルタは次のような再帰形式です。

y[n] = α*x[n]+（1-α）*y[n-1]

再帰形式である理由は、現在の出力が以前の出力にも依存するためです。つまり、システムにはメモリがあるということです。

定数アルファ（α）は、現在の入力に与える重みを、以前の出力に対してどれだけ設定するかを決定します。わかりやすくするために、式を次のように展開します。

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]

αが大きいほど、現在の入力が現在の出力に与える影響が大きくなります。システムが進化している場合、はるか過去の値は現在のシステムを示すものとはならないため、これは望ましいことです。一方、たとえば、システムが正常な状態から突然、一次的に変化した場合、これは望ましいことではありません。この場合、出力は、以前の出力の傾向に従うようにしたいからです。

それでは、MPU6050用のEMAフィルタのコードがどのように機能するか確認しましょう。

EMAフィルターコード：

コピー#include<wire.h>
 #include<mpu6050.h>

 MPU6050 mpu;

#define BUFFER_SIZE 11 // ウィンドウサイズ

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 接続に失敗しました！");
 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；
 }

//バッファに格納された加速度値を使用してEMAを計算
 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>

このコードを実行し、シリアルプロッタを確認すると、ウィンドウサイズ11、α値0.2を使用した場合、x、y、z軸の加速度に対して、粗い線と滑らかな線がペアになっているのがわかります（図3～図5）。

図3：X方向の生加速度値とフィルタ処理後の加速度値。（画像提供：Mustahsin Zarif氏）

図4：Y方向の生加速度値とフィルタ処理後の加速度値。（画像提供：Mustahsin Zarif氏）

図5：z方向の生加速度値とフィルタ処理後の加速度値。（画像提供：Mustahsin Zarif氏）

コードをさらにスマートにする

IIRフィルタがFIRフィルタに比べてコントローラに適している理由は、必要な加算と乗算の計算回数が大幅に少ないためです。しかし、このコードを実装すると、加算と乗算だけが実行されるわけではありません。新しい時間サンプルが入るたびにサンプルをシフトする必要があり、この処理は、内部で計算能力を必要とします。したがって、すべてのサンプルを各サンプリング時間間隔でシフトする代わりに、サーキュラバッファを活用することができます。

ここでは次のように、入力するデータサンプルのインデックスを記憶するポインタを用意します。それから、ポインタがバッファの最後の要素を指すたびに、次にポインタはバッファの最初の要素を指し、新しいデータで以前ここに格納されていたデータが置き換えられます。これは、もう不要になった最も古いデータだからです（図6）。この方法により、バッファ内の最も古いサンプルを追跡し、新しいデータを配列の最後の要素に格納するためにサンプルを毎回シフトする必要なく、そのサンプルを置き換えることができます。

図6：サーキュラバッファの例示図。（画像提供：Mustahsin Zarif氏）

これは、サーキュラバッファを使用したEMAフィルタのコードの例です。加速度センサの代わりにジャイロセンサで実行してみてください。また、係数もいろいろ試してみてください。

サーキュラバッファコードを使用したEMAフィルタ：

コピー#include<wire.h>

 #include<mpu6050.h>

 MPU6050 mpu;

#define BUFFER_SIZE 11 // ウィンドウサイズ

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 接続に失敗しました！");

 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; //サーキュラバッファの実装

 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フィルタの違いについて、特に計算効率に焦点を当てて説明しました。FIRからIIRへの移行で必要な演算数の削減という小さな例を取り上げることで、IIRフィルタの効率がアプリケーションのスケールアップ時にどれだけ向上するかを想像できます。これは、限られたハードウェアリソースで動作するリアルタイムアプリケーションにおいて重要です。

また、Arduino UnoとMPU6050 IMUを使用したサンプルプロジェクトでは、センサデータのノイズを低減しつつ、元の信号の挙動を捉えるために、指数移動平均フィルタを実装しました。最後に、効率向上の観点から、データのシフトを毎時間間隔で行う代わりにサーキュラバッファを採用したスマートなコードの例も紹介しました。

次回のブログでは、Red PitayaのFPGA機能を使用して、4タップFIRフィルタのデジタル回路を実装します。