Calculate Hierarchical Bayesian inverse filter
(VBMEG public function)
[syntax]
[VBfilt, Jinfo, bayes_parm, vb_parm] ...
= vbmeg_current_filter_z(proj_root, filt_parm);
[input]
proj_root: <<string>> VBMEG project root directory.
filt_parm: <<struct>> Parameters for calculation of inverse filter.
--- fields of filt_parm
bayesfile : <<string>> Current variance file (.bayes.mat).
filterfile: <optional> <<string>> inverse filter file (.vbfilt.mat). If
this field is not defined, inverse filter file is not
saved.
basisfile : <optional> <<string>> If this field is given, leadfield
file of this field is used for calculation of inverse filter
extra.basisfile : <optional> <<string>> If this field is given,
leadfield for extra dipole is given by this field
ix_area : <optional> <<vector>> Vertex indices for which inverse
filter is calculated (Z-current).
---
[output]
VBfilt: <<struct>> Inverse filter and related variables.
--- fields of VBfilt
a : Current variance (Njact,Nwindow,Ntask).
sx : Observation noise variance (Nwindow,Ntask).
Cov : Normalized sensor noise covariance matrix {Ntask}(Nch,Nch).
Hj : Model entropy (Nwindow,Ntask).
KW : VB inverse filter (Njact,Nch,Nwindow,Ntask).
KW_ext: VB inverse filter for extra-dipole.
SB_cov: Posterior sensor covariance matrix (Nch,Nch,Nwindow,Ntask).
ix_act: Vertex index corresponding to active current (Njact,1).
---
Jinfo : <<struct>> Information of cortical current.
--- fields of Jinfo
Lact : <<int>> Number of current direction.
ix_act : <<vector>> Vertex index corresponding to Z-current.
ix_act_ex: <<vector>> Vertex index corresponding to J-current.
Wact : <<matrix>> Spatial smoothing matrix of focal window.
---
[note]
J-current is calculated using Wact and Z-current
Jact = Wact * Zact
2008-10-29 M. Sato
2010-05-27 Taku Yoshioka
Extra-dipole support
2010-12-06 taku-y
[trivial] Jinfo.version = vb_version.
2010-12-07 taku-y
[trivial] Jinfo.version = vbmeg('version');
2017-10-23 Yusuke Takeda
Modified for dipoles with multiple orientations
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [VBfilt, Jinfo, bayes_parm, vb_parm] ... 0002 = vbmeg_current_filter_z(proj_root, filt_parm); 0003 % Calculate Hierarchical Bayesian inverse filter 0004 % (VBMEG public function) 0005 % 0006 % [syntax] 0007 % [VBfilt, Jinfo, bayes_parm, vb_parm] ... 0008 % = vbmeg_current_filter_z(proj_root, filt_parm); 0009 % 0010 % [input] 0011 % proj_root: <<string>> VBMEG project root directory. 0012 % filt_parm: <<struct>> Parameters for calculation of inverse filter. 0013 % --- fields of filt_parm 0014 % bayesfile : <<string>> Current variance file (.bayes.mat). 0015 % filterfile: <optional> <<string>> inverse filter file (.vbfilt.mat). If 0016 % this field is not defined, inverse filter file is not 0017 % saved. 0018 % basisfile : <optional> <<string>> If this field is given, leadfield 0019 % file of this field is used for calculation of inverse filter 0020 % extra.basisfile : <optional> <<string>> If this field is given, 0021 % leadfield for extra dipole is given by this field 0022 % ix_area : <optional> <<vector>> Vertex indices for which inverse 0023 % filter is calculated (Z-current). 0024 % --- 0025 % 0026 % [output] 0027 % VBfilt: <<struct>> Inverse filter and related variables. 0028 % --- fields of VBfilt 0029 % a : Current variance (Njact,Nwindow,Ntask). 0030 % sx : Observation noise variance (Nwindow,Ntask). 0031 % Cov : Normalized sensor noise covariance matrix {Ntask}(Nch,Nch). 0032 % Hj : Model entropy (Nwindow,Ntask). 0033 % KW : VB inverse filter (Njact,Nch,Nwindow,Ntask). 0034 % KW_ext: VB inverse filter for extra-dipole. 0035 % SB_cov: Posterior sensor covariance matrix (Nch,Nch,Nwindow,Ntask). 0036 % ix_act: Vertex index corresponding to active current (Njact,1). 0037 % --- 0038 % 0039 % Jinfo : <<struct>> Information of cortical current. 0040 % --- fields of Jinfo 0041 % Lact : <<int>> Number of current direction. 0042 % ix_act : <<vector>> Vertex index corresponding to Z-current. 0043 % ix_act_ex: <<vector>> Vertex index corresponding to J-current. 0044 % Wact : <<matrix>> Spatial smoothing matrix of focal window. 0045 % --- 0046 % 0047 % [note] 0048 % J-current is calculated using Wact and Z-current 0049 % Jact = Wact * Zact 0050 % 0051 % 2008-10-29 M. Sato 0052 % 2010-05-27 Taku Yoshioka 0053 % Extra-dipole support 0054 % 2010-12-06 taku-y 0055 % [trivial] Jinfo.version = vb_version. 0056 % 2010-12-07 taku-y 0057 % [trivial] Jinfo.version = vbmeg('version'); 0058 % 2017-10-23 Yusuke Takeda 0059 % Modified for dipoles with multiple orientations 0060 % 0061 % Copyright (C) 2011, ATR All Rights Reserved. 0062 % License : New BSD License(see VBMEG_LICENSE.txt) 0063 0064 vb_disp('--- Bayesian inverse filter calculation '); 0065 0066 % multiple bayse file for multiple conditions 0067 bayesfile = filt_parm.bayesfile; 0068 0069 if ~iscell(bayesfile), bayesfile = {bayesfile}; end; 0070 0071 Ntask = length(bayesfile); 0072 0073 if ~isempty(proj_root) 0074 for n=1:Ntask 0075 bayesfile{n} = [proj_root filesep bayesfile{n}]; 0076 end 0077 end 0078 0079 if Ntask > 1, 0080 vb_disp(['Multiple condition filters (N=' num2str(Ntask) ... 0081 ') are calculated']); 0082 end; 0083 0084 %%%%%% load VBMEG estimated result 0085 load(bayesfile{1}, 'bayes_parm','Model','vb_parm','Model_ext'); 0086 0087 Norient = vb_parm.Norient; % Number of current orientation component 0088 Norient_var = vb_parm.Norient_var;% Number of current orientation component 0089 % for variance estimation 0090 Ratio = Norient/Norient_var; % = (# of orientation) 0091 % /(# of orientation in variance) 0092 0093 %%%%%% check parameter of 'filt_parm' 0094 % Values of 'filt_parm' fields dominates over 0095 % those of 'bayes_parm' in bayesfile 0096 [bayes_parm, ix_area] = check_arg(bayes_parm, filt_parm); 0097 0098 if ~isempty(proj_root); 0099 bayes_parm_abs = vb_parm_absolute_path(proj_root, bayes_parm); 0100 else 0101 bayes_parm_abs = bayes_parm; 0102 end 0103 0104 %%%%%% Preparation of lead fields 0105 % Gact : leadfield of focal window 0106 % ix_act : Vertex index corresponding to active current Zact 0107 % ix_act_ex : Vertex index corresponding to active current Jact 0108 % Wact : Spatial smoothing matrix of focal window 0109 0110 % Focal window 0111 vb_disp('--- Preparation of lead field matrix '); 0112 0113 lf_parm.brainfile = bayes_parm_abs.brainfile; 0114 lf_parm.areafile = bayes_parm_abs.areafile; 0115 lf_parm.patch_norm = bayes_parm_abs.patch_norm; 0116 lf_parm.expand_spatial_filter = bayes_parm_abs.expand_spatial_filter; 0117 lf_parm.spatial_smoother = bayes_parm_abs.spatial_smoother; 0118 lf_parm.basisfile = bayes_parm_abs.basisfile; 0119 lf_parm.area_key = bayes_parm_abs.area_key; 0120 lf_parm.reduce = bayes_parm_abs.reduce; 0121 lf_parm.Rfilt = bayes_parm_abs.Rfilt; 0122 lf_parm.remove_area_key = []; 0123 0124 [Gact, ix_act, ix_act_ex, Wact, Lact] = ... 0125 vb_leadfield_preparation(lf_parm); 0126 0127 %%%%%% Area index in which current is calculated 0128 if ~isempty(ix_area), 0129 % Select vertex index 'ix_area' within the active current region 0130 [jx_area_ex, ix_area_ex] = vb_index2indexsub(ix_area, ix_act_ex); 0131 else 0132 jx_area_ex = 1:length(ix_act_ex); 0133 end 0134 0135 Wact = Wact(jx_area_ex,:); 0136 jx_act = find( sum(Wact, 1) > 0); 0137 Wact = Wact(:,jx_act); 0138 0139 % active index of Z-current 0140 ix_act = ix_act(jx_act); 0141 % active index of J-current 0142 ix_act_ex = ix_act_ex(jx_area_ex); 0143 0144 % # of active vertex 0145 Njact_area = length(jx_act); 0146 0147 %% Constant 0148 Nsession = length(Gact); % Number of sessions 0149 0150 if Nsession > 1, error('Multiple leadfield is not supported');end 0151 0152 G = Gact{1}; % Ga 0153 Gt = G'; 0154 0155 [Nj, Nwindow] = size(Model.a); 0156 Nch = size(Model.Cov{1},1); 0157 0158 % Extra-dipole support 0159 if ~isempty(Model_ext), 0160 vb_disp('Lead field matrix of extra-dipole '); 0161 Get = vb_load_basis(bayes_parm_abs.extra.basisfile{1}); 0162 Ge = Get'; 0163 end; 0164 0165 %%%%%% Estimated current variance 0166 a_inv = zeros(Nj, Nwindow, Ntask); 0167 sx = zeros(Nwindow, Ntask); 0168 Covs = cell(1,Ntask); 0169 0170 for n=1:Ntask 0171 load(bayesfile{n}, 'Model'); 0172 a_inv(:,:,n) = Model.a; 0173 sx(:,n) = Model.sx; 0174 Covs{n} = Model.Cov{1}; 0175 0176 % Extra-dipole support 0177 if ~isempty(Model_ext), 0178 e_inv(:,:,n) = Model_ext.a; 0179 end 0180 end 0181 0182 Hj = zeros(Nwindow,Ntask); 0183 KW = zeros(Njact_area*Ratio,Nch,Nwindow,Ntask); 0184 SB_cov = zeros(Nch,Nch,Nwindow,Ntask); 0185 0186 % Debug code 0187 %S2.a = a_inv(:,1); 0188 %S2.G = G; 0189 %S2.Covs = Covs; 0190 %save('tmp2.mat','S2'); 0191 %error('aho'); 0192 0193 % Index for dipoles with multiple orientations 0194 jx_act2=[]; 0195 for ratio=1:Ratio 0196 jx_act2=[jx_act2,jx_act+Nj*(ratio-1)]; 0197 end 0198 0199 if ~isempty(Model_ext), 0200 %%%%%% Time window loop 0201 for j=1:Nwindow 0202 for n=1:Ntask 0203 Cov = Covs{n}; % Sg 0204 0205 %%%% Calculate Inverse filter for current reconstruction 0206 A_z = a_inv(:,j,n); 0207 A_e = e_inv(:,j,n); 0208 0209 if Ratio > 1 0210 A_z = repmat(A_z,Ratio,1); 0211 end 0212 0213 % Inverse filter for each session 0214 % SB = Ga * a_inv *Ga' + Cov 0215 SB = G * vb_repmultiply(Gt, A_z) + Cov; % Sigma_B 0216 SB = SB + Ge * vb_repmultiply(Get, A_e); 0217 SB_inv = pinv( SB ); % Sigam_B^{-1} 0218 SB_cov(:,:,j,n) = SB_inv; 0219 0220 Hj(j,n) = - 0.5*( vb_log_det(SB) + Nch*log(sx(j,n)) ); 0221 0222 %%%% Calculate Inverse filter for current reconstruction 0223 GSB = Gt(jx_act2,:)*SB_inv; % Ga' * Sigma_B^{-1} 0224 KW(:,:,j,n) = vb_repmultiply(GSB, A_z(jx_act2)); 0225 GSB_ext = Get*SB_inv; 0226 KW_ext(:,:,j,n) = vb_repmultiply(GSB_ext, A_e); 0227 % SCVS(:,:,n) = SB_inv*Cov*SB_inv; 0228 end 0229 end % Timewindow loop 0230 else 0231 for j=1:Nwindow 0232 for n=1:Ntask 0233 Cov = Covs{n}; % Sg 0234 0235 %%%% Calculate Inverse filter for current reconstruction 0236 A_z = a_inv(:,j,n); 0237 if Ratio > 1 0238 A_z = repmat(A_z,Ratio,1); 0239 end 0240 0241 % Inverse filter for each session 0242 % SB = Ga * a_inv *Ga' + Cov 0243 SB = G * vb_repmultiply(Gt, A_z) + Cov; % Sigma_B 0244 SB_inv = inv( SB ); % Sigam_B^{-1} 0245 SB_cov(:,:,j,n) = SB_inv; 0246 0247 Hj(j,n) = - 0.5*( vb_log_det(SB) + Nch*log(sx(j,n)) ); 0248 0249 %%%% Calculate Inverse filter for current reconstruction 0250 GSB = Gt(jx_act2,:)*SB_inv; % Ga' * Sigma_B^{-1} 0251 KW(:,:,j,n) = vb_repmultiply(GSB, A_z(jx_act2)); 0252 % SCVS(:,:,n) = SB_inv*Cov*SB_inv; 0253 end 0254 end % Timewindow loop 0255 0256 KW_ext = []; 0257 end 0258 0259 VBfilt.a = a_inv(jx_act,:,:); 0260 VBfilt.sx = sx; 0261 VBfilt.Cov = Covs; 0262 0263 VBfilt.Hj = Hj; 0264 VBfilt.KW = KW; 0265 VBfilt.KW_ext = KW_ext; 0266 VBfilt.SB_cov = SB_cov; 0267 VBfilt.ix_act = ix_act; 0268 0269 %if isfield(fileinfo,'channel_id') 0270 % VBfilt.channel_id = fileinfo.channel_id; 0271 %end 0272 0273 % Current Info 0274 Jinfo.version = vbmeg('version'); 0275 Jinfo.curr_type = 'Z-current'; 0276 0277 Jinfo.Lact = Lact; 0278 Jinfo.Wact = Wact; 0279 Jinfo.ix_act = ix_act; 0280 Jinfo.ix_act_ex = ix_act_ex; 0281 0282 % ix_act : Vertex index corresponding to active current Zact 0283 % ix_act_ex : Vertex index corresponding to active current Jact 0284 % Wact : Spatial smoothing matrix of focal window 0285 % Jact = Wact * Zact 0286 0287 return 0288 0289 %%%% --------------- 0290 function [bayes_parm, ix_area] = ... 0291 check_arg(bayes_parm,filt_parm) 0292 0293 if isfield(filt_parm,'basisfile'), 0294 bayes_parm.basisfile = filt_parm.basisfile; 0295 end; 0296 if ~isfield(filt_parm,'ix_area'), 0297 ix_area = []; 0298 else 0299 ix_area = filt_parm.ix_area; 0300 end; 0301 if isfield(filt_parm,'extra'), 0302 if isfield(filt_parm.extra,'basisfile'), 0303 bayes_parm.extra.basisfile = filt_parm.extra.basisfile; 0304 end 0305 end 0306 0307 0308 return 0309 %%% --- Following fields are used for current calculation 0310 %if isfield(filt_parm,'megfile'), 0311 % bayes_parm.megfile = filt_parm.megfile ; 0312 %end; 0313 %if isfield(filt_parm,'twin_meg'), 0314 % bayes_parm.twin_meg = filt_parm.twin_meg ; 0315 %end; 0316 %if isfield(filt_parm,'Tperiod'), 0317 % bayes_parm.Tperiod = filt_parm.Tperiod ; 0318 %end; 0319 %if isfield(filt_parm,'Tnext'), 0320 % bayes_parm.Tnext = filt_parm.Tnext ; 0321 %end; 0322 % 0323 %if ~isfield(filt_parm,'trial_average'), 0324 % trial_average = ON; 0325 %else 0326 % trial_average = filt_parm.trial_average; 0327 %end; 0328 % 0329 %bayes_parm.trial_average = trial_average; 0330 % 0331 %if ~isfield(filt_parm,'tsubsamp'), 0332 % tsubsamp = []; 0333 %else 0334 % tsubsamp = filt_parm.tsubsamp; 0335 %end; 0336 %if ~isfield(filt_parm,'overlap_mode'), 0337 % overlap_mode = 0; 0338 %else 0339 % overlap_mode = filt_parm.overlap_mode; 0340 %end;