Wizard Notes

Python, JavaScript を使った音楽信号分析の技術録、作曲活動に関する雑記

移動平均と指数平滑移動平均の周波数特性(振幅,位相,群遅延)

はじめに

時系列信号の高周波成分を除去し低周波成分を抽出したり、時系列のトレンドをつかんだりするために、「移動平均」と「指数平滑移動平均」はよく使われる手法です。

音信号処理だけでなく、身近なところでは株価チャートでよく利用されます。

これらはローパスフィルタの一種ですが、具体的にはどのような特性と周波数応答を持つのかを解説したいと思います。

定義

移動平均(Moving Average)とは

移動平均は、一定期間(ウィンドウサイズ)内のデータポイントの平均値を計算し、その結果を平滑化したデータとして用いる手法です。例えば、株価の分析で用いる「5日移動平均」は、直近5日間の終値の平均を計算して、その平均値を毎日更新していくものです。

計算方法

移動平均の計算はシンプルで、次のようなFIRフィルタになります:

ここで、Mはウィンドウサイズ(移動平均を取る期間の長さ)です。この方法では、最新のデータポイントに対して過去のM個のデータの平均を計算します。

特徴と用途

移動平均の最大の特徴は、そのシンプルさと直感的な理解のしやすさです。これは、短期的なノイズを除去し、データの長期的なトレンドを強調するためによく使われます。しかし、ウィンドウサイズが固定されているため、直近のデータに反応しにくく、急な変化には追従しづらいという欠点があります。

指数平滑移動平均(Exponential Moving Average)とは

指数平滑移動平均(Exponential Moving Average, EMA)は、移動平均を改良したもので、直近のデータにより多くの重みを与えることで、データの変動により敏感に反応するように設計されています。これにより、最新の情報をより反映した平滑化が可能です。

計算方法

指数平滑移動平均はフィードバック構造をもつため、IIRフィルタとなります:

α は平滑化係数(0 <α ≤ 1)です。このαが小さくなるほど過去の値の影響が大きくなるとともに、現在の入力 x[n] の影響が大きくなります。すなわち、より高周波の信号を除去するように動作します。

特徴と用途

EMAの特徴は、新しいデータに対して素早く反応することです。直近のデータに重みを多く与えるため、急な変動にも敏感に反応できます。そのため、株価の短期的なトレンドの検出や、リアルタイムでのシステム制御など、最新のデータが特に重要な場面でよく使われます。

データの「滑らかさ」がαの値に依存するため、時系列信号の時間解像度(サンプリング周波数)を考慮し、適切な出力が得られるように α を選ぶ必要があります。

移動平均指数平滑移動平均の使い分け

移動平均指数平滑移動平均のどちらを使うべきかは、分析の目的によって異なります。もしデータのトレンドを長期的に見たい場合や、短期的なノイズを除去したい場合は、シンプルな移動平均が適しています。一方、データの急な変化をキャッチしたい場合や、最新の情報に迅速に反応する必要がある場合は、指数平滑移動平均の方が適しています。

なお、リアルタイムシステムとして考えると、単純な計算量やメモリの観点では、移動平均より指数平滑移動平均のほうが若干優れていると言えます。これは特にウィンドウサイズが大きい、すなわち急峻なローパスフィルタを作成しようとすると差が顕著になります。

周波数特性

サンプリング周波数を48 kHzとしてPythonで振幅特性,位相特性,群遅延を計算しました。

移動平均 (MA)

MAフィルタは時間領域で矩形窓なので、周波数領域ではsinc関数 (Lanczos window) のような形状と、それに基づいた特性になります。

線形位相であるため、群遅延は一定で (M-1)/2 となります。

ローパスフィルタとして使う場合には注意が必要で、高周波の帯域によっては十分に除去できない帯域がでてきます。また、プロットの通り、ウィンドウサイズMに応じて櫛形フィルタのような特性になります。したがって、零点となる周波数や、フィルタのサイドローブと除去したい成分の帯域の重なりを意識した利用が望ましいかと思います。

