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');
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 % 0059 % Copyright (C) 2011, ATR All Rights Reserved. 0060 % License : New BSD License(see VBMEG_LICENSE.txt) 0061 0062 vb_disp('--- Bayesian inverse filter calculation '); 0063 0064 % multiple bayse file for multiple conditions 0065 bayesfile = filt_parm.bayesfile; 0066 0067 if ~iscell(bayesfile), bayesfile = {bayesfile}; end; 0068 0069 Ntask = length(bayesfile); 0070 0071 if ~isempty(proj_root) 0072 for n=1:Ntask 0073 bayesfile{n} = [proj_root filesep bayesfile{n}]; 0074 end 0075 end 0076 0077 if Ntask > 1, 0078 vb_disp(['Multiple condition filters (N=' num2str(Ntask) ... 0079 ') are calculated']); 0080 end; 0081 0082 %%%%%% load VBMEG estimated result 0083 load(bayesfile{1}, 'bayes_parm','Model','vb_parm','Model_ext'); 0084 0085 %%%%%% check parameter of 'filt_parm' 0086 % Values of 'filt_parm' fields dominates over 0087 % those of 'bayes_parm' in bayesfile 0088 [bayes_parm, ix_area] = check_arg(bayes_parm, filt_parm); 0089 0090 if ~isempty(proj_root); 0091 bayes_parm_abs = vb_parm_absolute_path(proj_root, bayes_parm); 0092 else 0093 bayes_parm_abs = bayes_parm; 0094 end 0095 0096 %%%%%% Preparation of lead fields 0097 % Gact : leadfield of focal window 0098 % ix_act : Vertex index corresponding to active current Zact 0099 % ix_act_ex : Vertex index corresponding to active current Jact 0100 % Wact : Spatial smoothing matrix of focal window 0101 0102 % Focal window 0103 vb_disp('--- Preparation of lead field matrix '); 0104 0105 lf_parm.brainfile = bayes_parm_abs.brainfile; 0106 lf_parm.areafile = bayes_parm_abs.areafile; 0107 lf_parm.patch_norm = bayes_parm_abs.patch_norm; 0108 lf_parm.expand_spatial_filter = bayes_parm_abs.expand_spatial_filter; 0109 lf_parm.basisfile = bayes_parm_abs.basisfile; 0110 lf_parm.area_key = bayes_parm_abs.area_key; 0111 lf_parm.reduce = bayes_parm_abs.reduce; 0112 lf_parm.Rfilt = bayes_parm_abs.Rfilt; 0113 lf_parm.remove_area_key = []; 0114 0115 [Gact, ix_act, ix_act_ex, Wact, Lact] = ... 0116 vb_leadfield_preparation(lf_parm); 0117 0118 %%%%%% Area index in which current is calculated 0119 if ~isempty(ix_area), 0120 % Select vertex index 'ix_area' within the active current region 0121 [jx_area_ex, ix_area_ex] = vb_index2indexsub(ix_area, ix_act_ex); 0122 else 0123 jx_area_ex = 1:length(ix_act_ex); 0124 end 0125 0126 Wact = Wact(jx_area_ex,:); 0127 jx_act = find( sum(Wact, 1) > 0); 0128 Wact = Wact(:,jx_act); 0129 0130 % active index of Z-current 0131 ix_act = ix_act(jx_act); 0132 % active index of J-current 0133 ix_act_ex = ix_act_ex(jx_area_ex); 0134 0135 % # of active vertex 0136 Njact_area = length(jx_act); 0137 0138 %% Constant 0139 Nsession = length(Gact); % Number of sessions 0140 0141 if Nsession > 1, error('Multiple leadfield is not supported');end 0142 0143 G = Gact{1}; % Ga 0144 Gt = G'; 0145 0146 [Nj, Nwindow] = size(Model.a); 0147 Nch = size(Model.Cov{1},1); 0148 0149 % Extra-dipole support 0150 if ~isempty(Model_ext), 0151 vb_disp('Lead field matrix of extra-dipole '); 0152 Get = vb_load_basis(bayes_parm.extra.basisfile{1}); 0153 Ge = Get'; 0154 end; 0155 0156 %%%%%% Estimated current variance 0157 a_inv = zeros(Nj, Nwindow, Ntask); 0158 sx = zeros(Nwindow, Ntask); 0159 Covs = cell(1,Ntask); 0160 0161 for n=1:Ntask 0162 load(bayesfile{n}, 'Model'); 0163 a_inv(:,:,n) = Model.a; 0164 sx(:,n) = Model.sx; 0165 Covs{n} = Model.Cov{1}; 0166 0167 % Extra-dipole support 0168 if ~isempty(Model_ext), 0169 e_inv(:,:,n) = Model_ext.a; 0170 end 0171 end 0172 0173 Hj = zeros(Nwindow,Ntask); 0174 KW = zeros(Njact_area,Nch,Nwindow,Ntask); 0175 SB_cov = zeros(Nch,Nch,Nwindow,Ntask); 0176 0177 % Debug code 0178 %S2.a = a_inv(:,1); 0179 %S2.G = G; 0180 %S2.Covs = Covs; 0181 %save('tmp2.mat','S2'); 0182 %error('aho'); 0183 0184 if ~isempty(Model_ext), 0185 %%%%%% Time window loop 0186 for j=1:Nwindow 0187 for n=1:Ntask 0188 Cov = Covs{n}; % Sg 0189 0190 %%%% Calculate Inverse filter for current reconstruction 0191 A_z = a_inv(:,j,n); 0192 A_e = e_inv(:,j,n); 0193 0194 % Inverse filter for each session 0195 % SB = Ga * a_inv *Ga' + Cov 0196 SB = G * vb_repmultiply(Gt, A_z) + Cov; % Sigma_B 0197 SB = SB + Ge * vb_repmultiply(Get, A_e); 0198 SB_inv = pinv( SB ); % Sigam_B^{-1} 0199 SB_cov(:,:,j,n) = SB_inv; 0200 0201 Hj(j,n) = - 0.5*( vb_log_det(SB) + Nch*log(sx(j,n)) ); 0202 0203 %%%% Calculate Inverse filter for current reconstruction 0204 GSB = Gt(jx_act,:)*SB_inv; % Ga' * Sigma_B^{-1} 0205 KW(:,:,j,n) = vb_repmultiply(GSB, A_z(jx_act)); 0206 GSB_ext = Get*SB_inv; 0207 KW_ext(:,:,j,n) = vb_repmultiply(GSB_ext, A_e); 0208 % SCVS(:,:,n) = SB_inv*Cov*SB_inv; 0209 end 0210 end % Timewindow loop 0211 else 0212 for j=1:Nwindow 0213 for n=1:Ntask 0214 Cov = Covs{n}; % Sg 0215 0216 %%%% Calculate Inverse filter for current reconstruction 0217 A_z = a_inv(:,j,n); 0218 0219 % Inverse filter for each session 0220 % SB = Ga * a_inv *Ga' + Cov 0221 SB = G * vb_repmultiply(Gt, A_z) + Cov; % Sigma_B 0222 SB_inv = inv( SB ); % Sigam_B^{-1} 0223 SB_cov(:,:,j,n) = SB_inv; 0224 0225 Hj(j,n) = - 0.5*( vb_log_det(SB) + Nch*log(sx(j,n)) ); 0226 0227 %%%% Calculate Inverse filter for current reconstruction 0228 GSB = Gt(jx_act,:)*SB_inv; % Ga' * Sigma_B^{-1} 0229 KW(:,:,j,n) = vb_repmultiply(GSB, A_z(jx_act)); 0230 % SCVS(:,:,n) = SB_inv*Cov*SB_inv; 0231 end 0232 end % Timewindow loop 0233 0234 KW_ext = []; 0235 end 0236 0237 VBfilt.a = a_inv(jx_act,:,:); 0238 VBfilt.sx = sx; 0239 VBfilt.Cov = Covs; 0240 0241 VBfilt.Hj = Hj; 0242 VBfilt.KW = KW; 0243 VBfilt.KW_ext = KW_ext; 0244 VBfilt.SB_cov = SB_cov; 0245 VBfilt.ix_act = ix_act; 0246 0247 %if isfield(fileinfo,'channel_id') 0248 % VBfilt.channel_id = fileinfo.channel_id; 0249 %end 0250 0251 % Current Info 0252 Jinfo.version = vbmeg('version'); 0253 Jinfo.curr_type = 'Z-current'; 0254 0255 Jinfo.Lact = Lact; 0256 Jinfo.Wact = Wact; 0257 Jinfo.ix_act = ix_act; 0258 Jinfo.ix_act_ex = ix_act_ex; 0259 0260 % ix_act : Vertex index corresponding to active current Zact 0261 % ix_act_ex : Vertex index corresponding to active current Jact 0262 % Wact : Spatial smoothing matrix of focal window 0263 % Jact = Wact * Zact 0264 0265 return 0266 0267 %%%% --------------- 0268 function [bayes_parm, ix_area] = ... 0269 check_arg(bayes_parm,filt_parm) 0270 0271 if isfield(filt_parm,'basisfile'), 0272 bayes_parm.basisfile = filt_parm.basisfile; 0273 end; 0274 if ~isfield(filt_parm,'ix_area'), 0275 ix_area = []; 0276 else 0277 ix_area = filt_parm.ix_area; 0278 end; 0279 if isfield(filt_parm,'extra'), 0280 if isfield(filt_parm.extra,'basisfile'), 0281 bayes_parm.extra.basisfile = filt_parm.extra.basisfile; 0282 end 0283 end 0284 0285 0286 return 0287 %%% --- Following fields are used for current calculation 0288 %if isfield(filt_parm,'megfile'), 0289 % bayes_parm.megfile = filt_parm.megfile ; 0290 %end; 0291 %if isfield(filt_parm,'twin_meg'), 0292 % bayes_parm.twin_meg = filt_parm.twin_meg ; 0293 %end; 0294 %if isfield(filt_parm,'Tperiod'), 0295 % bayes_parm.Tperiod = filt_parm.Tperiod ; 0296 %end; 0297 %if isfield(filt_parm,'Tnext'), 0298 % bayes_parm.Tnext = filt_parm.Tnext ; 0299 %end; 0300 % 0301 %if ~isfield(filt_parm,'trial_average'), 0302 % trial_average = ON; 0303 %else 0304 % trial_average = filt_parm.trial_average; 0305 %end; 0306 % 0307 %bayes_parm.trial_average = trial_average; 0308 % 0309 %if ~isfield(filt_parm,'tsubsamp'), 0310 % tsubsamp = []; 0311 %else 0312 % tsubsamp = filt_parm.tsubsamp; 0313 %end; 0314 %if ~isfield(filt_parm,'overlap_mode'), 0315 % overlap_mode = 0; 0316 %else 0317 % overlap_mode = filt_parm.overlap_mode; 0318 %end;