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):
= alt.Data(values=sum([[{
source 'Real Part' : pt.real,
'Imaginary Part': pt.imag,
"n": i
for pt in zeros(i)] for i in range(1, n + 1)], []))
}
= alt.binding_range(min=1, max=n, step=1,name="n", )
slider = alt.selection_single(fields=['n'], bind=slider, init={'n': 1})
select_order 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()
32) plot_zeros(
はじめに
最も単純なローパスフィルタとして,移動平均フィルタが広く知られています. しかしながら,その詳細な性質を確認する機会はあまり無いのではと思います. そこで,移動平均フィルタの伝達関数,複素平面における零点の分布,周波数特性を示します. そして,ローパスフィルタとしての特性を再確認します.
やったこと
- 移動平均フィルタについて複素平面上で零点の分布をプロット
- 移動平均フィルタについて周波数特性をプロット
- 移動平均フィルタがローパスフィルタであり,直線位相であること確認
移動平均フィルタについて
移動平均フィルタとは,
複素平面上で零点を求める
前述の伝達関数
まず初めに問題を定義します.
ここで,両辺を
この時,
コンパニオン行列
と置くと,
以下に零点のプロットを示します.この結果より,以下を確認できます. * 零点は単位円上にのみ存在する *
周波数特性を求める
同様に周波数特性を求めます.周波数特性は
import math
def response(n):
= np.linspace(0, 1.0, 257, endpoint=True)
freq = np.sum(np.vander(np.exp(freq * 2 * np.pi * 1j), n), axis=1) / n
resp return resp, freq
def plot_freq_feature(n):
= alt.Data(values = sum([[{
source '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)], []))
}
= alt.Chart(source).encode(
base 'Normalized frequency:Q', axis=alt.Axis(title=None))
alt.X(
)
= base.mark_line(stroke='#57A44C', interpolate='monotone').encode(
amp 'gain:Q', scale=alt.Scale(domain=[-50, 0], clamp=True), axis=alt.Axis(title='Gain(dB)', titleColor='#57A44C')),
alt.Y(
)
= base.mark_line(stroke='#5276A7', interpolate='monotone').encode(
phase 'phase:Q', scale=alt.Scale(domain=[-math.pi, math.pi]), axis=alt.Axis(title='Phase(rad)', titleColor='#5276A7')),
alt.Y(
)
= alt.binding_range(min=1, max=n, step=1,name="n", )
slider = alt.selection_single(fields=['n'], bind=slider, init={'n': 1})
select_order return alt.layer(amp, phase).resolve_scale(
= 'independent'
y
).add_selection(select_order).transform_filter(select_order).interactive()
32) plot_freq_feature(
上図より以下のことが確認できます. - 低周波を通し,高周波をカットするローパスフィスタになっている - 複素平面上で零点の存在する周波数は利得が小さい(カットされている) - 直線位相となっている