"Analysis of the feature of moving average filter"
最も単純なローパスフィルタとして,移動平均フィルタが広く知られています. しかしながら,その詳細な性質を確認する機会はあまり無いのではと思います. そこで,移動平均フィルタの伝達関数,複素平面における零点の分布,周波数特性を示します. そして,ローパスフィルタとしての特性を再確認します.
移動平均フィルタとは, $n$ サンプル目の入力を$x_{n}$,出力を$y_{n}$,フィルタ長を$N$とすると以下のように表すことができます. $$y_{n} = \frac{1}{N}\sum_{m = 0}^{N} x_{n + m}$$ これは,連続する$N$サンプルの平均を出力とすることを意味しています. 今回は,零点と周波数特性の関係を調べたいので,上式にZ変換を行い周波数領域に持っていきます.
$$ \begin{aligned} \mathcal{Z}(x_{n+m}) & = z^{-m}X(z) \\\\ \mathcal{Z}(y_{n}) & = Y(z) \\\\ & = H(z)X(z) \\\\ & = \frac{1}{N} \left( \sum_{m = 0}^{N} z^{-m} \right) X(z) \\\\ \end{aligned} $$前述の伝達関数$H(z)$の根を求めることで,零点が複素平面上でどのように分布しているかを知ることができます.ここでは,$H(z)$の求根にはコンパニオン行列の固有値を計算する方法を用います. 求根方法の詳細についていかに述べます.
まず初めに問題を定義します.$H(z)$の根を求めることとは,以下の式を満たす$z$を求めることになります.
$$ \begin{aligned} H(z) & = \frac{1}{N} \sum_{m = 0}^{N} z^{-m}\\\\ &= 0 \end{aligned} $$ここで,両辺を$N$倍して以下の式を得ます $$ \begin{aligned} H'(z) & = \sum_{m = 0}^{N} z^{-m} \\\\ &= 0 \end{aligned} $$
この時,$H'(z)$に対するコンパニオン行列$C(H'(Z))$は以下のように求まります.
$$ \begin{aligned} C\left(H'(z)\right) & = \begin{pmatrix} 0 & 0 & \cdots & 0 & -1 \\ 1 & 0 & \cdots & 0 & -1 \\ 0 & 1 & \cdots & 0 & -1 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & -1 \\ \end{pmatrix} \end{aligned} $$コンパニオン行列$C\left(H'(z)\right)$には,以下の様な面白い性質があります.
$$ \begin{aligned} \left| zI = C\left(H'(z)\right) \right| &= H'(z) \end{aligned} $$ここで,左辺は$H'(z)$の固有多項式を表します.従って,
$$ \begin{aligned} \left| zI = C\left(H'(z)\right) \right| &= H'(z) \\ &= 0 \end{aligned} $$と置くと,$z$は$C\left(H'(z)\right)$の固有値に等しいことが解ります.
以下に零点のプロットを示します.この結果より,以下を確認できます.
import numpy as np
import altair as alt
def companion(n: int):
return np.hstack([np.vstack([np.zeros(n-1), np.eye(n-1,n-1)] ), -np.ones(n).reshape(-1,1)])
def zeros(n):
return np.linalg.eig(companion(n))[0]
def plot_zeros(n):
source = alt.Data(values=sum([[{
'Real Part' : pt.real,
'Imaginary Part': pt.imag,
"n": i
} for pt in zeros(i)] for i in range(1, n + 1)], []))
slider = alt.binding_range(min=1, max=n, step=1,name="n", )
select_order = alt.selection_single(fields=['n'], bind=slider, init={'n': 1})
return alt.Chart(source).mark_circle().encode(x=alt.X("Real Part:Q", scale=alt.Scale(domain=[-1.2, 1.2])), y=alt.Y("Imaginary Part:Q", scale=alt.Scale(domain=[-1.2, 1.2]))).add_selection(select_order).transform_filter(select_order).interactive()
plot_zeros(32)
同様に周波数特性を求めます.周波数特性は$H(z)$に対して$z \rightarrow \exp^{-i\omega}$と置換して$\left|H(z)\right|$,$\angle H(z)$を計算することで得られます.
import math
def response(n):
freq = np.linspace(0, 1.0, 257, endpoint=True)
resp = np.sum(np.vander(np.exp(freq * 2 * np.pi * 1j), n), axis=1) / n
return resp, freq
def plot_freq_feature(n):
source = alt.Data(values = sum([[{
'Normalized frequency': f,
'gain': 20 * math.log10(abs(r)),
'phase': math.atan2(r.imag, r.real),
'n': i
} for r, f in zip(*response(i))] for i in range(1, n+1)], []))
base = alt.Chart(source).encode(
alt.X('Normalized frequency:Q', axis=alt.Axis(title=None))
)
amp = base.mark_line(stroke='#57A44C', interpolate='monotone').encode(
alt.Y('gain:Q', scale=alt.Scale(domain=[-50, 0], clamp=True), axis=alt.Axis(title='Gain(dB)', titleColor='#57A44C')),
)
phase = base.mark_line(stroke='#5276A7', interpolate='monotone').encode(
alt.Y('phase:Q', scale=alt.Scale(domain=[-math.pi, math.pi]), axis=alt.Axis(title='Phase(rad)', titleColor='#5276A7')),
)
slider = alt.binding_range(min=1, max=n, step=1,name="n", )
select_order = alt.selection_single(fields=['n'], bind=slider, init={'n': 1})
return alt.layer(amp, phase).resolve_scale(
y = 'independent'
).add_selection(select_order).transform_filter(select_order).interactive()
plot_freq_feature(32)
上図より以下のことが確認できます.