移動平均フィルタの種類と比較:SMA・WMA・EMAのPython実装

移動平均フィルタ(SMA / WMA / EMA)の理論・設計・比較・Python実装(numpy.convolve, pandas.rolling, scipy.signal.lfilter)。畳み込みによる平滑化の数理、ノイズ除去性能・群遅延・周波数応答(sinc関数)の比較から、時系列データに最適な平滑化フィルタの選び方を整理します。

移動平均(Moving Average)は、時系列データの平滑化に最もよく使われるフィルタリング手法です。ノイズ除去、トレンド抽出、センサデータの前処理など幅広い場面で活用されます。

この記事では、代表的な3つの移動平均フィルタを数理的に定義し、Python実装と性能比較を行います。単に特性を並べるだけでなく、なぜその特性が生まれるのかを周波数応答・群遅延の導出から示し、実際にコードを実行して得られた具体的な数値・図表で定量的に比較します。

単純移動平均 (Simple Moving Average, SMA)

定義

SMAは、直近 \(N\) 個のデータ点の算術平均です。

\[ y_t = \frac{1}{N} \sum_{i=0}^{N-1} x_{t-i} \tag{1} \]

すべてのデータ点に等しい重み \(1/N\) を与えます。\(N\) が**窓幅(ウィンドウサイズ)**で、大きいほど平滑化が強くなります。

周波数応答の導出

SMAのインパルス応答 \(h[n] = 1/N\) (\(n = 0, \dots, N-1\) )に対する周波数応答は、有限等比級数として書けます。

\[ H(e^{j\omega}) = \frac{1}{N} \sum_{i=0}^{N-1} e^{-j i \omega} = \frac{1}{N} \cdot \frac{1 - e^{-jN\omega}}{1 - e^{-j\omega}} \tag{1a} \]

この式の分子・分母をそれぞれ \(e^{-jN\omega/2}\) 、\(e^{-j\omega/2}\) でくくり出すと、ディリクレ核(Dirichlet kernel)と呼ばれる形に整理できます。

\[ H(e^{j\omega}) = \frac{1}{N} \cdot \frac{\sin(N\omega/2)}{\sin(\omega/2)} \cdot e^{-j(N-1)\omega/2} \tag{1b} \]

この形から2つの重要な事実が読み取れます。

  1. 振幅特性は \(\left|\frac{1}{N}\frac{\sin(N\omega/2)}{\sin(\omega/2)}\right|\) で、\(\sin(N\omega/2) = 0\) となる \(\omega = 2\pi k/N\) (\(k = 1, 2, \dots\) 、\(N\) の倍数を除く)で厳密にゼロになります。これが特性表でいう「ゼロ点」です。
  2. 位相特性は \(-(N-1)\omega/2\) で、\(\omega\) に対して厳密に線形です。群遅延は位相の負の傾きとして定義されるため、
\[ \tau_g(\omega) = -\frac{d}{d\omega}\left(-\frac{(N-1)\omega}{2}\right) = \frac{N-1}{2} \tag{1c} \]

となり、すべての周波数で群遅延が定数 \((N-1)/2\) に一致します。これはSMAのインパルス応答が左右対称(\(h[n] = h[N-1-n]\) )であることに由来する線形位相FIRフィルタの一般的性質で、後述のWMA・EMAとの本質的な違いになります。

実際に数値検証すると、\(N=15\) のとき理論値 \((N-1)/2 = 7.0\) サンプルちょうどが、DC(\(\omega=0\) )だけでなく全周波数域で成立することが scipy.signal.group_delay で確認できます(詳細は後述の比較実験)。

特性

  • 長所: 実装が簡単。高周波ノイズの除去に効果的。位相が厳密に線形なので波形が歪まない
  • 短所: 急激な変化への応答が遅い(遅延が \(\frac{N-1}{2}\) サンプル)。窓幅分のデータをバッファとして保持する必要がある

Python実装

import numpy as np

def sma(x, window):
    """単純移動平均フィルタ"""
    kernel = np.ones(window) / window
    # 'valid' モードで出力(端の不完全な窓を除外)
    return np.convolve(x, kernel, mode='same')

