- Robinson S.E. and Vrba J. (1999) Functional neuroimaging by synthetic aperture Magnetometry (SAM). Recent advances in biomagnetism 302-305
- 眞溪歩 (2000) リードフィールドを用いたMEG逆問題の比較 日本生体磁気学会誌特別号 13(1):42-45
式を自分の解釈でscilabのcodeにしただけです。間違っているかもしれません。 天下に恥をさらしているような気もしますが、何かのたたき台にしていただければ幸いに存じます。
格子点の数をN、時間サンプル数をT、センサーの数をMとし、以下の行列を考えます。
ここで、
Q は脳内の格子点の２つの接平面方向の電流の大きさの時間的変化を、
B はセンサーの信号の時間的変化を、
H はB からQ を計算する導出行列を、
L は磁場導出行列であるとします。
H の定義は
です。ここで任意の格子点の１つの接平面方向について着目し、
とおくと、
となります。そこで、
という制約下で、
となるような
を求めます。なお 格子点における時間的電気活動の大きさの和が最小となる生理学的根拠はありません。 Lagrangeの未定係数法を使って解きます。
とおきます。次に
の偏微分の値がゼロになることを利用します。
すなわち
となります。
の結果を並べて整理すると
となることから
を求めると
となります。あとは簡単です。
お疲れ様でした。
但し、
が常にフルランクで逆行列が安定して求めらるわけではないので、以下のような Backus-Gilbertの適切化なるものを行うそうです。
ここで
μはRobinson先生の抄録では-1〜∞の値となっていますが、何がいいのかよくわかりません。
Σは脳磁図信号を特異値分解し、累積寄与率の低い成分で作成した雑音信号の、 各センサーにおける標本分散（不偏分散でない）からなる対角行列です。 雑音信号に用いる特異値の数は記載がないためよくわかりません。
行列の掛け算が分母になっていますが、(1×M)(M×M)(M×1)行・列の掛け算であり、 うまい具合に１つの数値に収まります。 脳内の任意の位置に設定した格子点・電流の向き毎この式を解けば、 目的とする行列H が求まり、Q も求まります。
SAMでは信号の評価量、雑音の評価量、そのz偏差を
と計算し、MRI画像などに重ね合わせているようです。
右正中神経電気刺激による体性感覚誘発磁場のデータで検討してみます。
脳内の格子点の設定は、MCEを行った際に得られる
xxx-bem_pointsetXXX.mat
を用いることにします。
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);
球の接平面方向のベクトルを表示します。 結構時間がかかります。
scf();set(gcf(),'visible','off');
P=R.P';
V2=R.V2';
V3=R.V3';
param3d(P(1,:),P(2,:),P(3,:));
h=gca();h.box='off';xset('use color',0);
h.children.mark_mode='on';h.children.mark_size=3;
PP=[P;P+0.01*V2];
for x=1:size(P,2);param3d(PP([1,4],x),PP([2,5],x),PP([3,6],x));end;
PP=[P;P+0.01*V3];
for x=1:size(P,2);param3d(PP([1,4],x),PP([2,5],x),PP([3,6],x));end;
h.box='off';h.axes_visible='off';h.isoview='on';
set(gcf(),'visible','on');
wm=1;// gradiometer : magnetometer
LFM=CalcVectorViewLFM(R,Ch,Center,wm);
// MEG data, noise data, variance.
meg=CalcSSP2(ave_0001(:,1:306));
meg=meg';time=ave_0001(:,$);
Cmeg=cov(meg');//something wrong! in intrinsic "mvvacov" function
[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;
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
for x=-15:15;mu=10^x;...
if rank(Cmeg+mu*Sigma)==size(Sigma,1);break;end;end;
mu
CS=inv(Cmeg+mu*Sigma);
H=zeros(LFM);
Zderivate=zeros(1,size(LFM,2));
for x=1:size(H,2);H(:,x)=CS*LFM(:,x)/(LFM(:,x)'*CS*LFM(:,x));...
Zderivate(x)=((LFM(:,x)'*Cmeg*LFM(:,x))/(LFM(:,x)'*Sigma*LFM(:,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,-2E-14,500,2E-14]);
xset('use color',0);set(gcf(),'visible','on');
上段が実際のセンサー波形で、下段が各格子点における波形です。
尚、固有値の累積寄与率90%の信号空間(固有値数8)を脳の信号とし、残りを雑音としました。 またBackus-Gilbertのパラメータμは逆行列が存在する最小の10n(nは整数10-9)としています。
格子点波形の分散が最小になるように解を求めるためか、低周波成分は消え、高周波成分が目立つような気がします。
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');
？？？ よくわかりません。
雑音信号空間を構成する固有値の数、Backus-Gilbertの適切化の係数といったparameterを、 どう設定すればよいのかはよくわかりません。 raw dataのような、もっとサンプル時間が多いdataだと信号空間の分散・共分散行列が小さくなり、いまとは違ったものが見えるのかもしれません。
そこで、Backus-Gilbert適切化parameterを計測データの共分散行列最大値の1万分の1にしてみました。
mu=max(Cmeg)/10000
CS=inv(Cmeg+mu*Sigma);
H=zeros(LFM);
Zderivate=zeros(1,size(LFM,2));
for x=1:size(H,2);H(:,x)=CS*LFM(:,x)/(LFM(:,x)'*CS*LFM(:,x));...
Zderivate(x)=((LFM(:,x)'*Cmeg*LFM(:,x))/(LFM(:,x)'*Sigma*LFM(:,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,-4E-9,500,4E-9]);
xset('use color',0);set(gcf(),'visible','on');
さっきよりはマシになった気がします。Backus-Gilbert適切化parameterは3.3333104・・・です。
累積寄与率は90%のままで、変化ないのでz偏差は変化しません。
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');
変化ありません。