#概要
今回は，ディジタルIIRフィルタの設計をiirdesign()を使ってやってみようと思います．
python×信号処理の参考文献はあまり見つけることができなかったので，困ったときになにか参考になれば幸いです．
(そもそもPythonで信号処理はあまりやらないかも…？)
#動作環境
- 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$ |
#iirdesignによる実装
ではまずは，iirdesign()を使って実装していきます．
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
b1, a1 = signal.iirdesign(fpass, fstop, gpass, gstop, ftype=l1, fs=fs)
n1 = len(a1)-1
f1, h1 = signal.freqz(b1, a1, fs=fs)
ax.semilogx(f1, 20*np.log10(abs(h1)), label=l1+'(n='+str(n1)+')')
# Chebychev1 filter
b2, a2 = signal.iirdesign(fpass, fstop, gpass, gstop, ftype=l2, fs=fs)
n2 = len(a2)-1
f2, h2 = signal.freqz(b2, a2, fs=fs)
ax.semilogx(f2, 20*np.log10(abs(h2)), label=l2+'(n='+str(n2)+')')
# Chebychev2 filter
b3, a3 = signal.iirdesign(fpass, fstop, gpass, gstop, ftype=l3, fs=fs)
n3 = len(a3)-1
f3, h3 = signal.freqz(b3, a3, fs=fs)
ax.semilogx(f3, 20*np.log10(abs(h3)), label=l3+'(n='+str(n3)+')')
# elliptic filter
b4, a4 = signal.iirdesign(fpass, fstop, gpass, gstop, ftype=l4, fs=fs)
print(b4, a4)
n4 = len(a4)-1
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
- フィルタ設計
iirdesign()を使います．引数は先に述べたとおり．返り値は伝達関数の分子と分母になります．
b1, a1 = signal.iirdesign(fpass, fstop, gpass, gstop, ftype=l1, fs=fs)
- 次数取得
グラフ描画には関係はないですが参考までに次数の取得をしておきます．
0乗まであるので-1をしておきます．
n1 = len(a1)-1
- 周波数とゲイン
ボード線図を描くため，周波数とゲインを取得します．ディジタルフィルターを扱っているので，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)+')')
実行結果
#まとめ
いかがだったでしょうか？iirdesign()を用いてフィルタ設計を行ってみました．
信号処理はまだまだ奥が深そうなので，時間があれば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