NumPyのconvolveを使うと、SMAは畳み込み演算として簡潔に書けます。

O(1)実装:再帰的な移動和と数値誤差の蓄積

式(1)を素朴に実装すると1ステップあたり \(O(N)\) の計算量がかかりますが、直前の和から差分だけ更新する再帰式を使えば \(O(1)\) にできます。

\[ S_t = S_{t-1} + x_t - x_{t-N}, \qquad y_t = \frac{S_t}{N} \tag{1d} \]

ここで窓内の和は

\[ S_t = \sum_{i=0}^{N-1} x_{t-i} \]

です。新しい値を足し、窓から外れた値を引くだけなので、組み込みシステムのようにCPU資源が限られる環境でよく使われます。

ただし、この再帰式は浮動小数点誤差が時間とともに蓄積するという落とし穴があります。実際に検証してみます。

import numpy as np

np.random.seed(0)
n = 2_000_000
window = 50
x = np.random.normal(0, 1, n).astype(np.float32)

# 再帰式(float32)で計算した移動平均
y_rec = np.zeros(n, dtype=np.float32)
s = np.float32(np.sum(x[:window], dtype=np.float64))
y_rec[window - 1] = s / window
for t in range(window, n):
    s = s + x[t] - x[t - window]
    y_rec[t] = s / window

# 都度再計算した厳密値(float64のcumsumで比較用に構築)
csum = np.cumsum(x.astype(np.float64))
y_direct = np.empty(n)
y_direct[: window - 1] = np.nan
y_direct[window - 1] = csum[window - 1] / window
y_direct[window:] = (csum[window:] - csum[:-window]) / window

valid = ~np.isnan(y_direct)
err = np.abs(y_rec[valid] - y_direct[valid])
print(f"最大誤差: {err.max():.3e}")
print(f"t=1000 の誤差: {err[1000 - window + 1]:.3e}")
print(f"t=2,000,000 の誤差: {err[-1]:.3e}")

実行結果:

最大誤差: 1.165e-05
t=1000 の誤差: 2.716e-07
t=2,000,000 の誤差: 9.250e-06

float32 で200万サンプルの再帰計算を続けると、都度再計算する厳密値との誤差が最終的に約 \(1.2 \times 10^{-5}\) まで蓄積することが確認できました。センサログを長時間動かし続ける組み込み用途では、一定間隔で窓内の和を直接再計算してドリフトをリセットする、あるいは float64 を使うといった対策が必要です。

加重移動平均 (Weighted Moving Average, WMA)

定義

WMAは、新しいデータほど大きな重みを付ける移動平均です。重みは線形に減少します。

\[ y_t = \frac{\sum_{i=0}^{N-1} (N - i) \cdot x_{t-i}}{\sum_{i=0}^{N-1} (N - i)} = \frac{2}{N(N+1)} \sum_{i=0}^{N-1} (N - i) \cdot x_{t-i} \tag{2} \]

正規化定数 \(\frac{2}{N(N+1)}\) は、重み \(N, N-1, \dots, 1\) の和が等差数列の和の公式

\[ \sum_{i=0}^{N-1} (N-i) = \sum_{k=1}^{N} k = \frac{N(N+1)}{2} \tag{2a} \]

から導かれます。最新のデータ点に重み \(N\) 、1ステップ前のデータ点に重み \(N-1\) 、…、\(N-1\) ステップ前のデータ点に重み \(1\) を付け、この和で正規化します。

なぜ線形位相にならないか

SMAが線形位相(群遅延一定)になる理由は、インパルス応答が左右対称 \(h[n] = h[N-1-n]\) だったことでした。WMAのインパルス応答は \(h[n] \propto (n+1)\) (新しいサンプルほど大きい)であり、\(h[0] \ne h[N-1]\) なので対称ではありません。したがってWMAは線形位相を持たず、群遅延は周波数に依存して変化します

