カットオフ周波数を決めれば、単純な計算だけでIIRフィルタの係数が求まる方法はないかと思っていたら、双一次変換と呼ばれる方法があったので紹介します。
今回の記事は以下の本を主に参考にしています。
双一次変換によるIIRフィルタ設計手順
双一次変換によるIIRフィルタの設計手順は以下です。
- アナログ基準ローバスフィルタを設計する。
- 周波数変換で所望のカットオフ周波数 にする。
- アナログフィルタに双一次変換を行い、伝達関数 を求める。
今回は、手始めに以下のような2次のIIRフィルタの伝達関数を作ります。
アナログ基準ローバスフィルタの設計
カットオフ角周波数が であるローバスフィルタ(LPF)を基準LPFと呼び、最初にアナログ基準LPFの伝達関数を設計します。
ちなみに、デジタル信号処理の世界では暗黙のうちに、周波数 をサンプリング周波数 で割った正規化周波数 を周波数 と呼びます。同様に、角周波数についても正規化角周波数 を角周波数 と呼びます。
アナログフィルタの設計では、 を代入すると以下の理想LPFの周波数特性を近似するような有理関数 を求めます。
有理関数というのは以下のように表せる関数です。
直接周波数特性を近似するのは難しい(周波数特性は複素数になる)ので、まず理想の2乗振幅特性を有理関数で近似します。アナログフィルタの代表であるバタワースフィルタでは以下ので近似させています。
図1をみると、確かに理想LPFに近い特性となっています。
2乗振幅特性は有理関数で表せたので、つぎに伝達関数 を求めていきます。
の平方根を求めて と置き換えればいいと思いますが、有理関数にならないので工夫が必要です。
が以下のように表せることを利用します。
は実数をフーリエ変換したときの性質となっています。この性質がないと周波数特性を逆フーリエ変換した結果が複素数になってしまいます。
に とおくと、
の極を求めてみると、
n が奇数のときは
n が偶数のときは
よって、n=2、n=3のときの の極の位置は図2のようになります。
という条件があるので、この条件を満たすように の極を選ぶ必要があります。加えて、安定なフィルタを設計するには s 平面の左半平面に極が位置する必要があるので、 の極は左半平面にある極全部となります。
したがって、n=2のときの伝達関数 は以下のようになります。
n=1 のときの伝達関数 も求めてみると以下のようになります。
これで、 アナログ基準 LPF の設計ができました。
周波数変換
フィルタのタイプは図3のように4つあります。
カットオフ角周波数が である各タイプのフィルタを設計したい場合は表のように s を変換すればよいです。
タイプ | 変換 |
---|---|
LPF | |
HPF | |
BPF | |
BRF |
ここで、 です。
を以上のLPF変換、HPF変換すると以下のようになります。
また、 を以上のBPF変換、BRF変換すると以下のようになります。
双一次変換
双一次変換は s 変数を以下の式で z 変数にする変換をいいます。
このとき、図4のようにアナログ角周波数 と デジタル角周波数 が対応するようになっています。
双一次変換の式に、、 を代入すると関係式が以下のようにわかります。
バタワースフィルタの伝達関数を双一次変換することで、 を得ることができます。ただし、双一次変換によりカットオフ周波数の位置もずれてしまうため、アナログフィルタ設計の時点で予め考慮する必要があります。
と の関係式より、デジタルフィルタのカットオフ角周波数を としたいときは、アナログフィルタの設計時ではカットオフ角周波数を以下の にすればよいです。
このことはプリワーピングと呼ばれています。
フィルタ係数の決定
双一次変換をバタワースフィルタの伝達関数に行った結果、各フィルタ係数は以下のようになります。 はプリワーピングしたカットオフ角周波数です。
プログラム
Python で記述したIIRフィルタの係数と周波数特性を求めるプログラムは以下です。
import numpy as np import scipy.signal as sg import matplotlib.pyplot as plt import sys # コマンドラインの引数を取得 args = sys.argv if len(args)!=4: print("\nusage: python main.py fs fc type") print(" fs : sampling frequency") print(" fc : cutoff frequency") print(" type : filter type (lpf, hpf, bpf, brf)") raise Exception("Argument error ") fs = int(args[1]) # サンプリング周波数 fc = float(args[2]) # カットオフ周波数 ftype = args[3] # フィルタの種類 LPF, HPF, BPF, BRF fb = 5000 # バンドパスの幅 omega_b = 2*np.pi*(fb/fs) # バンドパスの幅 omega_c = 2*np.pi*(fc/fs) # カットオフ正規化角周波数 omega_c = 2*np.tan(0.5*omega_c) # プリワーピング omega_c1 = omega_c-0.5*omega_b omega_c2 = omega_c+0.5*omega_b omega_o = np.sqrt(omega_c1*omega_c2) a = [0,0,0] b = [0,0,0] # フィルタ係数決定 if(ftype=="lpf"): b[0] = omega_c**2 b[1] = 2*omega_c**2 b[2] = omega_c**2 a[0] = omega_c**2+2*np.sqrt(2)*omega_c+4 a[1] = 2*omega_c**2-8 a[2] = omega_c**2-2*np.sqrt(2)*omega_c+4 elif(ftype=="hpf"): b[0] = 4 b[1] = -8 b[2] = 4 a[0] = omega_c**2+2*np.sqrt(2)*omega_c+4 a[1] = 2*omega_c**2-8 a[2] = omega_c**2-2*np.sqrt(2)*omega_c+4 elif(ftype=="bpf"): b[0] = 2*omega_b b[1] = 0 b[2] = -2*omega_b a[0] = omega_o**2+2*omega_b+4 a[1] = 2*omega_o**2-8 a[2] = omega_o**2-2*omega_b+4 else: # brf b[0] = omega_o**2+4 b[1] = 2*omega_o**2-8 b[2] = omega_o**2+4 a[0] = omega_o**2+2*omega_b+4 a[1] = 2*omega_o**2-8 a[2] = omega_o**2-2*omega_b+4 # 周波数特性の出力 w, h = sg.freqz(b, a) fig = plt.figure() ax = fig.add_subplot(111) ax.plot(w, np.abs(h), c="red") ax.set_ylim([0,1.2]) ax.set_xlabel("omega") ax.set_ylabel("Amplitude [dB]") fig.savefig(ftype+".png")
以下のようにコマンドを実行することで、サンプリング周波数 fs=48kHz、カットオフ周波数 fc=12kHz の LPF のフィルタ係数と振幅特性が求まります。BPF や BRF を指定したときはバンド幅は 5kHz で、第二引数の値がバンドの中心となります。
python main.py 48000 12000 lpf
周波数特性について
サンプリング周波数 fs=48kHz、カットオフ周波数 fc=12kHz を指定してさきほどのプログラムを実行した結果、図5のような特性になりました。
あまり急峻な特性ではないですが、2次の IIR フィルタなので仕方ないですかね。
白色雑音に対して各フィルタを畳み込んだときのスペクトログラムは図6のようになります。
やはり急峻ではないですが、上手く機能しているのがわかると思います。
おわりに
長くなりましたが、双一次変換によるIIRフィルタ設計についてまとめました。FPGAでフィルタを作成するときなどに利用しようと思います。