stacksize(50000000);
loadmatfile(xgetfile('.../xxx.mat'));
getf('CalcSSP2.sce');
MEG=CalcSSP2(ave_0001(:,1:306));
time=ave_0001(:,$);
scf();plot2d(time,MEG,rect=[-50,-600,500,600]);xset('use color',0);
体性感覚誘発磁場のoriginal波形です。
次にspecを使って固有値分解します。
CovMEG=mvvacov(MEG);//variance-covariance matrix
[Vpca,Dpca]=spec(CovMEG);//eigenvalue decomposition
[Dpca,index]=sort(diag(Dpca));//decreasing order
Vpca=Vpca(:,index);
PCA=Vpca(:,1)*Vpca(:,1)';
scf();plot2d(time,MEG*PCA,rect=[-50,-600,500,600]);xset('use color',0);
mvvacovは共分散行列を作る関数で、Matlabのcovに相当するのですが、 計算方法の違いによる、まるめ方の違いのためか、 微妙にmvvacovとcovの演算結果は異なるようです。
そこでMATLABのcov相当の関数 cov.sceを作成しました。
specはMatlabのeigに相当する関数です。 Dpcaが固有値からなる対角行列、 Vpca(:,n)が固有値Dpca(n,n)に対応する固有ベクトルです。 但しMatlabと異なり小さい固有値は複素数となり、 対応する固有ベクトルも複素数をふくむように計算されます。
図は最大固有値に対応する固有ベクトルをoriginalの脳磁場波形を掛け合わせて第1主成分を作り、グラフにしたものです。
PCA=zeros(306,306);
for x=1:3;PCA=PCA+Vpca(:,x)*Vpca(:,x)';end;
scf();plot2d(time,MEG*PCA,rect=[-50,-600,500,600]);xset('use color',0);
固有値の上位1〜3位に対応する固有ベクトルを用いて、第1〜3主成分を計算し、作成した図です。 original波形に似てきました。
PCA=zeros(306,306);
for x=1:7;PCA=PCA+Vpca(:,x)*Vpca(:,x)';end;
scf();plot2d(time,MEG*PCA,rect=[-50,-600,500,600]);xset('use color',0);
第1〜7主成分を用いた図です。original波形にかなり似ています。
sumDpca=sum(abs(Dpca));
for n=1:length(Dpca);if sum(abs(Dpca(1:n)))/sumDpca>0.99;break;end;end;
n
n=40
PCA=zeros(306,306);
for x=1:n;PCA=PCA+Vpca(:,x)*Vpca(:,x)';end;
scf();plot2d(time,MEG*PCA,rect=[-50,-600,500,600]);xset('use color',0);
固有値上位1〜40までで固有値の総和の99%以上となります。 第1〜40主成分で作図したグラフです。originalの波形とほぼ同じものになっています。
但し、当たり前ですが、行列のランクは40となってしまいます。
SSPはempty roomで得られた脳磁計のデータで「見た目」大きな波形の固有ベクトルを寄せ集めたものです。