Scipyにあるscipy.signal.iirdesignの使い方と、IIRフィルタの安定性評価の方法を整理します🐜
開発環境
サンプルプログラムを実行する開発環境は次の通りです。Python3であれば動くと思います。
| 項目 | 内容 |
|---|---|
| OS | Windows 10 |
| Python | 3.6.8 |
| Numpy | 1.18.5 |
| Scipy | 1.4.1 |
| Matplotlib | 3.2.1 |
サンプルプログラム
サンプルプログラム(sample_iirdesign.py)は以下の通りです。自由にコピペして頂いて構いません。ぜひ実際に動かしてみてください。
なお、フィルタの仕様をプロットする最新のスクリプトはGitHubで公開しています。
#!/usr/bin/python3 # -*- coding: utf-8 -*- import numpy as np from scipy import signal import matplotlib.pyplot as plt def plotiirba(b, a, fs, title='', worN=512, shape=100): # ------- F特計算 # --- ナイキスト周波数計算 fn = fs / 2.0 rad_to_freq = 1 / np.pi * fn # 周波数特性（振幅特性） w, h = signal.freqz(b, a, worN) x_freq_rspns = w * rad_to_freq y_freq_rspns = 20 * np.log10(abs(h) / 1.0 + 1.0e-06) # 群遅延 w, gd = signal.group_delay((b, a), worN) x_gd = w * rad_to_freq y_gd = gd # 単位円、ゼロ点、極 theta = np.linspace(0, 2*np.pi, 360) x_unit_circle = np.cos(theta) y_unit_circle = np.sin(theta) zeros, poles, k = signal.tf2zpk(b, a) # インパルス応答 x_impluse = signal.unit_impulse(shape) y_impluse = signal.lfilter(b, a, x_impluse) # ------- グラフ表示 fig = plt.figure(figsize=(8, 10), tight_layout=True) fig.suptitle(title) # --- 周波数特性（振幅特性） ax = plt.subplot2grid((3, 2), (0, 0), colspan=2) ax.plot(x_freq_rspns, y_freq_rspns) ax.set_ylabel('Amplitude [dB]') ax.set_xlabel('Frequency [Hz]') ax.grid() # --- 群遅延 ax = plt.subplot2grid((3, 2), (1, 0), colspan=2) ax.plot(x_gd, y_gd) ax.set_ylabel('Group delay [samples]') ax.set_xlabel('Frequency [Hz]') ax.grid() # --- Poles-zeros plot ax = plt.subplot2grid((3, 2), (2, 0)) # 単位円 ax.plot(x_unit_circle, y_unit_circle, c='k', ls=':') # ゼロ点 for z in zeros: x, y = np.real(z), np.imag(z) if np.sqrt(x**2+y**2) > 1: ax.plot(x, y, marker='o', c='r', fillstyle='none') else: ax.plot(x, y, marker='o', c='k', fillstyle='none') # 極 for p in poles: x, y = np.real(p), np.imag(p) if np.sqrt(x**2+y**2) > 1: ax.plot(x, y, marker='x', c='r') else: ax.plot(x, y, marker='x', c='k') ax.axis('equal') # 正方形にする ax.set_ylabel('Imaginary part') ax.set_xlabel('Real part') ax.grid() # --- インパルス応答 ax = plt.subplot2grid((3, 2), (2, 1)) ax.plot(y_impluse, label='System') ax.plot(x_impluse, label='Unit impulse') ax.set_ylabel('Magnitude') ax.set_xlabel('Samples') ax.legend() ax.grid() return fig class TestIirDesign(): def __init__(self): # --- 波形 self.fs = 250 # サンプリング周波数 self.sec = 3 # 時間 # 抽出したい波形 self.A = 1 # 振幅 self.f = 5 # 周波数 # ノイズ｜sin波 self.A_n = 0.1 self.f_n = 50 # --- IIRフィルタ self.wp = 10 # 通過域遮断周波数[Hz] self.ws = 50 # 阻止域遮断周波数[Hz] self.gpass = 0.5 # 通過域最大損失量[dB] self.gstop = 30 # 阻止域最小減衰量[dB] self.analog = False # ディジタルフィルタなのでFalse self.ftype = 'butter' # バタワース self.output = 'ba' # baまたはtf型 def run(self): # --- 波形 # 抽出したい波形 t = np.arange(0, self.sec, 1/self.fs) wave = self.A * np.sin(2.0 * np.pi * self.f * t) # ノイズ｜sin波 noise = self.A_n * np.sin(2.0 * np.pi * self.f_n * t) # 合成波 wave2 = wave + noise # --- IIRフィルタ b, a = signal.iirdesign(self.wp, self.ws, self.gpass, self.gstop, self.analog, self.ftype, self.output, self.fs) plotiirba(b, a, self.fs) # --- フィルタリング wave3 = signal.lfilter(b, a, wave2) fig = plt.figure(figsize=(10, 3), tight_layout=True) ax = fig.add_subplot() ax.plot(wave2, label='raw') ax.plot(wave3, label='filt') ax.legend() ax.grid() plt.show() if __name__ == '__main__': proc = TestIirDesign() proc.run()
実行結果はこちらです。今回は2枚グラフを出力します。
まず、設計したIIRフィルタです。1段目がIIRフィルタの振幅特性、2段目が群遅延、3段目が安定性を示しています。 1段目と2段目について、過去記事（ディジタル信号処理｜フィルタの周波数応答をPythonで実装）を参考にしてください。
3段目の左は、設計したIIRフィルタの零点と極を単位円と一緒にプロットしています。3段目の右は、設計したIIRフィルタのインパルス応答を示しています。詳細は後述します。
フィルタリングする前(raw)とした後(filt)のグラフです。ノイズが除去されて、きれいな正弦波が抽出できていることが確認できます。
解説
サンプルプログラムでの、波形データの条件は次の通りです。
| 項目 | 内容 |
|---|---|
| サンプリング周波数[Hz] | 250 |
| データ取得時間[s] | 10 |
| 抽出したい波形 | 振幅1, 周波数5[Hz] |
| ノイズ | 振幅0.1, 周波数50[Hz] |
計算方法などは、以前の記事（ディジタル信号処理｜周波数解析(FFT)をPythonで実装）を参照してください。
以降は、サンプルプログラムの一部を抜粋して解説します。
IIRフィルタ
ディジタル信号処理 ｜FIRとIIRの違いにある通り、IIRフィルタはFIRフィルタと比較して、少ない次数で急峻な振幅特性が得られます。しかし、発散する可能性があるため、必ず極配置などで安定性を評価する必要があります。 また、線形位相特性がない欠点もありますが、丁寧に設計すれば非常に強力なフィルタです。
scipy.signal.iirdesign
ScipyにはIIRフィルタを設計できる関数がいくつか用意されています。ここでは、通過域や阻止域など比較的細かく設計できるscipy.signal.iirdesign()を採用します。
scipy.signal.iirdesign()の入力値と共に、振幅特性の例を以下に図示します。
上記を見ながら、設計するIIRフィルタの仕様を以下とします。
| 項目 | 内容 |
|---|---|
| 通過域遮断周波数(wp)[Hz] | 10 |
| 阻止域遮断周波数(ws)[Hz] | 5 |
| 通過域最大損失量(gpass)[dB] | 0.5 |
| 阻止域最小減衰量(gstop)[dB] | 30 |
これらのパラメタを入力して、フィルタ設計をします。なお、今回はバタワース型を採用したので、通過域最大損失量はほぼゼロです。
上記のグラフより、おおよそ所望のフィルタを設計できていることがわかります。1
安定性評価
過去記事（ディジタル信号処理｜周波数応答とフィルタ設計）で言及しましたが、IIRフィルタが安定である必要十分条件は、伝達関数のすべての極が単位円の内側にあることです。これをBIBO(Bounded-Input Bounded-Output)安定性を有すると言います。
また、インパルス応答をみて評価する方法もあります。今回はこの2つの手法を実装しました。
極配置
以下では、単位円と、零点と極を計算しています。設計したIIRフィルタの係数bとaを、scipy.signal.tf2zpk()に入力することで、零点と極を計算できます。ソースをみると、np.roots()を使っていることがわかります。
# 単位円、ゼロ点、極 theta = np.linspace(0, 2*np.pi, 360) x_unit_circle = np.cos(theta) y_unit_circle = np.sin(theta) zeros, poles, k = signal.tf2zpk(b, a)
以下において、計算した極を✕でプロットしています。単位円の内部にある極は黒、外部にある極を赤でプロットすること（つまり不安定）で、安定性をひと目でわかるようにしています。
# 極 for p in poles: x, y = np.real(p), np.imag(p) if np.sqrt(x**2+y**2) > 1: ax.plot(x, y, marker='x', c='r') else: ax.plot(x, y, marker='x', c='k')
インパルス応答
以下では、入力する単位インパルスと、設計したフィルタのインパルス応答を計算しています。
# インパルス応答
x_impluse = signal.unit_impulse(shape)
y_impluse = signal.lfilter(b, a, x_impluse)
以下において、計算したインパルス応答をプロットしています。今回の例では応答がゼロに収束しているので、安定性があることがわかります。
ax = plt.subplot2grid((3, 2), (2, 1)) ax.plot(y_impluse, label='System') ax.plot(x_impluse, label='Unit impulse')
おわりに
IIRフィルタの設計と実装についてまとめました。
IIRフィルタを設計するには安定性評価が必須であり、その手間からこれまで避けてきましたが、今回整理したことで今後は問題なく設計できそうです。
上記の例ではバタワース型を採用しましたが、scipy.signal.iirdesignは他のタイプも簡単に設計できるので、ぜひ試してみてください。
参考になれば幸いです(^^)
参考文献
scipy.signal.iirdesign
Scipy公式にあるscipy.signal.iirdesignの解説ページです。関数の使い方はここで確認。IIRフィルタの設計～入力画面～ 石川高専 山田洋士 研究室ホームページ Digital Filter Design Services
DF-Designというディジタル信号処理におけるフィルタ設計ができるWebアプリを公開しています。信号処理入門 (計測・制御テクノロジーシリーズ)
ディジタル信号処理を調べるときは、まずこちらを開いてます。レベルも丁度良いのでオススメです。
- 速解 電子回路―アナログ回路の基礎と設計
アナログ信号処理の基礎がここにあります。
- 計算誤差があるので、厳密には異なることに注意しましょう。↩