これは特性表で「WMAの遅延はSMAより小さい」とだけ書くと不正確で、正しくは「WMAの群遅延はDC付近では小さいが周波数が上がるにつれて変化し、場合によっては負の値(信号を先取りするように見える値)も取りうる」という点に注意が必要です。実測値は後述の比較実験の図(群遅延の比較)で確認します。

特性

  • 長所: SMAよりも最新データへの追従が速い(DC付近の群遅延が小さい)
  • 短所: SMAと同様に窓幅分のバッファが必要。重みが線形のため、古いデータの影響が急に0になる(窓の端での不連続)。線形位相を持たないため、周波数によって群遅延が変化し波形形状がわずかに歪む

Python実装

def wma(x, window):
    """加重移動平均フィルタ"""
    weights = np.arange(1, window + 1, dtype=float)
    weights /= weights.sum()
    return np.convolve(x, weights[::-1], mode='same')

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

定義

EMAは再帰的に計算される移動平均で、過去のすべてのデータに指数関数的に減衰する重みを与えます。

\[ y_t = \alpha \cdot x_t + (1 - \alpha) \cdot y_{t-1} \tag{3} \]

\(\alpha\) (\(0 < \alpha \le 1\) )は平滑化係数です。\(\alpha\) が小さいほど平滑化が強くなります。

SMAの窓幅 \(N\) と対応させる場合、\(\alpha = \frac{2}{N + 1}\) がよく使われます。

閉形式(無限インパルス応答)の導出

式(3)の漸化式を初期時刻まで展開すると、EMAが実質的に無限長のインパルス応答を持つIIRフィルタであることが分かります。

\[ y_t = \alpha \sum_{i=0}^{t-1} (1-\alpha)^i \, x_{t-i} \; + \; (1-\alpha)^t \, y_0 \tag{3a} \]

第2項は \(t \to \infty\) で \(0\) に収束するため、定常状態では重み \(w_i = \alpha(1-\alpha)^i\) (\(i = 0, 1, 2, \dots\) )を持つ無限インパルス応答フィルタとみなせます。重みの総和が1になることは、等比級数の和の公式から確認できます。

\[ \sum_{i=0}^{\infty} \alpha (1-\alpha)^i = \alpha \cdot \frac{1}{1-(1-\alpha)} = 1 \tag{3b} \]

この重み分布の1次モーメント(重み付き平均インデックス)を計算すると、DC(\(\omega=0\) )付近の等価遅延が得られます。

\[ \tau_{DC} = \sum_{i=0}^{\infty} i \cdot \alpha(1-\alpha)^i = \frac{1-\alpha}{\alpha} \tag{3c} \]

\(N=15\) に対応する \(\alpha = 2/16 = 0.125\) を代入すると \(\tau_{DC} = 7.0\) となり、SMAの \((N-1)/2 = 7.0\) と一致するよう設計されていることが確認できます。ただし式(1c)のSMAの群遅延が全周波数で一定の7.0だったのに対し、EMAの群遅延はDCでのみ7.0で、周波数が上がるにつれて単調に減少し、高周波側では負の値に漸近するという違いがあります(WMAと同様、EMAも線形位相を持ちません)。数値的には scipy.signal.group_delay で計算した結果、\(\omega/\pi = 0\) で7.0、\(\omega/\pi = 0.5\) 付近で約 \(-0.45\) となることを確認しました(詳細は後述の図参照)。

特性

  • 長所: 過去データのバッファが不要(1ステップ前の値だけ保持すればよい)。メモリ効率が良い。新旧データの重みが滑らかに減衰する
  • 短所: パラメータ \(\alpha\) の選択が直感的でない場合がある。線形位相を持たないため群遅延が周波数依存

Python実装

def ema(x, alpha):
    """指数移動平均フィルタ"""
    y = np.zeros_like(x, dtype=float)
    y[0] = x[0]
    for t in range(1, len(x)):
        y[t] = alpha * x[t] + (1 - alpha) * y[t - 1]
    return y

EMAの周波数特性(ゲイン・位相特性)の詳細な解析、および深層学習における「重みEMA」への応用研究動向は 指数移動平均(EMA)フィルタの周波数特性 を参照してください。

エッジケース:外れ値は完全には忘れない

