Wizard Notes

音楽信号解析の技術録、音楽のレビューおよび分析、作曲活動に関する雑記です

SciPyを使ったデジタルフィルタの周波数応答プロット

この記事をシェアする

実装したIIR/FIRフィルタが正しく動作しているかどうか検証する必要があったため、

デジタルフィルタの周波数応答を見るスクリプトを試してみました。

基本的には scipy.signal.freqz() の example とほぼ同じで、

プロット部分を少し整理したスクリプトとなっています。

使用例

8次FIR (窓関数法) フィルタによるローパスフィルタ

f:id:Kurene:20200718103002p:plain

6次IIR (Butterworth) フィルタによるバンドパスフィルタ

f:id:Kurene:20200718102919p:plain

1次IIRフィルタによるオールパスフィルタ

f:id:Kurene:20200718103314p:plain

実装

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.unwrap(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()