L2ノルム2

通常のL2ノルムでは、脳表近くの格子点に双極子が分布するようです。 そこで深部の格子点に電流を振り分けるべく解を誘導するような重み付けをおこなってみます。 すなわち、

とします。 wにはいろいろな定義があると思います。定義は後ですることとします。
Lagrangeの未定乗数法を用い、以下のように関数Hを定義します。

q_iでHを偏微分した値がゼロになります。

まとめると、

となります。ここで

とおくと、

となります。あとはΛを求めた後、qを求めます。

ただし、ここでも

の逆行列が安定して求められるわけではありません。そこでHelsinkiに倣い特異値分解します。

右正中神経電気刺激による体性感覚誘発磁場のデータで検討してみます。
脳内の格子点の設定は、MCEを行った際に得られる
xxx-bem_pointsetXXX.mat
を用いることにします。
stacksize(50000000);
loadmatfile('.../xxx.mat');
loadmatfile('.../xxx-bem_pointsetXX.mat');
getf('.../CalcLFM.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);//unstable to provide complex number?
[Usvd,Ssvd,Vsvd]=svd(real(LFM));
n=120;//default
Usvd=Usvd(:,1:n);Ln=Ssvd(1:n,:)*Vsvd';
clear Vsvd;//recommended to save memory
Ssvd=diag(Ssvd);
sum(Ssvd(1:n))/sum(Ssvd)
ここまでは通常のL2ノルムと共通です。
次に重み付け行列

を、深部の格子点に電流をより多く振り分けるように、均質導体球モデルの中心から格子点までの距離の二乗の逆数します。
distance=2;
weight=sum((R.P-ones(size(R.P,1),1)*Center).^2,2).^(0.5*distance);
weight=diag([weight;weight]');//逆数の逆数
WMN=weight*Ln'*inv(Ln*weight*Ln')*Usvd';
// MEG data, noise data, variance.
meg=CalcSSP2(ave_0001(:,1:306));
meg=meg';time=ave_0001(:,$);
Q=WMN*meg;
これで各格子点の電流が求まりました。念のため元波形、計算波形、差分波形を見てみます。
scf();xset('use color',0);
subplot(311);plot2d(time,meg');
set(gca(),'data_bounds',[-50,-700;500,700]);
set(gca(),'tight_limits','on');
subplot(312);plot2d(time,Q'*LFM');
set(gca(),'data_bounds',[-50,-700;500,700]);
set(gca(),'tight_limits','on');
subplot(313);plot2d(time,meg'-Q'*LFM');
set(gca(),'data_bounds',[-50,-700;500,700]);
set(gca(),'tight_limits','on');

上段が元波形、中段が計算波形、下段が差分波形です。 元波形と計算波形が似ていることがわかります。通常のL2ノルムとの違うのだろうか？ とりあえず、これで準備OKです。
次に22msecと86msecの電流分布をみてみます。
t=min(find(time>=22));
//t=min(find(time>=86));
nq=size(Q,1)/2;
q=(Q(1:nq,t).^2+Q(nq+1:$,t).^2).^0.5;//unstable to cause complex number?
q=abs(q)'*1e+9;//nA
scf();set(gcf(),'visible','off');
ncolor=32;
set(gcf(),'color_map',jetcolormap(ncolor));
subplot(121);
P=R.P';
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*5)=[];PP=PP(:,n);...
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以下、赤が5nAです。なんかうまくいきませんでした。
距離の4乗で試してみました。変更点は上記の
distance=4
とするだけです。


距離の8乗で試してみました。


・・・なかなか理想どおりにいきません。
次数を30に下げてみます。
[Usvd,Ssvd,Vsvd]=svd(real(LFM));
n=30;//recommend in MCE
Usvd=Usvd(:,1:n);Ln=Ssvd(1:n,:)*Vsvd';
clear Vsvd;//recommended to save memory
Ssvd=diag(Ssvd);
sum(Ssvd(1:n))/sum(Ssvd)
distance=2;
weight=sum((R.P-ones(size(R.P,1),1)*Center).^2,2).^(0.5*distance);
weight=diag([weight;weight]');//逆数の逆数
WMN=weight*Ln'*inv(Ln*weight*Ln')*Usvd';
// MEG data, noise data, variance.
meg=CalcSSP2(ave_0001(:,1:306));
meg=meg';time=ave_0001(:,$);
Q=WMN*meg;
次に22msecと86msecの電流分布をみてみます。
t=min(find(time>=22));
//t=min(find(time>=86));
nq=size(Q,1)/2;
q=(Q(1:nq,t).^2+Q(nq+1:$,t).^2).^0.5;//unstable to cause complex number?
q=abs(q)'*1e+9;//nA
scf();set(gcf(),'visible','off');
ncolor=32;
set(gcf(),'color_map',jetcolormap(ncolor));
subplot(121);
P=R.P';
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*1)=[];PP=PP(:,n);...
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以下、赤が1nAです。次数の設定は正直よくわかんないです。