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)
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