周波数解析2
経時的に周波数の変化を調べることがあります。
その場合は時間をずらしながらfftを行います。
loadmatfile(xgetfile('*.mat'));
meg=CalXRnge(235)*double(raw_data(235,:));
で後頭葉の脳磁図データを得たとします。
以下のようなSCEファイル
sfft.sceを作成します。
次のように入力します
getf('.../sfft.sce');
[B,F,T]=sfft(meg,2^10,CalXRnge($),'hn',2^9);
B=B.*conj(B)/2^10;
F=F'*ones(T);
T=ones(F(:,1))*T;
scf();
ここでplot3d(T,F,B);
としたいのであるが、実はこれではうまく作図できない。
少し工夫して
T=[T(1,:);T;T($,:)];
F=[F(1,:);F;F($,:)];
B=[zeros(B(1,:));B;zeros(B($,:))];
plot3d(T,F,B);
set(gca(),'cube_scaling','on');
set(gca(),'rotation_angles',[50,20]);
set(gca(),'box','off');
xtitle('frequency analysis','[sec]','[Hz]');//Z軸は設定不能
h=gca();h.z_label.visible='off';//h.z_label.text=xxxも可能だが・・・
なんのことかよくわからないので、やり直し。
[B,F,T]=sfft(meg,2^10,CalXRnge($),'hn',2^9);
for x=1:length(F);if F(x)>=4;break;end;end;
for y=x:length(F);if F(y)>=80;break;end;end;
F=F(x:y);//4〜80Hzに絞る
B=B(x:y,:);//4〜80Hzに絞る
B=B.*conj(B)/2^10;
F=F'*ones(T);
T=ones(F(:,1))*T;
scf();
T=[T(1,:);T;T($,:)];
F=[F(1,:);F;F($,:)];
B=[zeros(B(1,:));B;zeros(B($,:))];
plot3d(T,F,B);
set(gca(),'cube_scaling','on');
set(gca(),'rotation_angles',[50,20]);
set(gca(),'box','off');
h=gca();h.data_bounds=[60,0,0;0,80,1.2e+6];//x軸反転
h.tight_limits='on';
xtitle('frequency analysis','[sec]','[Hz]');
h.z_label.text='fT^2/cm^2/Hz';
0〜30秒までは閉眼、30〜60秒では開眼しています。
10Hzのピークが30秒後には小さくなっています。60Hzの商用交流波のピークは開閉眼で変化ありません
またMATLABのdaspectに相当する関数はなさそうです。
MATLABのaxes propertyと比べ、SCILABのaxes propertyの3D表示の制御は不十分なままになっています。
z軸のラベルは中央固定です。