手軽にフィルタ設計

はじめに

　SciPyのSignalサブパッケージを用いて、IIRフィルタを設計する際の基礎的なメモです。今回は、バタワース特性を持つローパスフィルタを設計してみます。

　素人理解なので、間違いがありましたらご指摘いただければ幸いです。

IIRフィルタの設計[1-3]

　以下のような手順に沿って行われます。

既存のアナログフィルタから、要求する周波数特性を満たせそうなものを探す。
アナログフィルタの伝達関数 $G(s)$ に、双一次変換法を施しデジタルフィルタの伝達関数 $H(z)$ を得る。

　今回設計するバタワース特性を有するローパスフィルタの $H(z)$ は、下記のようになります（2次系）。

　 $\displaystyle H(z) = \frac{b_{0} + b_{1}z^{-1} + b_{2} z^{-2}}{1 + a_{1} z^{-1} +a_{2} z^{-2}}$

　ここで、 $a_{n}$ 及び $b_{n}$ はフィルタ系数という、カットオフ周波数とサンプリング周波数に依存した値になります。 $a_{n}$ 及び $b_{n}$ を含めた $H(z)$ の形が定まれば、これを差分方程式に直して入力信号に対するフィルタリング処理が可能になります。また、 $\displaystyle z = e^{j \omega T}$ とした時、 $\displaystyle |H(e^{j \omega T})|$ と $\displaystyle arg(H(e^{j \omega T}))$ がそれぞれ、振幅特性・位相特性を示します（ $\displaystyle H(e^{j \omega T})$ は、一般に複素関数。 $\omega$ と $T$ は、角周波数とサンプリング周期）。

　結局は、既存の $G(s)$ から任意のカットオフ周波数における、フィルタ係数を含めた $H(z)$ 全体を求めることになるのですが、まっとうにこの変換作業を行うのはなかなか骨の折れる作業となります。

SciPyでデジタルフィルタの伝達関数を求める[4, 5]

　Pythonの外部パッケージであるSciPyには信号処理用の関数がまとまったSignalというサブパッケージが用意されています。これを利用することで手軽に先の $\displaystyle H(e^{j \omega T})$ を求めることができます。

　今回使用するSignalサブパッケージの関数は、"iirfilter"と"freqz"の2つになります。前者でフィルタ係数を、後者で $\displaystyle H(e^{j \omega T})$ を求めることができます。詳細は、後述のソースコード内のコメント及び、公式のドキュメント[5]を参照してください。

　実際にこれらを使って、バタワース特性を持つローパスフィルタ（2次系）の振幅特性をプロットしてみたのが下記のソースコードになります。

以下、実行結果です（振幅特性のボード線図、フィルタ係数）。

f:id:aogiri_m2d:20170711013037p:plain

f:id:aogiri_m2d:20170711013117p:plain

　ちゃんと、カットオフ周波数 1 kHzから減衰しています。今回は試していませんが、"iirfilter"の返り値には、極と零点を指定することもできるみたいです。

　今回は2次系で行いましたが、高次系でも、同様に簡単にフィルタ係数を求めることができるのでありがたいですね。

参考文献

三谷政昭. 信号解析のための数学. 森北出版, 1990
大類重範. ディジタル信号処理. 日本理工出版会, 2001
鏡慎吾. "やる夫で学ぶディジタル信号処理".

http://www.ic.is.tohoku.ac.jp/~swk/lecture/yaruodsp/preface.html

(参照 2017-07-11)
中久喜健司. 科学技術計算のための Python入門. 技術評論社, 2016
The Scipy community. "Signal processing (scipy.signal)".

https://docs.scipy.org/doc/scipy-0.19.0/reference/signal.html (参照 2017-07-11)

情弱

情弱のVST作製メモ

手軽にフィルタ設計

はじめに

IIRフィルタの設計[1-3]

SciPyでデジタルフィルタの伝達関数を求める[4, 5]