SCEファイルの作成 2

MATLABの関数filtfiltはN×Mの配列、NチャンネルのMサンプルの変数がそのまま使えるのですが、 Scilabではどうも1チャンネル×Mサンプル以外認識しないようです。そこで以下のような関数SCEファイル Filtfilt.sceを作成します。


聴覚誘発磁界のデータを読み込みます。
loadmatfile('xxxAEF.mat');
SSP=eye(306,306)-PCA_001*PCA_0001'・・・-PCA_000N*PCA_000N';
meg=ave_0001(:,1:306)*SSP;
meg(:,3:3:306)=[];//magnetometerは除去
time=ave_0001(:,$);
fs=1000*(length(time)-1)/(max(time)-min(time));//sampling frequency
hz=iir(4,'sb','ellip',[55,65]/fs,[5,5]/fs);//交流波をとるため55〜65Hzの除去
getf('.../Filtfilt.sce');//関数Filtfilt組み込み
newmeg=Filtfilt(meg,hz);
scf();
subplot(211);plot2d(time,meg,rect=[-50,-80,500,120]);xset('use color',0);xgrid();
subplot(212);plot2d(time,newmeg,rect=[-50,-80,500,120]);xset('use color',0);xgrid();

これで交流波が消えました。

ついでにiirの部分も組み込みます。以下のような関数SCEファイル FiltfiltIIR.sceを作成します。

getf('.../FiltfiltIIR.sce');//FiltfiltIIR.sceの組み込み
newmeg=FiltfiltIIR(meg,time,4,'bp','ellip',[1 30],[0.1,0.1]);
とすれば次数4の楕円フィルターを用いた1-30Hzの帯域通過フィルターがかけられた newmegがえられます。ここではtimeを引数として 標本化周波数を求めています。またtime以降の引数はiirの引数に相当しますが、 標本化周波数で割り算を行う必要はありません。
scf();
subplot(211);plot2d(time,meg,rect=[-50,-80,500,120]);xset('use color',0);xgrid();
subplot(212);plot2d(time,newmeg,rect=[-50,-80,500,120]);xset('use color',0);xgrid();
とすれば下図が得られます。