事情があってIIRフィルタのシミュレーションの必要が生じました。ツールを探してみましたが、これが見つかりません。実はあるにはあるのです。MATLAB。MATLABは確かに安くなりましたが、個人用でもフリーではありませんし、しかもSimlinkがついてきません。SimlinkのついていないMATLABかぁ、とまだ購入はしていない状態です。
じゃぁ、Scilabを。となるわけですが、ScilabはMATLABと互換性がありません。両者に互換性があるなんてのは、私に言わせれば『C++とJavaには互換性がある』と言い放つくらい恐れを知らぬ蛮行ですよ。
ScilabのIIRシミュレーションは、Scilabに固有の多項式機能を使うため、MATLABとは全く異なる方法になります。それなりに理解できたので、ここにメモ代わりに記しておきます。
そもそもIIRフィルタは間違えやすい
もともとこんな事をし始めた理由は、実は記せません。お察しください。
しかしながら、とある事情により私は発狂寸前まで追い込まれ、頭をかきむしり、口汚くPCをののしり、机をひっくり返しそうになりました。そうなった理由の一つはIIRフィルタです。IIRフィルタと言えばBi-Quad構成ですが、実はこの極めて基本的なところにすら、IIRフィルタはおおきな陥穽を用意しています。
図1をご覧ください
図1 Direct Form II型 IIRフィルタ
これはDirect Form II型のIIRフィルタです。大変よく知られた形であり、すでに大学の授業で習うなど教科書事実化している知識ですので、何のおもしろみもありません。ところが、これには別の表現があります。
図2 Direct Form II型 IIRフィルタの別表現
「別表現も何も、係数の名前を変えただけじゃないか」という指摘がありそうです。全くその通りで、図1と図2は数学的に同じです。論じる価値すらありません。
ところが工学的には大違いです。なぜなら、あなたは設計者からこんなセリフとともに係数を受け取るからです。「大丈夫、MATLABで生成したフィルタだし、シミューレーションもやったから」
a0=1.074601
a1=-0.12691
a2=0.56373
b0=1.0
b1=-0.12691
b3=0.638331
これは、図1の係数でしょうか、図2の係数でしょうか。あなたが多少なりとも信号処理になじんでいるのなら、即座に図2だと見破るでしょう。しかし、「大丈夫大丈夫、信号処理って言ってもライブラリ呼ぶだけだから」と言われてのこのこ着いてきただけのプログラマであれば、これがどちらか知ることは不可能です。
おまけにIIRの関数はわからない
仮に上の係数を図2の型だと見抜いても、与えられたIIRフィルタ関数はこんな感じだったりします。
void init_iir( TFILTER *filter, float a0, float a1, …, float b2 );
正直言ってお手上げです。どうせ関数のHELPを読んでも「a0,a1 … b2 はそれぞれIIRフィルタの係数」程度しか書いていないに決まっています。
さらに、実装によっては図1のa1, a2や図2のb1, b2の符号を反転させなければならないことがあります。
設計ソフトとライブラリの間で道に迷うことは多い
以上の事はものすごく間抜けに思えますが、実のところその辺に数えられないほど転がっています。
「MATLABで設計すれば簡単だよ」
「IIRフィルタなんて教科書に載っているんだから誰でも実装できる」
こんなことを言う人たちが作り上げる作業フローは、ほぼ1/2の確率であなたをものすごく嫌なトラブルに巻き込みます。そして、確実にそれを繰り返します。
結局、入り口が悪いのか出口が悪いのか全くわからない状態では問題は解決できません。ですので、とっかかりとして、まずは与えられた数列がIIRフィルタの係数としてどのようなものか知る必要があります。それはつまり、図1のものなのか、図2のものなのか、そして係数が反転していないかを知る必要があると言うことです。
ScilabによるIIRフィルタのシミュレーション
Scilabを使うとIIRフィルタのシミュレーションが可能です。つまり、与えられたIIRフィルタの係数から、そのフィルタの周波数特性をプロットできます。しかしその方法はMATLABとは全く違います。
Scilabの場合、IIRフィルタの伝達関数を多項式で表し、次にその多項式の周波数特性を解析してプロットします。
図1のフィルタの伝達関数は次のように表されます。
b0 + b1*z^-1 + b2*z^-2 h = -------------------------------
1 + a1*z^-1 + a2*z^-2
ここでz^-1をpと置き換えれば、何の変哲も無い二次の多項式になります。
b0 + b1*p + b2*p^2 h = -------------------------------
1 + a1*p + a2*p^2
これをScilabで実行すると次のようになります。
-->h = poly( [ 2 3 4 ], "p", "coeff" ) / poly( [1 5 6], "p", "coeff" ); h = 2
2 + 3p + 4p
-----------
2
1 + 5p + 6p
poly関数は第二引数の文字を変数とした多項式を作ります。ここではpの多項式です。第三変数に”coeff”を指定すると、関数に「第一引数のベクトルが、各項の係数だ」と宣言したことになります。第一引数に渡したベクトルからどのような多項式が作られるかは、上の例から自明です。なお、例にあるように、多項式同士の四則も可能です。
いったん伝達関数を作ったら、次にそれを可視化します。周波数特性を調べるには、frmag関数を使います。
-->[xm, fr] = frmag(h, 1024); -->plot(fr,xm);
これを実行すると次のような図が表示されます。
図3 frmagによる周波数特性の解析
このままでは見づらいので、縦軸を対数に、横軸を周波数で表示します。横軸は単にサンプル周波数を掛ければ済みます。
-->plot(fr*44100,20*log10(xm));
図4 Bode線図表示
このようにして得られた特性が、与えられたフィルタの特性と一致することを確認しておけば、問題解析をするときの強固な土台となります。
まとめ
IIRフィルタは教科書に掲載されているような基本的な知識で動作しますが、そもそも教科書によって表現方法が異なるため、実装時に不整合が起きることが多いです。そのような場合、一体何が間違いで何が間違っていないのかわからなくなります。しかし、Scilabによるシミュレーションを行うことで、与えられた形式の解釈に余地がないようにしてしまえば、そこを基準として一連のフローのどこに問題があるのか明確にすることが出来ます。