Home > vbmeg > functions > estimation > dipole > vb_dipole_estimate_cov.m

vb_dipole_estimate_cov

PURPOSE ^

---

SYNOPSIS ^

function [Vindx, J ,emin] = vb_dipole_estimate_cov(bexp, basis, L, Cov)

DESCRIPTION ^

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

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

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

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Vindx, J ,emin] = vb_dipole_estimate_cov(bexp, basis, L, Cov)
0002 % ---
0003 % function [Vindx, J ,emin] = vb_dipole_estimate(bexp, basis, L, Cov)
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 % Cov  : Noise covariance matrix  : Cov(sensor,sensor)
0015 % ---
0016 %
0017 % Copyright (C) 2011, ATR All Rights Reserved.
0018 % License : New BSD License(see VBMEG_LICENSE.txt)
0019 
0020 if nargin<3,
0021   L = 2;
0022 end;
0023 
0024 basis    = basis';        % basis(Nsensor, Ndipole)
0025 
0026 [N , M] = size(basis);
0027 Npoint    = M/L;
0028 
0029 if nargin<4,
0030   SS       = eye(N,N);
0031 else
0032   SS       = pinv(Cov);        % Inverce of Noise covariance matrix
0033 end;
0034 
0035 err     = zeros(Npoint,1);
0036 JJ        = zeros(Npoint,L);
0037 
0038 id        = 0:(L-1);
0039 id        = Npoint*id; % [ 0 , Npoint] for L=2
0040 
0041 Bnorm   = sum(sum(bexp.^2));
0042 
0043 for n = 1:Npoint,
0044   % Lead Field for the n-th dipole
0045   G  = basis(: , n + id);
0046     
0047   % Min. ( bexp - G*J )^2
0048   GG = inv( G' * SS * G ); % (L x L) matrix inverse
0049     
0050   % Diploe current
0051   J    = GG * (G' * SS * bexp);
0052   berr = bexp - G * J;
0053     
0054   % Reconstruction error
0055   err(n)  = sum( berr' * SS * berr )/Bnorm;
0056   JJ(n,:) = sqrt(sum( J.^2, 1));
0057 end;
0058 
0059 [ emin, Vindx ] = min( err );
0060 J = JJ(Vindx,:);

Generated on Mon 22-May-2023 06:53:56 by m2html © 2005