移動平均とローパスフィルタ

デジタル制御ではセンサのAD変換値を平滑化するために移動平均を使います.この移動平均を連続系で考えてみるために考察します.

サンプリング周波数$ f_s $[Hz], 移動平均数 N個の移動平均は、ローパスフィルタと同じような役割をします.

カットオフ周波数$ f_c $ は、おおよそ下記で表されます.

\[ f_c = \frac{0.443}{\sqrt{N^2 - 1}} \cdot f_s \]

ステップ応答の比較

まずステップ応答を比べてみます.

計算条件:カットオフ周波数44.3Hzで合わせたフィルタ


① サンプリング周波数 1kHz, 10個の移動平均
② サンプリング周波数 2kHz, 20個の移動平均
③ カットオフ周波数 44.3Hzの1次ローパスフィルタ $ \frac{2 \pi \cdot 44.3}{s+(2  \pi \cdot 44.3)} $
④ カットオフ周波数 44.3Hz、減衰比0.7の2次ローパスフィルタ $ \frac{(2 \pi \cdot 44.3)^2}{s^2+2 \cdot 0.7 \cdot (2 \pi \cdot 44.3)s+(2 \pi \cdot 44.3)^2} $

応答としては2次のローパスフィルタに近いと考えたらよいのでしょうか?

ボード線図の比較

次はボード線図で比べてみます.Scilabを使って①、②、③、④を比較した図です.
(考え方は[2]を参考にしました)

ボード線図の見どころは

  1. 先にあげた公式はゲインが-3.0dBになる周波数を求めている
    (4つの曲線はすべて44.3Hz, -3.0dBを通る)
  2. 移動平均のフィルタとしての特性は、移動平均数 ÷ サンプリング周波数で決まる.
  3. 0~100Hzの範囲での移動平均は1次のローパスフィルタより強く、2次のローパスフィルタに似ている
    (先のステップ応答で見られた通り)
  4. 100~500Hzまで見てみると、ところどころ減衰が悪い周波数があり、その周波数では1次LPF程度のレベルで考えたほうが良い

というところだと思います.

100 Hz前後までのボード線図

もう少しひいて1kHz前後までのボード線図


係数の理屈

パターン①

参考文献[1]や[2]を参照して、
サンプリング周波数 $ f_s $で無次元化した周波数 $ \omega = 2 \pi \frac{f}{f_s} $ に対して、 移動平均の伝達関数は

\[ H(\omega) = \frac{1}{N} \frac{e^{-i \omega N/2}}{e^{-i \omega/2}} \frac{\sin\left(\frac{\omega N}{2}\right)}{\sin\left(\frac{\omega}{2}\right)} \]

これからゲインを計算すると

\[ |H(\omega)| = \frac{1}{N} \left|\frac{\sin\left(\frac{\omega N}{2}\right)}{\sin\left(\frac{\omega}{2}\right)}\right| \]
です。これを$ \omega = 0 $ 周りでテイラー展開して(1/N*Sin(x∗N/2)/Sin(x/2) expands at x=0 - Wolfram|Alpha)
\[ |H(\omega)| \approx 1 + \frac{1}{24} (1 - N^2) \omega^2 + O(\omega^4) \]

するとカットオフ周波数 $ f_c $は$ \omega_c = 2 \pi \frac{f_c}{f_s} $ とおいて

\[ \begin{array}{rrl} & \frac{1}{\sqrt{2}} & = 1 + \frac{1}{24} (1 - N^2) \omega^2 \\ \Leftrightarrow & \frac{1}{\sqrt{2}} &= 1 +\frac{1}{24} (1 - N^2) \left(2 \pi \frac{f_c}{f_s}\right)^2 \\ \Leftrightarrow & f_c^2 &= \frac{24 (\sqrt{2} - 1)}{4 \sqrt{2} \pi^2 (N^2 - 1)} f_s^2 \\ \Leftrightarrow & f_c & \approx \frac{0.4219 \dots}{ \sqrt{N^2 - 1}} f_s \end{array} \]
ここで得た0.4219..に対して、テーラー展開の展開点 $ \omega = 0 $ から離れてしまうことによる誤差を勘案して調節した秘伝の係数が 0.44294.. ということらしい。 よく見つけるわ…。

パターン②

ある定義域 $ 0 < \omega \leq \pi/2N $ で

\[ \begin{array}{rl} |H(\omega)|&= \frac{1}{N} \left|\frac{\sin\left(\frac{\omega N}{2}\right)}{\sin\left(\frac{\omega}{2}\right)}\right| \\ &\approx \frac{1}{N \omega}{\sqrt{2(1-\cos N \omega)}} \end{array} \]

近似状況の例: plot 1/N*Sin(N*x/2)/Sin(x/2) and 1/(x*N)*sqrt(2-2*cos(x*N)) at N = 4 between x = 0 and 6 - Wolfram|Alpha

したがって、

\[ \begin{array}{rrl} & \frac{1}{\sqrt{2}} &= \frac{1}{N \omega}{\sqrt{2(1-\cos N \omega)}} \\ \Rightarrow & N \omega &\approx 2.78311 \dots \\ \Leftrightarrow & f_c &\approx \frac{2.78311 \dots / 2 \pi}{N} f_s \\ &&= \frac{0.42946 \dots}{N} f_s \end{array} \]

参考文献

[2] デジタルフィルタを計算する
[3] filter design - 3dB-Cut off frequency of moving average - Signal Processing Stack Exchange
より良い近似を求めて論争が繰り広げられていますが、0.4429..意外に良いんですよね…。

No comments:

Post a Comment