指数移動平均（EMA）フィルタの周波数特性

2 min read

指数移動平均 (Exponential Moving Average; EMA)

指数移動平均（EMA）は、時系列データを平滑化する代表的な手法の一つです。現在の観測値と1ステップ前の平滑化された値に重み付けをして足し合わせることで、新しい平滑化値を計算します。

$$ y_t = (1 - \alpha) y_{t-1} + \alpha x_t $$

ここで、

  • $y_t$: 時刻 $t$ におけるEMAの値
  • $x_t$: 時刻 $t$ における元の観測値
  • $\alpha$: 平滑化係数（または平滑化率、$0 < \alpha \le 1$）

この式は、最近のデータに大きな重みを置きつつ、古いデータも完全に切り捨てずに考慮するという特徴があります。元の記事の式は $\beta$ を使っていますが、一般的に $\alpha$ が使われることが多いので、ここでは $\alpha$ に統一します。元の記事の $\beta$ は、ここでいう $1-\alpha$ に相当します。

周波数特性：ゲイン特性と位相特性

フィルタの周波数特性は、入力信号の周波数成分がフィルタによってどのように変化するかを示します。これは、フィルタの伝達関数 $G(s)$ を用いて分析されます。

上記のEMAの式をZ変換（離散時間システムの解析に用いられる）すると、伝達関数 $G(z)$ は次のようになります。

$$ Y(z) = (1 - \alpha) z^{-1} Y(z) + \alpha X(z) $$ $$ G(z) = \frac{Y(z)}{X(z)} = \frac{\alpha}{1 - (1 - \alpha) z^{-1}} $$

ここで、$z = e^{j\omega T}$ （$T$ はサンプリング周期、ここでは $T=1$ と仮定）とすると、周波数応答 $G(j\omega)$ が得られます。

$$ G(j\omega) = \frac{\alpha}{1 - (1 - \alpha) e^{-j\omega}} $$

この周波数応答から、ゲイン特性 $|G(j\omega)|$ と位相特性 $\angle G(j\omega)$ を導出できます。

ゲイン特性

$$ |G(j\omega)| = \frac{\alpha}{\sqrt{1 - 2(1 - \alpha) \cos \omega + (1 - \alpha)^2}} $$

ゲイン特性

位相特性

$$ \angle G(j\omega) = \arctan\left(\frac{(1 - \alpha) \sin \omega}{1 - (1 - \alpha) \cos \omega}\right) $$ 元の記事の式は符号が逆になっていますが、これは$\arctan(Y/X)$の計算において$X$と$Y$の符号の組み合わせを考慮する必要があるためです。一般的な定義では上記のとおりです。

位相特性

考察:

  • $\alpha$ が小さい場合（例: $\alpha=0.1$）：フィルタの周波数帯域が狭くなり、より強い平滑化効果が得られます。しかし、位相遅れが大きくなります。
  • $\alpha$ が大きい場合（例: $\alpha=0.9$）：平滑化能力は低下しますが、位相遅れは小さくなります。

これは、平滑化の度合いと応答速度（位相遅れ）の間にトレードオフがあることを示しています。

プログラム

以下のPythonコードは、EMAフィルタのゲイン特性と位相特性を計算し、プロットするものです。

import math
import numpy as np # numpyをインポート
import matplotlib.pyplot as plt

# 平滑化係数αのリスト
ALPHAS = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
# 周波数範囲 (rad/s)
FREQ_MAX = 300 # 300 rad/sまで
NUM_POINTS = 500 # プロット点数

def get_gain(alpha, omega):
    """EMAフィルタのゲインを計算する"""
    denominator = np.sqrt(1 - 2 * (1 - alpha) * np.cos(omega) + (1 - alpha)**2)
    return alpha / denominator

def get_phase(alpha, omega):
    """EMAフィルタの位相を計算する (ラジアン)"""
    numerator = (1 - alpha) * np.sin(omega)
    denominator = 1 - (1 - alpha) * np.cos(omega)
    return np.arctan2(numerator, denominator) # atan2を使用することで正しい象限の角度が得られる

# プロット用のデータ格納リスト
frequencies = [[] for _ in range(len(ALPHAS))]
gains = [[] for _ in range(len(ALPHAS))]
phases = [[] for _ in range(len(ALPHAS))]

# 各α値について計算
for alpha_idx, alpha_val in enumerate(ALPHAS):
    for i in range(NUM_POINTS):
        omega = i * (np.pi / NUM_POINTS) # 0からπまでをプロット (ナイキスト周波数まで)
        # omega = i * (FREQ_MAX / NUM_POINTS) # もしFREQ_MAXまでプロットしたい場合
        
        frequencies[alpha_idx].append(omega)
        gains[alpha_idx].append(get_gain(alpha_val, omega))
        phases[alpha_idx].append(np.degrees(get_phase(alpha_val, omega))) # 位相を度に変換

# ゲイン特性のプロット
plt.figure(figsize=(10, 6))
plt.xlabel('Frequency (rad/s)')
plt.ylabel('Gain') # dB表示にする場合は20*log10(Gain)
plt.title('EMA Filter Gain Characteristics')
for alpha_idx, alpha_val in enumerate(ALPHAS):
    plt.plot(frequencies[alpha_idx], gains[alpha_idx], label=f"alpha={alpha_val}")
plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')
plt.grid(True)
plt.tight_layout()
plt.show()

# 位相特性のプロット
plt.figure(figsize=(10, 6))
plt.xlabel('Frequency (rad/s)')
plt.ylabel('Phase (degrees)')
plt.title('EMA Filter Phase Characteristics')
for alpha_idx, alpha_val in enumerate(ALPHAS):
    plt.plot(frequencies[alpha_idx], phases[alpha_idx], label=f"alpha={alpha_val}")
plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')
plt.grid(True)
plt.tight_layout()
plt.show()

参考