Home > functions > estimation > bayes > vb_invfilter_z.m

vb_invfilter_z

PURPOSE ^

Inverse filter for Z-current reconstruction

SYNOPSIS ^

function [KW,KW_ext] = vb_invfilter_z(a_inv,G,Cov,jx_act,e_inv,Gext)

DESCRIPTION ^

 Inverse filter for Z-current reconstruction

 --- Syntax
 [KW,KW_ext] = vb_invfilter_z(a_inv,G,Cov,jx_act,e_inv,Gext)

 --- Input
 a_inv : A    :  variance of active current                    
 G     : Ga : leadfield of active current (smoothed) G0*Wa  
 Cov   : Sensor noise covariance matrix

 --- Optional input
 jx_act : selected vertex index for current estimation
          relative index within a_inv
 e_inv  : 

 ---Output
 KW    : Inverse filter for active current
 KW_ext: Inverse filter for 

 --- Current reconstruction
 Zact = (KW  * Bt);
 Original vertex current can be obtained by applying spatial smoothing filter
 Jact = Wact * Zact
 Note: smoothing filter should not be applied to extra dipole currents. 

 --- History
 2008-7-1 M.Sato
 2008-08-19 Taku Yoshioka (support extra dipole)

 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 [KW,KW_ext] = vb_invfilter_z(a_inv,G,Cov,jx_act,e_inv,Gext)
0002 % Inverse filter for Z-current reconstruction
0003 %
0004 % --- Syntax
0005 % [KW,KW_ext] = vb_invfilter_z(a_inv,G,Cov,jx_act,e_inv,Gext)
0006 %
0007 % --- Input
0008 % a_inv : A    :  variance of active current
0009 % G     : Ga : leadfield of active current (smoothed) G0*Wa
0010 % Cov   : Sensor noise covariance matrix
0011 %
0012 % --- Optional input
0013 % jx_act : selected vertex index for current estimation
0014 %          relative index within a_inv
0015 % e_inv  :
0016 %
0017 % ---Output
0018 % KW    : Inverse filter for active current
0019 % KW_ext: Inverse filter for
0020 %
0021 % --- Current reconstruction
0022 % Zact = (KW  * Bt);
0023 % Original vertex current can be obtained by applying spatial smoothing filter
0024 % Jact = Wact * Zact
0025 % Note: smoothing filter should not be applied to extra dipole currents.
0026 %
0027 % --- History
0028 % 2008-7-1 M.Sato
0029 % 2008-08-19 Taku Yoshioka (support extra dipole)
0030 %
0031 % Copyright (C) 2011, ATR All Rights Reserved.
0032 % License : New BSD License(see VBMEG_LICENSE.txt)
0033 
0034 % Nact   # of estimated variance in focal area
0035 % Njact  # of current in focal area
0036 % # of active current = Njact = Ratio * (# of a_inv)
0037 
0038 if nargin<5, e_inv = []; end
0039 
0040 Njact = size(G,2);
0041 Nact  = length(a_inv);
0042 Ratio = Njact/Nact;
0043 
0044 if Ratio > 1
0045   a_inv  = repmat(a_inv , [Ratio 1]);
0046 end
0047 
0048 if isempty(e_inv), 
0049   AGt = vb_repmultiply(G', a_inv);   % = A*Ga'
0050   EGt = [];                       % E*Gext'
0051   SB_inv = pinv( G * AGt + Cov ); % SB^(-1) = [Ga*A*Ga' + Cov]^{-1}
0052 else
0053   AGt = vb_repmultiply(G', a_inv);    % = A*Ga'
0054   EGt = vb_repmultiply(Gext', e_inv); % E*Gext'
0055   SB_inv = pinv( G * AGt + Gext * EGt + Cov ); 
0056 end
0057 
0058 if exist('jx_act','var'), 
0059   if Ratio > 1
0060     jx_act = repmat(jx_act(:),1,Ratio) ...
0061              + repmat([0:Ratio-1]*Nact,length(jx_act),1);
0062   end
0063   jx_act = jx_act(:);
0064   AGt = AGt(jx_act,:);
0065 end
0066 
0067 %%%% Inverse filter
0068 KW  = AGt * SB_inv;            % A*Ga'*(SB)^{-1}
0069 if ~isempty(EGt), KW_ext = EGt * SB_inv; 
0070 else KW_ext = []; end

Generated on Tue 27-Aug-2013 11:46:04 by m2html © 2005