今回は、前回までのローパスフィルタ、ハイパスフィルタ、バンドパス/ストップフィルタの周波数特性を表示する方法について紹介する。
1. フィルタ周波数特性
前回までに設計したフィルタの各周波数に対する利得、位相ずれの周波数特性をPython, Scipyの関数を使って表示することができる。
基本的にはフィルタ設計で得られたフィルタ伝達関数の係数b,aを用いてw, h = signal.freqs(b, a)とすると角周波数(2π f):wと周波数応答:hを返す。
周波数応答hは複素数になっており、絶対値をとると(abs関数用いて)利得、角度をとると(angle関数用いて)位相ずれを算出することができる。
設計時との違いはフィルタ伝達関数の係数b,aを算出する時にsignal.butter()関数の3番目のanalogフラグをTrueに設定し、通過域端周波数と阻止域端周波数はサンプリング周波数を用いた正規化周波数ではなく、角周波数で設定する。(ここではfreqsの出力も角周波数のため、通常の周波数のように扱っても辻褄があう）
scipy.signal.buttord — SciPy v1.7.0 Manual
scipy.signal.freqs — SciPy v1.7.0 Manual
2. ローパスフィルタの周波数特性
以下コードを用いてローパスフィルタの周波数特性を算出する。freqs()はなにも指定しなければ、200点の周波数データを出力するが、np.logspace()とすることで対数周波数軸で特性を出力する。
np.logspace(1, 3, 100)は対数軸上の101から103で100点出力する。
import numpy as np from scipy import signal from matplotlib import pyplot as plt fp = 300 # 通過域端周波数[Hz] fs = 500 # 阻止域端周波数[Hz] gpass = 3 # 通過域端最大損失[dB] gstop = 40 # 阻止域端最小損失[dB] N, Wn = signal.buttord(fp, fs, gpass, gstop,True) #オーダーとバターワースの正規化周波数を計算 b, a = signal.butter(N, Wn, "low",True) #フィルタ伝達関数の分子と分母を計算 w, h = signal.freqs(b, a, np.logspace(1, 3, 100))
利得を表示するにはhの絶対値をとって20*log10によるデシベルに変換する。位相ずれはhのangleと位相飛びをつなげるためにunwrap関数を利用し、最後に180/np.piをかけてラジアンから角度に変換している。
fig, ax1 = plt.subplots() ax1.set_title('Digital filter frequency response') ax1.plot(w, 20 * np.log10(abs(h)), 'b') ax1.set_ylabel('Amplitude [dB]', color='b') ax1.set_xlabel('Frequency [rad/sample]') ax2 = ax1.twinx() angles = np.unwrap(np.angle(h)) ax2.plot(w, angles*180/np.pi, 'g') ax2.set_ylabel('Angle (Degrees)', color='g') ax2.grid()
実行すると以下の300Hzのカットオフ、600Hzで-40dBで設計したバタワースローパスフィルタの周波数特性が得られる。
3. ハイパス、バンドパス/ストップフィルタの周波数特性
ローパスと同じコードで通過域、阻止域とbutter()関数でhigh,band,bandstopとしてハイパス、バンドパス、バンドストップフィルタの周波数特性をそれぞれ計算する。
ハイパスフィルタ
fp =500 # 通過域端周波数[Hz] fs = 300 # 阻止域端周波数[Hz]
バンドパスフィルタ
fp =np.array([1000, 3000] #通過域端周波数[Hz] fs =np.array([500, 6000]) #阻止域端周波数[Hz]
バンドストップフィルタ
fp = np.array([500, 6000]) #通過域端周波数[Hz] fs = np.array([1000, 3000]) #阻止域端周波数[Hz]
3. まとめ
今回はそれぞれのフィルタの周波数特性を計算・表示する方法を紹介した。
