周波数フィルターの作成

かなり細かい制御ができそうですが、 とりあえず無限インパルス応答(Infinite Impulse Response: IIR)の伝達関数を iirという関数を用いて作りました。



伝達関数=iir(次数,フィルターの種類,アナログフィルターの種類,周波数,許容範囲? フィルターの種類 アナログフィルターの種類 周波数 'lp' 低域通過フィルター 'hp' 高域通過フィルター 'bp' 帯域通過フィルター 'sb' 帯域阻止フィルター 'butt' butterworthフィルター 'cheb1' Chebyshev I型フィルター 'cheb2' Chebychev II型フィルター 'ellip' 楕円フィルター 0〜0.5の値 下限 上限


とりあえず2〜40Hzの帯域通過フィルターを以下のように作成しました。
time=ave_0001(:,$); // 時間
fs=1000*(length(time)-1)/(time($)-time(1)); //標本化周波数
hz=iir(4,'bp','ellip',[2,40]/fs,[0.5,5]/fs);
//4次の楕円フィルター 周波数幅0.5、5Hz？
伝達関数hzは変数zの有理多項式となっています。
hz=
0.0175310-0.0760732z+0.1534536z²・・・
--------------------------------------
0.3178548-2.827968z+11.132419z²・・・
このhzの周波数応答は
[hw,w]=frmag(hz,2^10);
scf(1); //新しいfigure作成　MATLABのfigure(1)
// 注 SCILABでfigure(1)としない方が無難です。
plot2d(w*fs,hw,rect=[0,0,100,1]); //0-100Hzで表示
xgrid(); xtitle('magnitude of frequency response','Hz'); //文字追加
h=gca(); //ハンドル取得
h.title.font_size=3; //タイトルを大きく
h.x_label.font_size=3; //Hzを少し大きく

で調べることができますが、位相の遅延などを調べるにはrepfreqを使います。
このrepfreqでは伝達関数ではなく状態空間表現を引数とします。
と書いたものの、私自身は状態空間表現が何であるかは理解していません。 使って感覚で理解するということにして
hzs=tf2ss(hz); //状態空間表現に変換
[w,hw]=repfreq(hzs); //MATLABのfreqzに相当。
scf(2); //新しいfigure作成
subplot(211);plot2d(w*fs,abs(hw),rect=[0,0,fs/2,1]);
xgrid();xtitle('Magnitude','Hz');
h=gca();h.title.font_size=3;h.x_label.font_size=2;
[db,phi]=dbphi(hw);
subplot(212);plot2d(w*fs,phi,rect=[0,-300,fs/2,300]);
xgrid();xtitle('Phase','Hz');
h=gca();h.title.font_size=3;h.x_label.font_size=2;
とすれば下図が表示されます。


実はこんな面倒なことをしなくても
clf();
bode(tf2ss(hz)); //ボード図をプロット
とすればボード図が得られます。

但しこのbodeは値を返しません。そこで
h=gcf();//ハンドル取得
mag=h.children(2).children(2).children.data;//上のaxes.agregation.polyplot.data
phi=h.children(1).children(2).children.data;//下のaxes.agregation.polyplot.data
w=mag(:,1);
mag=10^(mag(:,2)/20);
phi=phi(:,2);
とすれば値が得られます。