PythonでIIRフィルタ

Pythonで信号処理シリーズ。
IIRフィルタ編。

ホントはscipyのsignalモジュールにFIRフィルタとかIIRフィルタを実現するためのクラスだかメソッドだかがあるので、それを使えばいいと言う話なのだけど。
キチンと仕組みを理解してフィルタを使う、プログラミングの勉強の一環として、個人的な趣味、などの理由からフィルタをPythonで実装してみようという趣旨。


前置きはこのあたりにして、本題。


前の記事にも書いたけどIIRフィルタは出力のフィードバックをそのシステム内部に含むデジタルフィルタ。

出力のフィードバックを含むというのは、プログラムの観点から言うと「フィルタ自身が過去の出力を保存しており、適宜それを使用する仕組みを持っている」ということを意味している。

これに対して、FIRフィルタは「過去の入力を保存しておき、それを適宜使用する」といったものだった。
このFIRフィルタにフィードバック機構を備え付けたものがIIRフィルタだ。
イメージとしては

FIRフィルタ+フィードバック機構=IIRフィルタ

といったような感じ（と、自分は捉えてる）。


それを踏まえてココを参考にIIRフィルタを以下のような形で設計してみた。





中央の加算器を挟んで左半分がFIR部分、右側がIIR部分となっている。
今回はIIRフィルタの影響のみを調べたいので、FIRフィルタの係数はインプットをそのまま出力する形となっている。


んでソースは下のようなものに。

#coding:utf-8

from scipy import *
from matplotlib.pylab import *

def IIR(array):
    output = []
    
    #係数格納行列
    FIRCoefficient = [1,0,0]
    IIRCoefficient = [0.7,0.2,0.1]
    
    for i in range(0,len(array)):
        temp = 0

    for i in range(len(array)):
        temp = 0
        for j in range(len(FIRCoefficient)):
            temp += array[i - j] * FIRCoefficient[j]
            if len(output) > j:
                temp += output[i - j - 2] * IIRCoefficient[j]
            
        output.append(temp)

    return output

Fs = 100.
t = arange(0,1,1/Fs)
f = 5

inp = []

#入力信号生成(正弦波+ノイズ)
for i in range(len(t)):
    inp.append(sin(2 * pi * f * t[i]) + rand())

outp = IIR(inp)

plot(t,inp,"b")
plot(t,outp,"r")
xlabel("time(s)")
show()

出力結果がコチラ。





…発振してしまっている。

どうやらIIRフィルタの係数設定はかなり難しいみたい。
FIRフィルタのときは直感的に係数を決めて、目で見て分かるくらいの効果があるフィルタを作ることができた。
でもIIRの場合、そううまくはいかないようで。

係数の設定や構成など、設計段階で考慮しなければいけないことが多いみたいだ。
ということで次回はデジタルフィルタ設計の勉強をすることにするか。