IIRフィルタ
Pythonで音声信号処理(2011/05/14)
今まで作ってきたFIRフィルタ(2011/10/23)に続いて、今回は、IIR(Infinite Impulse Response)フィルタをためしてみよう。IIRフィルタでもローパス、ハイパスなどのようなフィルタが作れる。FIRフィルタよりタップ数を小さくすることができるためより高速なフィルタになるのが利点のようだ。
IIRフィルタの定義式は、
x(n)は入力信号でy(n)は出力信号。FIRフィルタとの違いは、過去の出力信号がフィードバックされる点。第1項は、FIRフィルタの式そのものだが、第2項で y(n-j) と過去の「出力信号」が使われている。FIRフィルタと同じようにZ変換してがしがし計算していくと、IIRフィルタの伝達関数 H(z) は、
となる。A(z)とB(z)はフィルタ係数のZ変換の式がそのまま入る。
IIRフィルタの設計
IIRフィルタの設計は、FIRフィルタのときのように単に逆フーリエ変換するだけではできないようで、双一次変換という手法を使う。まず、変数がsのアナログフィルタを設計し、変数がzのデジタルフィルタに変換するという方法のようだ。バターワースフィルタ H(s) というアナログフィルタを用いるのが一般的とのこと。天下りだが、二次のバターワースフィルタ(バターワース特性)は、
と変数sの関数で表される。高次のバターワースフィルタを用いるほど理想特性に近いフィルタになる。高次のバターワースフィルタの式は参考文献を参照。このバターワースフィルタの式からローパスフィルタ(LPF)、ハイパスフィルタ(HPF)、バンドパスフィルタ(BPF)、バンドストップフィルタ(BSF)をそれぞれ作るには、変数sを次のように置き換える。
式をよく見るとLPFとHPF、BPFとBSFは逆数の関係にあるようだ。これでアナログフィルタができた。次に、アナログフィルタの変数sをデジタルフィルタの変数zに置き換える。
Wikipediaだと2/Tという比例定数が付いているのだけどこれはなくてもよいのかな?参考文献だとついてない。上の式の周波数fはアナログフィルタの周波数であることに注意。デジタルフィルタの周波数とは違うようでアナログフィルタの周波数f_aとデジタルフィルタの周波数f_dの間には下の関係がある。
カットオフ周波数はデジタルフィルタの周波数で表すのが普通なので、フィルタに周波数を指定するときは上の式でアナログフィルタの周波数に変換する必要がある。で、上の式をもとにガシガシ計算してローパスフィルタ、ハイパスフィルタ、バンドパスフィルタ、バンドストップフィルタのH(z)を計算してみた。
ローパスフィルタ
がしがし計算すると最終的にH(z)は次式になる。
ここで、フィルタ係数は
となる。見るのもおぞましい式(笑)だけど導出はそんなに難しくない。普通に代入してzの次数ごとに項をまとめるだけ。最後に、分母の定数項が1になるように分子分母を割るのがミソ。フィルタ係数の分母は全部同じ、またb(0) = b(2)となっていることもわかる。また、IIRの式を見ればわかるようにa(0)は使わない(コードでは使わないけど1.0を代入)。
この式と一番はじめに書いた伝達関数の式を見比べるとP=2、Q=2であり、それぞれのフィルタ係数a(i)、b(i)がすでに求まっている!二次のバターワースフィルタからはじめたのでP=2、Q=2となり、2個前までの入力信号 x(n), x(n-1), x(n-2)、出力信号 y(n-1), y(n-2) が使われるみたいね。三次のバターワースフィルタも計算してみたらP=3、Q=3となって3個前までの信号を見てた。FIRフィルタの係数の数は桁違いに多かったのとは対照的だなぁ。
実際にプログラムを書いて試してみよう。
#coding:utf-8 import wave import struct import numpy as np from pylab import * """IIRフィルタ""" def fft(x, fs): start = 0 N = 512 # FFTのサンプル数 # 最低、N点ないとFFTできないので0.0を追加 for i in range(N): x.append(0.0) X = np.fft.fft(x[start:start+N]) freqList = np.fft.fftfreq(N, d=1.0/fs) amplitudeSpectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in X] # 波形を描画 subplot(211) plot(range(start, start+N), x[start:start+N]) axis([start, start+N, -1.0, 1.0]) xlabel("time [sample]") ylabel("amplitude") # 振幅スペクトルを描画 subplot(212) n = len(freqList) / 2 # FFTの結果は半分まで見ればOK plot(freqList[:n], amplitudeSpectrum[:n], linestyle='-') axis([0, fs/2, 0, 2]) xlabel("frequency [Hz]") ylabel("amplitude spectrum") show() def createLPF(fc): """IIR版ローパスフィルタ、fc:カットオフ周波数""" a = [0.0] * 3 b = [0.0] * 3 denom = 1 + (2 * np.sqrt(2) * np.pi * fc) + 4 * np.pi**2 * fc**2 b[0] = (4 * np.pi**2 * fc**2) / denom b[1] = (8 * np.pi**2 * fc**2) / denom b[2] = (4 * np.pi**2 * fc**2) / denom a[0] = 1.0 a[1] = (8 * np.pi**2 * fc**2 - 2) / denom a[2] = (1 - (2 * np.sqrt(2) * np.pi * fc) + 4 * np.pi**2 * fc**2) / denom return a, b def createHPF(fc): """IIR版ハイパスフィルタ、fc:カットオフ周波数""" a = [0.0] * 3 b = [0.0] * 3 denom = 1 + (2 * np.sqrt(2) * np.pi * fc) + 4 * np.pi**2 * fc**2 b[0] = 1.0 / denom b[1] = -2.0 / denom b[2] = 1.0 / denom a[0] = 1.0 a[1] = (8 * np.pi**2 * fc**2 - 2) / denom a[2] = (1 - (2 * np.sqrt(2) * np.pi * fc) + 4 * np.pi**2 * fc**2) / denom return a, b def createBPF(fc1, fc2): """IIR版バンドパスフィルタ、fc1、fc2:カットオフ周波数""" a = [0.0] * 3 b = [0.0] * 3 denom = 1 + 2 * np.pi * (fc2 - fc1) + 4 * np.pi**2 * fc1 * fc2 b[0] = (2 * np.pi * (fc2 - fc1)) / denom b[1] = 0.0 b[2] = - 2 * np.pi * (fc2 - fc1) / denom a[0] = 1.0 a[1] = (8 * np.pi**2 * fc1 * fc2 - 2) / denom a[2] = (1.0 - 2 * np.pi * (fc2 - fc1) + 4 * np.pi**2 * fc1 * fc2) / denom return a, b def createBSF(fc1, fc2): """IIR版バンドストップフィルタ、fc1、fc2:カットオフ周波数""" a = [0.0] * 3 b = [0.0] * 3 denom = 1 + 2 * np.pi * (fc2 - fc1) + 4 * np.pi**2 * fc1 * fc2 b[0] = (4 * np.pi**2 * fc1 * fc2 + 1) / denom b[1] = (8 * np.pi**2 * fc1 * fc2 - 2) / denom b[2] = (4 * np.pi**2 * fc1 * fc2 + 1) / denom a[0] = 1.0 a[1] = (8 * np.pi**2 * fc1 * fc2 - 2) / denom a[2] = (1 - 2 * np.pi * (fc2 - fc1) + 4 * np.pi**2 * fc1 * fc2) / denom return a, b def iir(x, a, b): """IIRフィルタをかける、x:入力信号、a, b:フィルタ係数""" y = [0.0] * len(x) # フィルタの出力信号 Q = len(a) - 1 P = len(b) - 1 for n in range(len(x)): for i in range(0, P + 1): if n - i >= 0: y[n] += b[i] * x[n - i] for j in range(1, Q + 1): if n - j >= 0: y[n] -= a[j] * y[n - j] return y def save(data, fs, bit, filename): """波形データをWAVEファイルへ出力""" wf = wave.open(filename, "w") wf.setnchannels(1) wf.setsampwidth(bit / 8) wf.setframerate(fs) wf.writeframes(data) wf.close() def H2(f, a, b): nume = b[0] + b[1] * np.exp(-2j * np.pi * f) + b[2] * np.exp(-4j * np.pi * f) deno = 1 + a[1] * np.exp(-2j * np.pi * f) + a[2] * np.exp(-4j * np.pi * f) val = nume / deno return np.sqrt(val.real ** 2 + val.imag ** 2) if __name__ == '__main__': wf = wave.open("whitenoise.wav", "r") fs = wf.getnframes() x = wf.readframes(wf.getnframes()) x = frombuffer(x, dtype="int16") / 32768.0 # ディジタルフィルタを設計 fc1_digital = 1000.0 # 正規化したカットオフ周波数 fc2_digital = 3000.0 # 正規化したカットオフ周波数 fc1_analog = np.tan(fc1_digital * np.pi / fs) / (2 * np.pi) fc2_analog = np.tan(fc2_digital * np.pi / fs) / (2 * np.pi) a, b = createLPF(fc1_analog) # a, b = createHPF(fc1_analog) # a, b = createBPF(fc1_analog, fc2_analog) # a, b = createBSF(fc1_analog, fc2_analog) print a, b # フィルタをかける y = iir(x, a, b) # 伝達関数Hの周波数特性 amp = [] for f in range(0, fs/2): f = float(f) / fs amp.append(H2(f, a, b)) plot(range(0, fs/2), amp) show() # 正規化前のバイナリデータに戻す y = [int(v * 32767.0) for v in y] y = struct.pack("h" * len(y), *y) # 音声を保存 save(y, fs, 16, "whitenoise2.wav")
残りのフィルタの計算結果は、式書くの面倒なので省略・・・導出して、コードと比較してみるとよいかも?参考文献にはすべての式が載っている。IIRの周波数特性を表示してみると
ちゃんとカットオフ周波数1000Hzで減衰していることがわかる。FIRよりなだらかだけど2次のフィルタならこんなもので十分なのかな?何よりFIRに比べてフィルタ係数が圧倒的に少ない(FIRでは数百個だったけどIIRはたったの5個)ので計算も高速のようだ。音声も高音のシャーがなくなってゴーになっているのが確認できた。