式(3b)からわかるとおり、EMAの重みは理論上 \(i \to \infty\) でしか0になりません。つまり、1つの外れ値の影響は指数的に減衰するものの、いつまでも完全にはゼロになりません。これに対しSMA・WMAは窓 \(N\) サンプルを過ぎれば影響が厳密に0になります。実際に単発のスパイク(\(x_{20}=100\) 、それ以外は0)を入力して確認します。

n = 100
x = np.zeros(n)
x[20] = 100.0
window = 15
alpha = 2 / (window + 1)

y_sma = sma(x, window)
y_ema = ema(x, alpha)

print("SMAの非ゼロ区間:", np.nonzero(y_sma)[0].min(), "-", np.nonzero(y_sma)[0].max())
for dt in [1, 5, 15, 30, 50, 79]:
    idx = 20 + dt
    print(f"t=20+{dt:>2}: EMA残差={y_ema[idx]:.3e}  理論値 alpha*(1-alpha)^dt*100={alpha*(1-alpha)**dt*100:.3e}")

実行結果:

SMAの非ゼロ区間: 13 - 27
t=20+ 1: EMA残差=1.094e+01  理論値 alpha*(1-alpha)^dt*100=1.094e+01
t=20+ 5: EMA残差=6.411e+00  理論値 alpha*(1-alpha)^dt*100=6.411e+00
t=20+15: EMA残差=1.687e+00  理論値 alpha*(1-alpha)^dt*100=1.687e+00
t=20+30: EMA残差=2.276e-01  理論値 alpha*(1-alpha)^dt*100=2.276e-01
t=20+50: EMA残差=1.575e-02  理論値 alpha*(1-alpha)^dt*100=1.575e-02
t=20+79: EMA残差=3.278e-04  理論値 alpha*(1-alpha)^dt*100=3.278e-04

SMA(\(N=15\) )はスパイクの影響が \(t=13\) 〜\(27\) の15サンプルにきっちり収まって消える一方、EMAは79サンプル後(\(t=99\) )でもまだ \(3.3 \times 10^{-4}\) の残差が残っており、理論式 \(\alpha(1-\alpha)^{\Delta t} \times 100\) と完全に一致しました。異常値検知など「過去の影響を確実に断ち切りたい」用途ではSMA/WMAの方が安全です。

エッジケース:NaNの伝播

式(3)の再帰計算は、1つの NaN が入力されるとそれ以降すべての出力がNaNになるという重大な弱点を持ちます。SMAは窓が過ぎれば回復しますが、EMAは回復しません。

import numpy as np

np.random.seed(1)
x_nan = np.random.normal(0, 1, 30)
x_nan[10] = np.nan

y_ema_nan = ema(x_nan, alpha)
y_sma_nan = sma(x_nan, window)

print("EMA (t=9..14):", np.round(y_ema_nan[9:15], 4))
print("t=10以降すべてNaN:", np.all(np.isnan(y_ema_nan[10:])))
print("SMAがNaNになる区間:", np.where(np.isnan(y_sma_nan))[0].min(), "-", np.where(np.isnan(y_sma_nan))[0].max())

実行結果:

EMA (t=9..14): [0.3398    nan    nan    nan    nan    nan]
t=10以降すべてNaN: True
SMAがNaNになる区間: 3 - 17

EMAは t=10 以降のすべての値がNaNになり、リアルタイムシステムでは欠測値を事前に補間するか、NaNを検出したら\(y_{t-1}\) を保持する(更新をスキップする)などの防御的な実装が必須です。SMAは窓(この例では \(N=15\) なのでインデックス3〜17)を過ぎれば自動的に回復します。

3つのフィルタの比較

比較実験

ノイズを含む信号に対して、3つのフィルタを適用して比較します。以下のコードを実際に実行し、得られた図をそのまま掲載しています。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
n = 200
t = np.linspace(0, 4 * np.pi, n)

# 元の信号 + ノイズ
signal = np.sin(t) + 0.5 * np.sin(5 * t)
noise = np.random.normal(0, 0.5, n)
observed = signal + noise

