移動平均フィルタの零点と周波数特性

信号処理

Jupyter

Published

October 10, 2020

はじめに

最も単純なローパスフィルタとして，移動平均フィルタが広く知られています．しかしながら，その詳細な性質を確認する機会はあまり無いのではと思います．そこで，移動平均フィルタの伝達関数，複素平面における零点の分布，周波数特性を示します．そして，ローパスフィルタとしての特性を再確認します．

やったこと

移動平均フィルタについて複素平面上で零点の分布をプロット
移動平均フィルタについて周波数特性をプロット
移動平均フィルタがローパスフィルタであり，直線位相であること確認

移動平均フィルタについて

移動平均フィルタとは， \(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)\)の固有値に等しいことが解ります．

以下に零点のプロットを示します．この結果より，以下を確認できます． * 零点は単位円上にのみ存在する * \(z=-1\)以外の零点は全て複素共役になっている

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)

上図より以下のことが確認できます． - 低周波を通し，高周波をカットするローパスフィスタになっている - 複素平面上で零点の存在する周波数は利得が小さい（カットされている） - 直線位相となっている

参考

[1]

Eigenvalue-polynomials. http://web.mit.edu/18.06/www/Spring17/Eigenvalue-Polynomials.pdf.