# M が奇数のとき
>python test_ma.py
Calculated Group Delay at w=0 for M=1: -0.00 samples
Calculated Group Delay at w=0 for M=3: 1.00 samples
Calculated Group Delay at w=0 for M=5: 2.00 samples
Calculated Group Delay at w=0 for M=15: 7.00 samples

# M が偶数のとき
>python test_ma.py
Calculated Group Delay at w=0 for M=2: 0.50 samples
Calculated Group Delay at w=0 for M=4: 1.50 samples
Calculated Group Delay at w=0 for M=8: 3.50 samples
Calculated Group Delay at w=0 for M=16: 7.50 samples

指数平滑移動平均 (EMA)

移動平均とは異なり、一般的なローパスフィルタ的な特性を持ちます。 非線形位相なフィルタであることに注意が必要です。

Pythonコード

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

def plot(b, a, label, fs=48000, worN=2048, amp_ylim=[-40,0]):
    # 周波数応答と群遅延を計算
    w, h = signal.freqz(b, a, worN=worN)  # 正規化周波数 (ラジアン/サンプル)
    angles = np.unwrap(np.angle(h))  # 位相のアンラッピング
    gd = -np.diff(angles) / np.diff(w)  # 群遅延の計算 (サンプル単位)
    
    # 振幅応答のプロット
    amp_dB = 10 * np.log10(np.abs(h+1e-19)**2)
    
    plt.subplot(3, 1, 1)
    plt.plot(w / np.pi * (fs / 2), amp_dB, label=label)  # Hzでプロット
    plt.ylabel('Power [dB]', color='b')
    plt.xlabel('Frequency [Hz]')
    plt.grid()
    plt.ylim(amp_ylim)
    plt.legend()

    # 位相応答のアンラッピングとプロット
    plt.subplot(3, 1, 2)
    plt.plot(w / np.pi * (fs / 2), angles, label=label)  # Hzでプロット
    plt.ylabel('Unwrapped Angle [rad]', color='g')
    plt.grid()
    plt.axis('tight')
    plt.xlabel('Frequency [Hz]')
    plt.legend()
    
    # 群遅延のプロット
    plt.subplot(3, 1, 3)
    plt.plot(w[:-1] / np.pi * (fs / 2), gd, label=label)  # 周波数をHzでプロット
    plt.legend()
    plt.ylabel('Group Delay [samples]')
    plt.xlabel('Frequency [Hz]')
    plt.grid()
    plt.ylim([0, 10])
    
    return w, amp_dB, angles, gd
    
plt.figure(figsize=(8,8))
M_list = [2, 4, 8, 16]
#M_list = [1, 3, 5, 15]
for M in M_list:
    b = np.ones(M) / M 
    a = [1]
    w, amp_dB, angles, gd = plot(b, a, label=f"M={M}", fs=48000)
    print(f'Calculated Group Delay at w=0 for M={len(b)}: {gd[0]:.2f} samples')

plt.tight_layout()
plt.show()
plt.close()

plt.figure(figsize=(8,8))
alpha_list = np.arange(0.9, 0.0, -0.15)**2
for alpha in alpha_list:
    b = 1.0*alpha
    a = np.array([1.0, -(1.0-alpha)])
    w, amp_dB, angles, gd = plot(b, a, label=f"alpha={alpha:1.2f}", fs=48000, amp_ylim=[-20,0])

plt.tight_layout()
plt.show()
plt.close()

まとめ

移動平均指数平滑移動平均は、それぞれの特徴を理解した上で、目的に応じて使い分けることが重要です。データのトレンドを捉えたり、ノイズを除去したりするためのローパスフィルタとしての特性として、どちらの手法が適しているかをしっかり考えた上で選択する必要があります。

関連

www.wizard-notes.com