# フィルタ適用
window = 15
alpha = 2 / (window + 1)

y_sma = sma(observed, window)
y_wma = wma(observed, window)
y_ema = ema(observed, alpha)

# プロット
fig, axes = plt.subplots(2, 2, figsize=(12, 8), sharex=True)

axes[0, 0].plot(t, observed, alpha=0.4, label='Observed')
axes[0, 0].plot(t, signal, 'k--', label='True signal')
axes[0, 0].set_title('Original')
axes[0, 0].legend()

axes[0, 1].plot(t, observed, alpha=0.3)
axes[0, 1].plot(t, y_sma, label=f'SMA (N={window})')
axes[0, 1].set_title('SMA')
axes[0, 1].legend()

axes[1, 0].plot(t, observed, alpha=0.3)
axes[1, 0].plot(t, y_wma, label=f'WMA (N={window})')
axes[1, 0].set_title('WMA')
axes[1, 0].legend()

axes[1, 1].plot(t, observed, alpha=0.3)
axes[1, 1].plot(t, y_ema, label=f'EMA (α={alpha:.3f})')
axes[1, 1].set_title('EMA')
axes[1, 1].legend()

for ax in axes.flat:
    ax.grid(True, alpha=0.3)
    ax.set_ylabel('Amplitude')

axes[1, 0].set_xlabel('Time')
axes[1, 1].set_xlabel('Time')
plt.tight_layout()
plt.show()

時間領域での比較(観測信号・真の信号・SMA・WMA・EMA)

定量評価:ノイズ除去性能と追従性

