vb_invfilter_calc_z

 Inverse filter for (Z-)current reconstruction
   [KW, KW0, Var, R] = ...
      vb_invfilter_calc(a_inv, v_inv, G, Gb, GG, Cov)

 --- Input
 a_inv : A    :  variance of active current                    
 v_inv : tau  : variance of active current                   
 G     : Ga : leadfield of active current (smoothed) G0*Wa  
 Gb    : leadfield of background current (smoothed) G0*Wb
 GG    =  Gb*Gb'
 Cov   : Sensor noise covariance matrix
 ---Output
 KW   : Inverse filter for active current
 KW0  : Inverse filter for background current
 Var  : posterior variance
   When active area is empty,
   posterior variance for background area is calculated
 R    : Estimation gainn = diag(Resolution matrix)
 --- Current reconstruction
 Zact = (KW  * Bt);
 Zbck = (KW0 * Bt);
 Original vertex current can be obtained by applying spatial smoothing filter
 Jact = Wact * Zact
 Jbck = Wbck * Zbck

 2007-3-5 M.Sato

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

0001 function [KW, KW0, Var, R] = ...
0002    vb_invfilter_calc_z(a_inv, v_inv, G, Gb, GG, Cov)
0003 % Inverse filter for (Z-)current reconstruction
0004 %   [KW, KW0, Var, R] = ...
0005 %      vb_invfilter_calc(a_inv, v_inv, G, Gb, GG, Cov)
0006 %
0007 % --- Input
0008 % a_inv : A    :  variance of active current
0009 % v_inv : tau  : variance of active current
0010 % G     : Ga : leadfield of active current (smoothed) G0*Wa
0011 % Gb    : leadfield of background current (smoothed) G0*Wb
0012 % GG    =  Gb*Gb'
0013 % Cov   : Sensor noise covariance matrix
0014 % ---Output
0015 % KW   : Inverse filter for active current
0016 % KW0  : Inverse filter for background current
0017 % Var  : posterior variance
0018 %   When active area is empty,
0019 %   posterior variance for background area is calculated
0020 % R    : Estimation gainn = diag(Resolution matrix)
0021 % --- Current reconstruction
0022 % Zact = (KW  * Bt);
0023 % Zbck = (KW0 * Bt);
0024 % Original vertex current can be obtained by applying spatial smoothing filter
0025 % Jact = Wact * Zact
0026 % Jbck = Wbck * Zbck
0027 %
0028 % 2007-3-5 M.Sato
0029 %
0030 % Copyright (C) 2011, ATR All Rights Reserved.
0031 % License : New BSD License(see VBMEG_LICENSE.txt)
0032 
0033 % Nact   # of estimated variance in focal area
0034 % Njact  # of current in focal area
0035 % Njbck  # of current in global area
0036 
0037 error('This function is a canditate to be removed, but invoked.');
0038 
0039 [Nch, Njact] = size(G);
0040 Njbck = size(Gb,2);
0041 Nact  = length(a_inv);
0042 Ratio = Njact/Nact;  % # of active current = Njact = Ratio * (# of a_inv)
0043 
0044 A_inv  = repmat(a_inv , [Ratio Nch]); % [Njact x Nch]
0045 AGt    = A_inv .* G';          % A*Ga'
0046 GGA    = G * AGt + v_inv * GG; % Ga*A*Ga' + tau*Gb*Gb'
0047 SB_inv = inv( GGA + Cov );     % SB^(-1) = [Ga*A*Ga' + tau*Gb*Gb'+Cov]^{-1}
0048 %%%% Inverse filter
0049 KW  = AGt * SB_inv;            % A*Ga'*(SB)^{-1}   : active
0050 KW0 = v_inv .* (Gb' * SB_inv); % tau*Gb'*(SB)^{-1} : back
0051 
0052 if nargout < 3, return; end;
0053 
0054 %%%% Posterior variance calculation
0055 
0056 % Posterior variance
0057 if Njact > 0
0058     if Njbck == Njact
0059         % PP = (A*Ga') + v_inv*(Gb')
0060         PP = (AGt + v_inv .* Gb' );
0061         % A + tau - (A*Ga'+tau*Gb')*SB_inv*(...)'
0062         Var = repmat(a_inv , [Ratio 1]) + v_inv - sum((PP*SB_inv).*PP ,2) ;
0063         R   = sum((KW + KW0) .* G' ,2);
0064     else
0065         % A - (A*Ga')*SB_inv*(...)'
0066         Var = repmat(a_inv , [Ratio 1]) - sum((AGt*SB_inv).*AGt ,2) ;
0067         R   = sum(KW .* G' ,2);
0068     end
0069 elseif Njbck > 0
0070     % tau - (tau*Gb')*SB_inv*(...)'
0071     Var = v_inv - v_inv^2 * sum(Gb.*(SB_inv*Gb) ,1)';
0072     R   = sum(KW0 .* Gb' ,2);
0073 end
0074