vb_dipole_estimate

 ---
 function [ Vindx, J ,emin ] = vb_dipole_estimate(bexp, basis, L)
 単一ダイポール推定

 Vindx : 推定されたダイポールのインデックス
       : ダイポールの位置座標 = V(Vindx,:) 
 J     : ダイポール電流強度
 emin  : 復元誤差

 bexp  : MEG-data (Nsensor ,Time) 
 basis : Lead Field (Ndipole, Nsensor) 
 L = 2 (電流方向)
 ---

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function [Vindx, J ,emin] = vb_dipole_estimate(bexp, basis, L)
0002 % ---
0003 % function [ Vindx, J ,emin ] = vb_dipole_estimate(bexp, basis, L)
0004 % 単一ダイポール推定
0005 %
0006 % Vindx : 推定されたダイポールのインデックス
0007 %       : ダイポールの位置座標 = V(Vindx,:)
0008 % J     : ダイポール電流強度
0009 % emin  : 復元誤差
0010 %
0011 % bexp  : MEG-data (Nsensor ,Time)
0012 % basis : Lead Field (Ndipole, Nsensor)
0013 % L = 2 (電流方向)
0014 % ---
0015 %
0016 % Copyright (C) 2011, ATR All Rights Reserved.
0017 % License : New BSD License(see VBMEG_LICENSE.txt)
0018 
0019 if nargin<3,
0020   L = 2;
0021 end;
0022 
0023 basis = basis';        % basis(Nsensor, Ndipole)
0024 
0025 [N , M] = size(basis);
0026 Npoint    = M/L;
0027 
0028 err     = zeros(Npoint,1);
0029 JJ    = zeros(Npoint,size(bexp,2));
0030 
0031 id    = 0:(L-1);
0032 id    = Npoint*id; % [ 0 , Npoint] for L=2
0033 
0034 Bnorm   = sum(sum(bexp.^2));
0035 
0036 for n = 1:Npoint,
0037   % Lead Field for the n-th dipole
0038   G  = basis(: , n + id);
0039     
0040   % Min. ( bexp - G*J )^2
0041   GG = inv( G'*G ); % (L x L) matrix inverse
0042     
0043   % Diploe current
0044   J  = GG * (G' * bexp);    % (L x T)
0045     
0046   % Reconstruction error
0047   err(n)  = sum(sum( (bexp - G * J).^2 ))/Bnorm;
0048   JJ(n,:) = sqrt(sum( J.^2, 1));
0049 end;
0050 
0051 [ emin, Vindx ] = min( err );
0052 J = JJ(Vindx,:);