図だけでは違いが分かりにくいので、真の信号(signal)に対する誤差を定量化します。畳み込みの端効果を避けるため、窓の半分(edge = window // 2)を除いた区間で評価します。

edge = window // 2

def rmse(y):
    return np.sqrt(np.mean((y[edge:-edge] - signal[edge:-edge]) ** 2))

print(f"フィルタなし: RMSE={np.sqrt(np.mean((observed[edge:-edge]-signal[edge:-edge])**2)):.4f}")
print(f"SMA:          RMSE={rmse(y_sma):.4f}")
print(f"WMA:          RMSE={rmse(y_wma):.4f}")
print(f"EMA:          RMSE={rmse(y_ema):.4f}")

実行結果:

フィルタなし: RMSE=0.4686
SMA:          RMSE=0.2713
WMA:          RMSE=0.3227
EMA:          RMSE=0.4152

同一窓幅・対応する \(\alpha\) の条件では、SMA(RMSE 0.2713)が最もノイズ除去性能が高く、EMA(RMSE 0.4152)が最も低い結果になりました。これは、この信号に含まれる高周波成分(\(0.5\sin(5t)\) )に対して、EMAの単調ななだらかな減衰特性よりも、SMAのシャープな阻止特性(後述のゼロ点)の方がこの実験条件では有利に働いたためです。ただし、この優劣は信号の周波数成分と窓幅・\(\alpha\) の組み合わせに強く依存する点に注意してください(後述の周波数特性の比較を参照)。

実測遅延:mode='same' の落とし穴

特性表では「SMAの遅延は \((N-1)/2\) サンプル」としばしば書かれますが、この記事のコードのように np.convolve(x, kernel, mode='same') を使うと、実際にはその遅延が観測されません。理由を確認します。

from scipy.signal import lfilter

# 因果的(リアルタイム実装可能)なSMA: 過去と現在のサンプルのみ使用
b = np.ones(window) / window
y_causal = lfilter(b, [1.0], observed)

# 本記事のコードが使う mode='same' のSMA
y_same = np.convolve(observed, np.ones(window) / window, mode='same')

shift = (window - 1) // 2
diff = np.max(np.abs(y_causal[shift:] - y_same[: len(y_causal) - shift]))
print(f"因果的SMAを{shift}サンプル前にずらした系列と mode='same' の最大差: {diff:.2e}")

実行結果:

因果的SMAを7サンプル前にずらした系列と mode='same' の最大差: 0.00e+00

差はちょうど0でした。つまり mode='same' の出力は、因果的な(未来のデータを使わない)SMAの出力を \((N-1)/2\) サンプルだけ時間的に前にずらしたものと完全に一致します。言い換えると、mode='same' は窓の中心を出力位置に合わせるために未来のサンプルを使っており、**非因果的(non-causal)**です。バッチ処理・オフライン解析では問題になりませんが、株価のリアルタイム表示やセンサのオンライン処理のように「その時点までのデータしか使えない」用途で mode='same' をそのまま使うと、実際には未来の情報が混入した結果になってしまいます。リアルタイム用途では scipy.signal.lfilter や式(1d)の再帰式など、因果的な実装を使う必要があります。

これは、先ほどの相関ベースの遅延測定で「SMA/WMAの実測遅延が理論値より小さく見える」現象の直接の原因でもあります。理論値 \((N-1)/2=7\) は因果的フィルタの遅延であり、mode='same' はその遅延を意図的に相殺してしまうため、真の信号との相互相関のピーク位置がほぼラグ0に現れます。

特性比較表

特性SMAWMAEMA
重みの分布均一線形減衰指数減衰
メモリ使用\(O(N)\) (窓幅分のバッファ)\(O(N)\)\(O(1)\) (前回値のみ)
計算量(1ステップ)\(O(N)\) (式(1d)の差分法で\(O(1)\) )\(O(N)\)\(O(1)\)
過去データの影響窓外で完全に0窓外で完全に0指数的に減衰するが理論上0にはならない
群遅延全周波数で一定 \((N-1)/2=7.0\)周波数依存(DC付近 約4.7、高域で負)周波数依存(DC \((1-\alpha)/\alpha=7.0\) 、高域で負)
線形位相あり(対称インパルス応答)なしなし
ステップ応答直線的に上昇曲線的に上昇指数的に収束
外れ値1点の影響\(N\) サンプル後に完全消失\(N\) サンプル後に完全消失指数減衰するが理論上は無限に残る
NaN 1点の伝播窓を過ぎれば回復窓を過ぎれば回復それ以降すべてNaNになる

(群遅延の数値は \(N=15\) 、\(\alpha=2/(N+1)=0.125\) での実測値)

周波数特性の比較

SMAの伝達関数は以下のとおりです。

\[ G_{SMA}(z) = \frac{1}{N} \sum_{i=0}^{N-1} z^{-i} = \frac{1}{N} \cdot \frac{1 - z^{-N}}{1 - z^{-1}} \tag{4} \]

SMAのゲイン特性にはゼロ点(式(1b)より \(\omega/\pi = 2k/N\) で厳密にゲイン0になる点)が存在するのに対し、WMA・EMAのゲインは単調に近い形で減少します。実際に scipy.signal.freqz で計算すると、次のような数値が得られます。

from scipy.signal import freqz

b_sma = np.ones(window) / window
w, h_sma = freqz(b_sma, [1.0], worN=4096)

w_wma = np.arange(1, window + 1, dtype=float)
w_wma /= w_wma.sum()
_, h_wma = freqz(w_wma[::-1], [1.0], worN=4096)

_, h_ema = freqz([alpha], [1, -(1 - alpha)], worN=4096)

# k=1 のゼロ点ちょうどの周波数
om_zero = 2 * np.pi * 1 / window
_, h_sma_zero = freqz(b_sma, [1.0], worN=[om_zero])
_, h_wma_zero = freqz(w_wma[::-1], [1.0], worN=[om_zero])
_, h_ema_zero = freqz([alpha], [1, -(1 - alpha)], worN=[om_zero])

print(f"ゼロ点 omega/pi = {om_zero/np.pi:.4f} でのゲイン:")
print(f"  SMA: {20*np.log10(abs(h_sma_zero[0])):.1f} dB")
print(f"  WMA: {20*np.log10(abs(h_wma_zero[0])):.1f} dB")
print(f"  EMA: {20*np.log10(abs(h_ema_zero[0])):.1f} dB")

実行結果:

ゼロ点 omega/pi = 0.1333 でのゲイン:
  SMA: -319.8 dB
  WMA: -10.4 dB
  EMA: -10.3 dB

SMAは設計どおりこの周波数を完全に(数値誤差レベルまで)遮断できるのに対し、WMA・EMAは同じ周波数で約-10dBしか減衰できません。ノイズの周波数があらかじめ分かっている場合、窓幅 \(N\) を調整してSMAのゼロ点をその周波数に合わせるのが最も効果的であることが定量的に確認できます。

-3dB遮断周波数(正規化周波数、\(\times\pi\) rad/sample)を比較すると、SMA 0.0592、WMA 0.0718、EMA 0.0426となり、同じ窓幅・対応する \(\alpha\) の条件ではEMAが最も低域まで通過帯域が狭い(強く平滑化する)ことも分かります。

周波数応答の比較(SMAのゼロ点 vs WMA/EMAのなだらかな減衰)

前節で導出したとおり、SMAだけがすべての周波数で群遅延が一定(線形位相)であるのに対し、WMA・EMAは周波数によって群遅延が変化します。実際に scipy.signal.group_delay で計算すると、下図のようにSMAの群遅延はDC〜ナイキスト周波数まで完全に水平な直線(7.0で一定)になりますが、WMA・EMAはDC付近で7.0未満の値から出発し、高周波側では負の値に転じることが確認できます。

群遅延の比較(SMAのみ線形位相で群遅延が一定)

各フィルタの周波数特性の詳細な比較(バターワース・チェビシェフとの比較を含む)は ローパスフィルタの設計と比較 も参照してください。

実務上の落とし穴・エッジケースまとめ

これまでの導出・実験で確認した注意点を整理します。

  • 偶数窓幅では出力が半サンプルずれる: SMA・WMAは窓の中心が出力位置になるため、\(N\) が奇数なら整数サンプル位置に中心が来ますが、\(N\) が偶数だと中心が2サンプルの中間(半サンプル)にずれます。ステップ応答で確認すると、\(N=15\) (奇数)では立ち上がり中の値が [0.4, 0.467, 0.533, 0.6, ...] ときれいに \(1/15\) 刻みになりますが、\(N=16\) (偶数)では [0.375, 0.438, 0.5, 0.562, ...] と半サンプル分ずれた位置が中心(0.5)になります。窓幅は奇数を選ぶと解析上シンプルです。
  • mode='same' は非因果的: 前述のとおり、np.convolve(..., mode='same') は \((N-1)/2\) サンプル分未来を先取りしています。リアルタイム処理では lfilter や再帰式を使う。
  • EMAは外れ値を完全には忘れない: 理論上無限に影響が残る(指数減衰)。異常値の影響を確実に断ち切りたい場合はSMA/WMAが安全。
  • EMAはNaN 1つで以降すべて壊れる: 欠測値の事前補完、またはNaN検出時の更新スキップが必須。
  • 再帰的SMA(式(1d))は浮動小数点誤差が蓄積する: 長時間運用するシステムでは定期的な再計算やfloat64の使用を検討する。

最新の研究動向:深層学習の時系列予測におけるSMA/EMA分解

移動平均は古典的な信号処理の道具ですが、2023年以降のTransformer全盛の時系列予測研究でも重要な役割を果たしています。

Zeng, Chen, Zhang, Xu (2023) の “Are Transformers Effective for Time Series Forecasting?"(AAAI 2023, Vol. 37, pp. 11121–11128, arXiv:2205.13504 )が提案した DLinear は、入力系列を「トレンド成分」と「季節成分(残差)」に分解してから、それぞれを単純な線形層で予測するという、驚くほどシンプルな構成で当時の多くのTransformer系SOTAモデルを上回り、大きな議論を呼びました。この分解には本記事のSMAとまったく同じ考え方(AvgPool1d によるパディング付き移動平均)が使われており、トレンド成分と季節成分はそれぞれ次のように定義されます。

\[ X_t = \text{AvgPool}(\text{Padding}(X)), \qquad X_s = X - X_t \]

これに続き、Stitsyuk & Choi (2025) の “xPatch: Dual-Stream Time Series Forecasting with Exponential Seasonal-Trend Decomposition”(AAAI 2025, pp. 20601–20609, arXiv:2412.17323 )は、DLinearが使うSMAベースの分解をEMAベースの分解に置き換えることでさらに精度を改善できることを報告しています。論文は「SMAは単純すぎるトレンド信号しか作れず、季節成分が複雑になりすぎる」と指摘し、直近のデータをより重視するEMAに置き換えることで、DLinearに対して統一設定でMSE 2.46%・MAE 2.34%、ハイパーパラメータ探索設定ではMSE 5.29%・MAE 3.81%の改善を報告しています。

これはまさに本記事で導出した「SMAはゼロ点による特定周波数の完全遮断に優れる一方、EMAは滑らかで直近重視の応答を持つ」という周波数特性上のトレードオフが、現代の深層学習アーキテクチャの設計判断としてそのまま現れている興味深い例です。なお、EMAを深層学習モデルの重み(勾配降下で更新されるパラメータ)に対して適用し学習を安定化させる「重みEMA」という別の応用については、 指数移動平均(EMA)フィルタの周波数特性 で扱っている研究動向(Morales-Brotons et al., 2024)を参照してください。

用途別の選択ガイド

用途推奨フィルタ理由
センサデータの前処理EMAメモリ効率が良く、リアルタイム処理に適している
株価・金融データのトレンド分析SMA / EMASMAはテクニカル分析の標準指標、EMAは短期トレンドの追従に優れる
組み込みシステムEMA\(O(1)\) のメモリと計算量で実装可能
オフライン信号処理SMA / WMAバッチ処理なのでメモリ制約が緩く、窓幅ベースの制御が直感的
ノイズの周波数が既知の場合SMA窓幅を調整してゼロ点を目的周波数に合わせられる(実測で-320dB遮断を確認)
異常値の影響を確実に断ち切りたい場合SMA / WMAEMAは外れ値の影響が理論上無限に残るため不向き

おすすめ書籍

はじめて学ぶディジタル・フィルタと高速フーリエ変換(三上直樹、CQ出版)

ディジタルフィルタと FFT の原理を基礎から丁寧に解説した定番の入門書です。本記事で扱った処理の理論的背景を体系的に学べます。

※ 上記は Amazon アソシエイトのリンクです。

まとめ

  • SMA は直近 \(N\) サンプルを均一な重みで平均する最も単純なFIRフィルタで、インパルス応答が対称なため線形位相(群遅延が全周波数で一定 \((N-1)/2\) )を持つ。ゼロ点で特定周波数を完全に遮断できる(実測で-320dB)
  • WMA は線形に減衰する重みでSMAよりDC付近の追従が速いが、線形位相を持たず群遅延は周波数依存
  • EMA は \(O(1)\) メモリの再帰フィルタで組み込み・リアルタイム向きだが、外れ値の影響が理論上無限に残り、NaN1つで以降すべての出力が壊れるという弱点がある
  • np.convolve(..., mode='same') は非因果的(未来のサンプルを使用)であり、理論上の遅延 \((N-1)/2\) を相殺してしまうため、リアルタイム実装では lfilter や再帰式を使う必要がある
  • 深層学習の時系列予測でも、SMAベースの分解(DLinear, 2023)からEMAベースの分解(xPatch, 2025)への進化という形で、本記事のSMA/EMAのトレードオフがそのまま現れている

関連記事

参考文献

  • Smith, S. W. (1997). The Scientist and Engineer’s Guide to Digital Signal Processing. California Technical Publishing.
  • Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
  • Zeng, A., Chen, M., Zhang, L., & Xu, Q. (2023). Are Transformers Effective for Time Series Forecasting? Proceedings of the AAAI Conference on Artificial Intelligence, 37, 11121–11128. arXiv:2205.13504
  • Stitsyuk, A., & Choi, J. (2025). xPatch: Dual-Stream Time Series Forecasting with Exponential Seasonal-Trend Decomposition. Proceedings of the AAAI Conference on Artificial Intelligence, 20601–20609. arXiv:2412.17323