dataがN次元配列のとき、fft(data,-1)はN次元フーリエ変換 MATLABのfft(fft(fft(data,[],1),[],2),[],3)等の値を返します。 MATLABのfft(data,[],dim)に相当するのが、fft(data,-1,dim,1) の筈なのですが、どうも返ってくる値が意味不明なので私は用いていません。
安静閉眼時の後頭葉の周波数を調べます。
loadmatfile(xgetfile('*.mat')); // raw fiffを読み込む
for x=1:306;SnsName=string(ChanName(x));...//セル配列を文字列に
if grep(SnsName,'MEG 2113');break;end;end; //MEG2113のチャンネル番号を取得
sfreq=CalXRnge($);//標本化周波数
wl=2^10;//波形長さ
meg=CalXRnge(x)*double(raw_data(x,((1:wl)+2000))); //MEG2113チャンネルの波形
time=((1:wl)+2000-1)/sfreq;
scf();plot2d(time,meg,rect=[time(1),-800,time($),300]);
set(gca(),'grid',[0,0]);set(gca(),'tics_color',2);//青グリッド(default状態)
xtitle(sprintf('MEG2113 %0.1f - %0.1f sec',time(1),time($)),...
'[sec]','[fT/cm]');//タイトル表示
この波形の周波数を調べます。
data=fft(meg,-1);
data=data(1:(wl/2+1));
data=data.*conj(data)/wl;//power値
Hz=(0:wl/2)/wl*sfreq;
scf();
plot2d(Hz,data);
set(gca(),'data_bounds',[0,10;100,1e+7]);//;に注意
set(gca(),'log_flags','nl');//log表示
set(gca(),'tight_limits','on');
set(gca(),'grid',[0,0]);set(gca(),'tics_color',2);//グリッドの色を青(default状態)に
xtitle('frequency with eyes closed','[Hz]','[fT^2/cm^2/Hz]');
60Hzのピークは60Hzの商用交流波です。10Hzにピークがあり、後頭葉のα波と考えられます。 シールドルームの性能のためか4Hz以下では雑音が多いように感じます。
Fourier変換は周期的な波形を想定しており、窓関数を用いて波形の端を整えます。
win=window('hn',wl);//hanning窓
megwin=meg.*win;//hanning窓をかける。
scf();plot2d(time',[meg;megwin]',rect=[time(1),-800,time($),300]);//転置
h=gca();h.children.children(1).foreground=5;//Axex.Agregation.Polylineの色を赤に
h.grid=[0,0];h.tics_color=2;//青グリッド
xtitle('effect of hanning window','[sec]','[fT/cm]');
黒がoriginalの波形、赤がhanning窓をかけた後の波形です。 赤の波形に対してFourier変換をおこないます。
data2=fft(megwin,-1);
data2=data2(1:(wl/2+1));data2=data2.*conj(data2)/wl;
scf();
plot2d(Hz',[data;data2]',rect=[0,10,100,1e+7]);//転置
h=gca();h.children.children(1).foreground=5;
h.grid=[0,0];h.tics_color=2;//青グリッド
h.log_flags='nl';//対数表示
xtitle('frequency with eyes closed processed hanning window','[Hz]','[fT^2/cm^2/Hz]');
60Hz近傍の商用交流波の影響が小さくなっています。