実装したIIR/FIRフィルタが正しく動作しているかどうか検証する必要があったため、
デジタルフィルタの周波数応答（周波数特性）をプロットするスクリプトを作りました。
基本的には scipy.signal.freqz() の example とほぼ同じで、
プロット部分を少し整理したスクリプトとなっています。
使用例
8次FIR (窓関数法) フィルタによるローパスフィルタ
6次IIR (Butterworth) フィルタによるバンドパスフィルタ
1次IIRフィルタによるオールパスフィルタ
実装
import numpy as np import matplotlib.pyplot as plt from scipy import signal ''' Reference: https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.signal.freqz.html ''' # 振幅プロットの最小値 ylim_min_dB = -120 # デジタルフィルタの係数を算出 ## FIR-lowpass by Window function method #b = signal.firwin(80, 0.5, window=('kaiser', 8)) #a = 1.0 ## IIR-Bandpass by Butterworth filter #b, a = signal.butter(6, [np.pi/16, np.pi/8], btype='band') ## All pass filter lam = 0.5 b, a = np.array([lam, -1]), np.array([1, -lam]) # 周波数応答を算出 w, h = signal.freqz(b, a) # 振幅をデシベル(dB)に変換 amp_dB = 20 * np.log10(np.abs(h)/np.max(np.abs(h))) angles = np.angle(h) # プロット fig = plt.figure() plt.title('Digital filter frequency response') # 振幅応答のプロット ax1 = fig.add_subplot(111) plt.plot(w, amp_dB, 'b') plt.ylabel('Amplitude [dB]', color='b') plt.xlabel('Frequency [rad/sample]') ax1.set_ylim(np.minimum(ylim_min_dB, np.min(amp_dB)), 0) # 位相応答のプロット ax2 = ax1.twinx() plt.plot(w, angles, 'g') plt.ylabel('Angle (radians)', color='g') plt.grid() plt.axis('tight') plt.tight_layout() plt.show()
