格子点？の取得は、 皮髄境界mesh作成 を参考にしてください。
FTPを使ってPCへ転送後、Fiff2MatFileのmeshボタンを押し、MATファイル化します。
これで準備OKです。
stacksize(50000000);
loadmatfile(xgetfile('*_mesh.mat'));
loadmatfile('.../xxx.mat');//SEF
getf('.../CalcLFM.sce');
getf('.../CalcSSP2.sce');
Ch=ConversionCoordination(MEG2HEAD,ChInform);
Center=[0,0,40]/1000;
以下のような SCEファイルを作成しました。
getf('.../CalcVectorViewLFMRadial.sce');
LFM=CalcVectorViewLFMRadial(NodesPts,Ch,Center,1);
MCEに倣い、ランクの次元を30として特異値分解します。 但し、転置して行＞列とし、economy modeのsvd(LFM,0)を使います。
LFM=LFM';
[V,S,U]=svd(LFM,0);
clear LFM;//to save memory
n=30;
Un=U(:,1:n);
Ln=(V*S(:,1:n));
LFM=Ln*inv(Ln'*Ln)*Un';
// MEG data, noise data, variance.
meg=CalcSSP2(ave_0001(:,1:306));
meg=meg';time=ave_0001(:,$);
22msec、86msecの電流源推定を行います。
t=min(find(time>=22));
//t=min(find(time>=86));
q=LFM*meg(:,t);
q=abs(q)'*1e+9;//nA
scf();set(gcf(),'visible','off');
ncolor=32;
set(gcf(),'color_map',jetcolormap(ncolor));
subplot(121);
P=NodesPts(1:3,:);
param3d(P(1,:),P(2,:),P(3,:));
h=gca();h.box='off';
h.children.mark_mode='on';h.children.mark_size=1;
h.children.foreground=1;
nq=size(NodesPts,2);
for nn=1:ncolor-1;...
PP=P;n=1:nq;n(q < nn/ncolor*0.08)=[];PP=PP(:,n);...
if size(PP,2)==0;break;end;
param3d(PP(1,:),PP(2,:),PP(3,:));...
h.children(1).mark_mode='on';h.children(1).mark_size=round(nn/32*5+3);...
h.children(1).foreground=nn+1;...
end;
h.rotation_angles=[70,200];
h.axes_visible='off';h.isoview='on';
subplot(122);Matplot((ncolor:-1:1)');
hh=gca();hh.axes_visible='off';hh.box='off';
hh.tight_limits='on';hh.axes_bounds=[0.75,0.2,0.05,0.6];
h.axes_bounds=[0,0,0.9,1];h.margins=[0,0,0,0];
set(gcf(),'visible','on');
青が0nA、赤が0.08nAです。
でもやってみます。なお、γを最大特異値の二乗の1%とします。
clear LFM Ln Sn;// to save memory;
LFM=V*(S*U');clear V;
V=LFM';U=zeros(306,306);
U(1:153,:)=V(1:153,:)*LFM;
U(154:$,:)=V(154:$,:)*LFM;
clear V;U=inv(U+S(1,1).^2/100*eye(306,306));
LFM=LFM*U;
22msec、86msecの電流源推定を行います。
t=min(find(time>=22));
//t=min(find(time>=86));
q=LFM*meg(:,t);
q=abs(q)'*1e+9;//nA
scf();set(gcf(),'visible','off');
ncolor=32;
set(gcf(),'color_map',jetcolormap(ncolor));
subplot(121);
P=NodesPts(1:3,:);
param3d(P(1,:),P(2,:),P(3,:));
h=gca();h.box='off';
h.children.mark_mode='on';h.children.mark_size=1;
h.children.foreground=1;
nq=size(NodesPts,2);
for nn=1:ncolor-1;...
PP=P;n=1:nq;n(q < nn/ncolor*0.08)=[];PP=PP(:,n);...
if size(PP,2)==0;break;end;
param3d(PP(1,:),PP(2,:),PP(3,:));...
h.children(1).mark_mode='on';h.children(1).mark_size=round(nn/32*5+3);...
h.children(1).foreground=nn+1;...
end;
h.rotation_angles=[70,200];
h.axes_visible='off';h.isoview='on';
subplot(122);Matplot((ncolor:-1:1)');
hh=gca();hh.axes_visible='off';hh.box='off';
hh.tight_limits='on';hh.axes_bounds=[0.75,0.2,0.05,0.6];
h.axes_bounds=[0,0,0.9,1];h.margins=[0,0,0,0];
set(gcf(),'visible','on');
青が0nA、赤が0.08nAです。
-50〜500ミリ秒のGIF画像を作成し、動画化しました。
scf();
ncolor=32;P=NodesPts(1:3,:);nq=size(NodesPts,2);
set(gcf(),'color_map',jetcolormap(ncolor));
for t=1:length(time);...
set(gcf(),'visible','off');
q=abs(LFM*meg(:,t))'*1e+9;...
subplot(121);...
param3d(P(1,:),P(2,:),P(3,:));...
h=gca();h.box='off';...
h.children.mark_mode='on';h.children.mark_size=1;...
h.children.foreground=1;...
for nn=1:ncolor-1;...
PP=P;n=1:nq;n(q < nn/ncolor*0.08)=[];PP=PP(:,n);...
if size(PP,2)==0;break;end;
param3d(PP(1,:),PP(2,:),PP(3,:));...
h.children(1).mark_mode='on';h.children(1).mark_size=round(nn/32*5+3);...
h.children(1).foreground=nn+1;...
end;
h.rotation_angles=[70,200];...
h.axes_visible='off';h.isoview='on';...
subplot(122);Matplot((ncolor:-1:1)');...
hh=gca();hh.axes_visible='off';hh.box='off';...
hh.tight_limits='on';hh.axes_bounds=[0.75,0.2,0.05,0.6];...
h.axes_bounds=[0,0,0.9,1];h.margins=[0,0,0,0];...
xstring(-6,-2,sprintf('%0.1fmsec',time(t)));...
set(gce(),'font_size',3);...
set(gcf(),'visible','on');...
xs2gif(get(gcf(),'figure_id'),sprintf('L2norm_%d.gif',t));...
clf();...
end;
GIF動画(3.4MB) MPEG1動画(2.6MB)
等価電流双極子の数が1〜3個程度であればgoodness of fitnessを用いて、脳活動の信号空間と雑音空間を区別しています。 逆に考えると、双極子の数がセンサー数より多い場合では、信号空間のランクの決定がgoodness of fitnessと同じようなことをしているのだと思います。