式をプログラムしただけです。妥当かどうかわかりません。
を解とするadaptive beamformerにもいろいろあって、
1) 最小分散法(minimum variance) |
2) 磁束導出行列正則化最小分散法(minimum variance with normalized lead field) |
3) 重み付け正則化最小分散法(weight-normalized minimum variance) |
があるそうです。 なお変数は特に指定しない限り、adaptive beamformerのHP と同じものを使っています。また日本語は適当につけました。
１番目の最小分散法はadaptive beamformerのHP で作りました。そこで ２番目の磁束導出行列正則化最小分散法に挑戦してみます。
stacksize(50000000);
loadmatfile('.../xxx.mat');
loadmatfile('.../xxx-bem_pointsetXX.mat');
getf('.../CalcLFM.sce');
getf('.../cov.sce');
getf('.../CalcSSP2.sce');
Ch=ConversionCoordination(MEG2HEAD,ChInform);
Center=[0,0,40]/1000;
R=CalcTangentialVectors(pointset.points,Center);
wm=1;// gradiometer : magnetometer
LFM=CalcVectorViewLFM(R,Ch,Center,wm);
ここで磁束導出行列正則化最小分散法の
を計算します。
nLFM=LFM./(ones(size(LFM,1),1)*sum(LFM,1));
ここを変えただけです。以下同様。
// MEG data, noise data, variance.
meg=CalcSSP2(ave_0001(:,1:306));
meg=meg';time=ave_0001(:,$);
Cmeg=cov(meg');
[Usvd,Ssvd,Vsvd]=svd(meg);
EigenValue=diag(Ssvd).^2;
EigenValue=EigenValue./sum(EigenValue);
for x=1:size(EigenValue,1);if sum(EigenValue(1:x))>0.90;break;end;end;
Ssvd(1:x)=0;
x
Sigma=Usvd(:,x+1:$)*Ssvd(x+1:$,:)*Vsvd';
Sigma=sum((Sigma-mean(Sigma,2)*ones(1,size(Sigma,2))).^2,2)/size(Sigma,2);
Sigma=diag(Sigma);
// Adaptive beamforming using minimum variance withnormalized lead field
mu=max(Cmeg)/10000
CS=inv(Cmeg+mu*Sigma);
H=zeros(nLFM);
Zderivate=zeros(1,size(nLFM,2));
for x=1:size(H,2);H(:,x)=CS*nLFM(:,x)/(nLFM(:,x)'*CS*nLFM(:,x));...
Zderivate(x)=((nLFM(:,x)'*Cmeg*nLFM(:,x))/(nLFM(:,x)'*Sigma*nLFM(:,x))).^0.5;...
end;
Zderivate=abs(Zderivate);
SAM.Q=H'*meg;
SAM.P=R.P;
SAM.Sigma=Sigma;
clf();set(gcf(),'visible','off');
subplot(211);plot2d(time,meg',rect=[-50,-600,500,600]);
subplot(212);plot2d(time,SAM.Q',rect=[-50,-200,500,200]);
xset('use color',0);set(gcf(),'visible','on');
最小分散法とは波形が異なるようです。・・・単位は良くわかりません。
scf();set(gcf(),'visible','off');
ncolor=32;
set(gcf(),'color_map',jetcolormap(ncolor));
subplot(121);
P=R.P';N=size(P,2);
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;
Z=(Zderivate(1:N).^2+Zderivate((N+1):2*N).^2).^0.5;
for nn=1:ncolor-1;...
PP=P;n=1:length(Z);n(Z < nn/ncolor*15)=[];PP=PP(:,n);
param3d(PP(1,:),PP(2,:),PP(3,:));...
h.children(1).mark_mode='on';...
h.children(1).mark_size=round(nn/ncolor*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');
zは変化していないようです。
parameterの設定は良くわかりません。経験則なのだと思います。