#概要
今回は，ディジタルIIRフィルタの設計をiirfilter()を使ってやってみようと思います．
python×信号処理の参考文献はあまり見つけることができなかったので，困ったときになにか参考になれば幸いです．
(そもそもPythonで信号処理はあまりやらないかも…？)
[【信号処理】Pythonでディジタルフィルタの設計をやってみた！(２)]
(https://qiita.com/trami/items/521f2fc562d8bbc1fba8)
#動作環境
- Windows10(64bit)
- Python 3.7.2
ディジタルフィルタを作成する関数
Pythonでディジタルフィルタを作成するには色々な方法がありそうです．
以下の公式ドキュメントに信号処理で扱う関数がまとめられています．
https://docs.scipy.org/doc/scipy/reference/signal.html
今回はフィルタの設計なので，Filter designとMatlab-style IIR filter designとMatlabのところになります．
関数がたくさんありすぎて，どれを使えばよいのか悩みますが，IIRフィルタの設計であれば，iirdesignとiirfilterの2つかなと思います．
FIRフィルタであれば，firwin()やfirwin2()，kaiserord()などでしょうか？今回は扱いませんが…
-iirdesign(wp, ws, gpass, gstop, ftype, fs)
通過域遮断周波数，阻止域遮断周波数，通過域最大損失量，阻止域最小減衰量をもとにフィルタを設計したいときに使う．
-iirfilter(N, Wn, rp, rs, btype, ftype, fs)
フィルタの次数と遮断周波数をもとにフィルタを設計したいときに使う．
という感じで使い分けます．
フィルタ仕様
今回は設計したいフィルタの仕様は以下とします．
| サンプリング周波数 | $F_S = 20kHz$ |
| 通過域遮断周波数 | $F_p = 3.4kHz$ |
| 阻止域遮断周波数 | $F_s = 4.6kHz$ |
| 通過許容減衰量 | $A_p = 0.4dB$ |
| 阻止域最小減衰量 | $A_s = 30dB$ |
#iirfilterによる実装
では実装していきます．
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
# パラメータ設定
fs = 20 * 10**3 # サンプリング周波数[Hz]
fpass = 3.4 * 10**3 # 通過遮断周波数[Hz]
fstop = 4.6 * 10**3 # 阻止域遮断周波数[Hz]
gpass = 0.4 # 通過域最大損失量[dB]
gstop = 30 # 阻止域最小減衰量[dB]
# 正規化
fn = fs/2 # ナイキスト周波数
wp = fpass/fn # 正規化通過遮断周波数(無次元)
ws = fstop/fn # 正規化阻止域遮断周波数(無次元)
l1, l2, l3, l4 = 'butter', 'cheby1', 'cheby2', 'ellip'
fig, ax = plt.subplots()
# ButterWorth filter
n1, wn1 = signal.buttord(wp, ws, gpass, gstop, fs=fs)
b1, a1 = signal.iirfilter(n1, wn1, rp=gpass, rs=gstop, btype='low', ftype=l1)
f1, h1 = signal.freqz(b1, a1, fs=fs)
ax.semilogx(f1, 20*np.log10(abs(h1)), label=l1+'(n='+str(n1)+')')
# Chebychev1 filter
n2, wn2 = signal.cheb1ord(wp, ws, gpass, gstop, fs=fs)
b2, a2 = signal.iirfilter(n2, wn2, rp=gpass, rs=gstop, btype='low', ftype=l2)
f2, h2 = signal.freqz(b2, a2, fs=fs)
ax.semilogx(f2, 20*np.log10(abs(h2)), label=l2+'(n='+str(n2)+')')
# Chebychev2 filter
n3, wn3 = signal.cheb2ord(wp, ws, gpass, gstop, fs=fs)
b3, a3 = signal.iirfilter(n3, wn3, rp=gpass, rs=gstop, btype='low', ftype=l3)
f3, h3 = signal.freqz(b3, a3, fs=fs)
ax.semilogx(f3, 20*np.log10(abs(h3)), label=l3+'(n='+str(n3)+')')
# elliptic filter
n4, wn4 = signal.ellipord(wp, ws, gpass, gstop, fs=fs)
b4, a4 = signal.iirfilter(n4, wn4, rp=gpass, rs=gstop, btype='low', ftype=l4)
f4, h4 = signal.freqz(b4, a4, fs=fs)
ax.semilogx(f4, 20*np.log10(abs(h4)), label=l4+'(n='+str(n4)+')')
plt.grid(which="both")
ax.set_title('Digital low pass filter frequency response')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Amplitude [dB]')
ax.legend(loc=0)
plt.show()
#解説
それでは，解説をしていきます．
- 正規化
ディジタルフィルタを扱う際には，正規化というものを行う必要があるようです．
正規化とは，Pythonではナイキスト周波数で割ってやって無次元量にする操作のようです．
一般には，サンプリング周波数で割るようことを指すようなので注意が必要です．
fn = fs/2
wp = fpass/fn
ws = fstop/fn
- 次数と遮断周波数の取得
今回は，次数と遮断周波数を引数に渡す必要があるため，フィルタの仕様から次数と遮断周波数を求めます．
使う関数はbutteord()です．これはmatlab由来の関数のようです．
n1, wn1 = signal.buttord(wp, ws, gpass, gstop, fs=fs)
- フィルタ設計
iirfilter()を使います．引数は先に述べたとおり．返り値は伝達関数の分子と分母になります．
rp, rsはButterWorthフィルタでは必要ないですが，そのほかでは渡してもいいです．
b1, a1 = signal.iirfilter(n1, wn1, btype='low', ftype=l1)
- 周波数とゲイン
ボード線図を描くため，周波数とゲインを取得します．ディジタルフィルターを扱っているので，freqz()を使います．引数は先に述べたとおり．返り値は伝達関数の分子と分母になります．
引数は，伝達関数の分子と分母です．
アナログフィルタであれば，freqs()を使います．
f1, h1 = signal.freqz(b1, a1, fs=fs)
- 描画
最後に描画をします．x軸を対数スケールにしたいので，semilogx()を使います．
ax.semilogx(f1, 20*np.log10(abs(h1)), label=l1+'(n='+str(n1)+')')
実行結果
#まとめ
いかがだったでしょうか？iirfilter()を用いてフィルタ設計を行ってみました．
信号処理はまだまだ奥が深そうなので，時間があればFIRフィルタなどもやってみたいです．
私自身，信号処理の知識が少し薄いので間違いなどあれば指摘していただけると助かります．
参考文献
https://docs.scipy.org/doc/scipy/reference/signal.html
http://gsmcustomeffects.hatenablog.com/entry/2019/01/16/225309
https://org-technology.com/posts/low-pass-filter.html