scilabで遊ぼう (5)

今回は、ディジタル・フィルタにも使われる離散時間フィルタをつくります。

scilabは、zの式でもフィルタを記述できる

連続時間フィルタは、ラプラス変換後のsの式を使って記述しました。 同様に離散時間フィルタも、z変換後のzの式を使って記述することが出来ます。 フィルタを作成する関数は、またしても"syslin"です。

-->z=poly(0,'z');
 
-->sys=syslin('d',(1+z^(-1)+z^(-2)+z^(-3)+z^(-4))/5)
 sys  =

                     2      3      4
    0.2 + 0.2z + 0.2z + 0.2z + 0.2z
    -------------------------------
                    4
                   z

ここで例題として作成したのは、5項の移動平均フィルタです。 連続時間フィルタの時と同様にボード線図を描かせることができます。

-->scf(1);clf;bode(sys,1e-2,1e2);
WARNING: calfrq: Frequencies beyond Nyquist frequency are ignored.

100Hzまでのグラフを描かせるつもりだったのですが、「ナイキスト周波数を超えた周波数は無視します。」とメッセージが表示され、グラフは、0.5Hzで打ち切られてしまいました。

これを見ると、0.1Hz付近にカットオフ周波数のあるローパス・フィルタのような特性になっています。

scilabで、ステップ応答をみる

連続時間でのシミュレーションは、"csim"で行いましたが、離散時間でのシミュレーションは、"flts"を使います。

-->x=ones(1,101);

入力は、"1"ばかりで構成される101要素のベクトルです。

-->y=flts(x,sys);

この入力をフィルタ"sys"と共に"flts"に与えてシミュレーションを行います。

-->scf(1);clf;plot(t,x,t,y);
 
-->ax1=gca();ax1.grid=[33,33];

このように、5秒後に100%に達するするどい立ち上がり特性が得られました。

scilabで、正弦波応答をみる

離散時間シミュレーションの場合も入力に正弦波を与えて応答をみることができます。

-->t=[0:1:100];
 
-->x=sin(0.05*2*%pi*t);
 
-->y=flts(x,sys);
 
-->scf(1);clf;plot(t,x,t,y);
 
-->ax1=gca();ax1.grid=[33,33];

周波数が0.05Hzの時には、振幅は90%程度です。

-->x=sin(0.1*2*%pi*t);
 
-->y=flts(x,sys);
 
-->scf(1);clf;plot(t,x,t,y);
 
-->ax1=gca();ax1.grid=[33,33];

一方、0.1Hzの時には、60%になっています。 これも広義のローパス・フィルタになっているようです。

離散時間フィルタにおける周波数の意味

ボード線図で周波数特性をみた時に判明したように、このフィルタはサンプリング周期が1秒であると想定されています。 もっと高い周波数に対応するディジタルフィルタは、どうやったらできるのでしょうか。 実は、離散時間フィルタだけを設計する時には、サンプリング時間の絶対値には意味がないのです。

ボード線図で表示されている「周波数」は、サンプリング周波数で「正規化」された相対周波数を示しているに過ぎません。 サンプリング周波数を千倍速くすると、周波数特性も単純に千倍高くなります。 しかしながら、別のアナログ・フィルタと組み合わせた時の周波数特性を考えたときには、周波数の絶対値を気にしなくてはなりませんので、どこかでサンプリング周波数の概念を入れなくてはなりません。 これは、次回以降の課題になります。

付録 : 「scilab で遊ぼう」索引

• scilabで遊ぼう (1)
• scilabで遊ぼう (2)
• scilabで遊ぼう (3)
• scilabで遊ぼう (4)
• scilabで遊ぼう (5)
• scilabで遊ぼう (6)
• scilabで遊ぼう (7)
• scilabで遊ぼう (8)
• scilabで遊ぼう (9)
• scilabで遊ぼう (10)
• scilabで遊ぼう (11)
• scilabで遊ぼう (12)
• scilabで遊ぼう (13)