vb_invfilter_calc

 Calculate Inverse filter for current reconstruction
   [KW, KW0, Var] = ...
      vb_invfilter_calc(a_inv, v_inv, G, GbW, GG, Cov, Wact, Wbck)

 --- Input
 a_inv : A    :  variance of active current                    
 v_inv : tau  : variance of active current                   
 G     : Ga : leadfield of active current (smoothed) G0*Wa  
 GbW   =  Gb*Wb' = G0*Wb*Wb'
 Gb    : leadfield of background current (smoothed) G0*Wb
 GG    =  Gb*Gb'
 Wact  : Wa : spatial smoothing matrix of active current    
 Wbck  : Wb : spatial smoothing matrix of background current
 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

 --- Current reconstruction
 Jact = (KW  * Bt);
 Jbck = (KW0 * Bt);

 2006-9-1 M.Sato

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

0001 function [KW, KW0, Var] = ...
0002    vb_invfilter_calc(a_inv, v_inv, G, GbW, GG, Cov, Wact, Wbck)
0003 % Calculate Inverse filter for current reconstruction
0004 %   [KW, KW0, Var] = ...
0005 %      vb_invfilter_calc(a_inv, v_inv, G, GbW, GG, Cov, Wact, Wbck)
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 % GbW   =  Gb*Wb' = G0*Wb*Wb'
0012 % Gb    : leadfield of background current (smoothed) G0*Wb
0013 % GG    =  Gb*Gb'
0014 % Wact  : Wa : spatial smoothing matrix of active current
0015 % Wbck  : Wb : spatial smoothing matrix of background current
0016 % Cov   : Sensor noise covariance matrix
0017 % ---Output
0018 % KW   : Inverse filter for active current
0019 % KW0  : Inverse filter for background current
0020 % Var  : posterior variance
0021 %   When active area is empty,
0022 %   posterior variance for background area is calculated
0023 %
0024 % --- Current reconstruction
0025 % Jact = (KW  * Bt);
0026 % Jbck = (KW0 * Bt);
0027 %
0028 % 2006-9-1 M.Sato
0029 %
0030 % Copyright (C) 2011, ATR All Rights Reserved.
0031 % License : New BSD License(see VBMEG_LICENSE.txt)
0032 
0033 % Nvact  # of vertex in focal area
0034 % Nvex   # of extended vertex for focal area
0035 % Nvbck  # of extended vertex for global area
0036 
0037 % Nact   # of estimated variance in focal area
0038 % Njact  # of current in focal area
0039 % Njbck  # of current in global area
0040 
0041 error('This function is a canditate to be removed, but invoked.');
0042 
0043 [Nvex, Nvact] = size(Wact);
0044 [Nch,  Njact] = size(G);
0045 
0046 Nvbck = size(Wbck,1);
0047 Njbck = size(GbW,2);
0048 Nact  = length(a_inv);
0049 
0050 Ratio = Njact/Nact;  % Njact = Ratio * Nact
0051 Lact  = Njact/Nvact; % Njact = Lact  * Nvact
0052 
0053 if Nvbck > 0,
0054     Lbck  = Njbck/Nvbck; % Njbck = Lbck  * Nvbck
0055 else
0056     Lbck  = Lact;
0057 end
0058 
0059 A_inv  = repmat(a_inv , [Ratio Nch]); % [Njact x Nch]
0060 AGt    = A_inv .* G';          % A*Ga'
0061 GGA    = G * AGt + v_inv * GG; % Ga*A*Ga' + tau*Gb*Gb'
0062 SB_inv = inv( GGA + Cov );     % SB^(-1) = [Ga*A*Ga' + tau*Gb*Gb'+Cov]^{-1}
0063 % Wa*A*Ga'
0064 WAGt  = Wact * reshape(AGt, [Nvact Lact*Nch]); % [Nvex x Lact*Nch]
0065 WAGt  = reshape(WAGt, [Nvex*Lact Nch]);        % [Nvex*Lact x Nch]
0066 
0067 %switch Lact
0068 %case 1,
0069 %    WAGt  = Wact*AGt;
0070 %case 2,
0071 %    WAGt  = [Wact*AGt(1:Nvact,:); Wact*AGt(Nvact+(1:Nvact),:)];
0072 %case 3,
0073 %    WAGt  = [Wact*AGt(1:Nvact,:); Wact*AGt(Nvact+(1:Nvact),:);  ...
0074 %             Wact*AGt(2*Nvact+(1:Nvact),:)];
0075 %end
0076 
0077 %%%% Inverse filter
0078 KW  = WAGt * SB_inv;            % Wa*A*Ga'*(SB)^{-1}   : active
0079 KW0 = v_inv .* (GbW' * SB_inv); % tau*Wb*Gb'*(SB)^{-1} : back
0080 
0081 if nargout < 3, return; end;
0082 
0083 %%%% Posterior variance calculation
0084 % diagonals of WAW = diag(Wa*A*Wa')
0085 WAW  = sum( (Wact*spdiags(a_inv,0,Nvact,Nvact)).*Wact ,2 );
0086 
0087 % Posterior variance
0088 if Nvex > 0
0089     if Nvbck == Nvex
0090         % PP = (Wa*A*Ga') + v_inv*(Wb*Gb')
0091         PP = (WAGt + v_inv .* GbW' );
0092         % Wa*A*Wa' + tau*Wb*Wb' - (Wa*A*Ga'+tau*Wb*Gb')*SB_inv*(...)'
0093         Var = repmat(WAW , [Lact, 1]) ...
0094             + v_inv * repmat(sum(Wbck.*Wbck,2),[Lact,1]) ...
0095             - sum((PP*SB_inv).*PP ,2) ;
0096     else
0097         % Wa*A*Wa' - (Wa*A*Ga')*SB_inv*(...)'
0098         Var  = repmat(WAW , [Lact, 1]) ...
0099              - sum((WAGt*SB_inv).*WAGt ,2) ;
0100     end
0101 elseif Nvbck > 0
0102     % tau*Wb*Wb' - (tau*Wb*Gb')*SB_inv*(...)'
0103     Var  = v_inv * repmat(sum(Wbck.*Wbck,2),[Lbck, 1]) ...
0104          - v_inv^2 * sum(GbW.*(SB_inv*GbW) ,1)' ;
0105 end
0106