Description of vbmeg_current_filter_z

0001 function [VBfilt, Jinfo, bayes_parm, vb_parm] ...
0002           = vbmeg_current_filter_z_meeg(proj_root, filt_parm);
0003 % Calculate Hierarchical Bayesian inverse filter for MEG+EEG
0004 % (VBMEG public function)
0005 % The scale of inverse filter (KW) is original digit. (not division by bsnorm_meeg)
0006 % Usage: J = VBfilt.KW*[B; E];  %B and E are MEG and EEG data.
0007 %
0008 % [syntax]
0009 % [VBfilt, Jinfo, bayes_parm, vb_parm] ...
0010 %    = vbmeg_current_filter_z(proj_root, filt_parm);
0011 %
0012 % [input]
0013 % proj_root: <<string>> VBMEG project root directory.
0014 % filt_parm: <<struct>> Parameters for calculation of inverse filter.
0015 % --- fields of filt_parm
0016 %  bayesfile : <<string>> Current variance file (.bayes.mat).
0017 %  filterfile: <optional> <<string>> inverse filter file (.vbfilt.mat). If
0018 %              this field is not defined, inverse filter file is not
0019 %              saved.
0020 %  basisfile_meg : <optional> <<string>> If this field is given, leadfield
0021 %              file of this field is used for calculation of inverse filter.
0022 %  basisfile_eeg : <optional> <<string>> same as above.
0023 %  extra.basisfile : <optional> <<string>> If this field is given,
0024 %              leadfield for extra dipole is given by this field
0025 %  ix_area   : <optional> <<vector>> Vertex indices for which inverse
0026 %              filter is calculated (Z-current).
0027 % ---
0028 %
0029 % [output]
0030 % VBfilt: <<struct>> Inverse filter and related variables.
0031 % --- fields of VBfilt
0032 %  sx    : Observation noise variance (Nwindow,Ntask).
0033 %  Cov   : Normalized sensor noise covariance matrix {Ntask}(Nch,Nch).
0034 %  KW    : VB inverse filter  (Njact,Nch,Nwindow,Ntask).
0035 %  KW_ext: VB inverse filter for extra-dipole.
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 % 2017-03-13 Yusuke Fujiwara
0060 %  MEG+EEG
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 %%%%%% check parameter of 'filt_parm'
0088 % Values of 'filt_parm' fields dominates over
0089 %   those of 'bayes_parm' in bayesfile
0090 [bayes_parm, ix_area] = check_arg(bayes_parm, filt_parm);
0091             
0092 if ~isempty(proj_root);
0093   bayes_parm_abs = vb_parm_absolute_path(proj_root, bayes_parm);
0094 else
0095   bayes_parm_abs = bayes_parm;
0096 end
0097 
0098 %%%%%% Preparation of lead fields
0099 % Gact   : leadfield of focal window
0100 % ix_act : Vertex index corresponding to active current Zact
0101 % ix_act_ex : Vertex index corresponding to active current Jact
0102 % Wact   : Spatial smoothing matrix of focal window
0103 
0104 % Focal window
0105 vb_disp('--- Preparation of lead field matrix ');
0106 
0107 %MEG leadfield
0108 lf_parm.basisfile = [proj_root filesep bayes_parm_abs.basisfile_meg];
0109 
0110 lf_parm.brainfile = bayes_parm_abs.brainfile;
0111 lf_parm.areafile = bayes_parm_abs.areafile;
0112 lf_parm.patch_norm = bayes_parm_abs.patch_norm;
0113 lf_parm.expand_spatial_filter = bayes_parm_abs.expand_spatial_filter;
0114 lf_parm.spatial_smoother = bayes_parm_abs.spatial_smoother;
0115 lf_parm.area_key = bayes_parm_abs.area_key;
0116 lf_parm.reduce = bayes_parm_abs.reduce;
0117 lf_parm.Rfilt = bayes_parm_abs.Rfilt;
0118 lf_parm.remove_area_key = [];
0119 
0120 [Gact_b, ix_act, ix_act_ex, Wact, Lact]  =  ...
0121      vb_leadfield_preparation(lf_parm);
0122 
0123  
0124 %EEG leadfield
0125 lf_parm.basisfile = [proj_root filesep bayes_parm_abs.basisfile_eeg]; 
0126 Gact_e = vb_leadfield_preparation(lf_parm);
0127 
0128 for n = 1:length(Gact_b)
0129     Gact{n} = [Gact_b{n}./Model.bsnorm_meg; Gact_e{n}./Model.bsnorm_eeg];
0130     Gact0{n} = [Gact_b{n}; Gact_e{n}];
0131 end
0132 
0133 
0134 %%%%%% Area index in which current is calculated
0135 if ~isempty(ix_area),
0136   % Select vertex index 'ix_area' within the active current region
0137   [jx_area_ex, ix_area_ex] = vb_index2indexsub(ix_area, ix_act_ex);
0138 else
0139   jx_area_ex = 1:length(ix_act_ex);
0140 end
0141 
0142 Wact   = Wact(jx_area_ex,:);
0143 jx_act = find( sum(Wact, 1) > 0);
0144 Wact   = Wact(:,jx_act);
0145 
0146 % active index of Z-current
0147 ix_act = ix_act(jx_act);
0148 % active index of J-current
0149 ix_act_ex = ix_act_ex(jx_area_ex);
0150 
0151 % # of active vertex
0152 Njact_area = length(jx_act);
0153 
0154 %% Constant
0155 Nsession = length(Gact);     % Number of sessions
0156 
0157 if Nsession > 1, error('Multiple leadfield is not supported');end
0158 
0159 G = Gact{1};    % Ga
0160 Gt = G';
0161 
0162 [Nj, Nwindow] = size(Model.a);
0163 Nch = size(Model.Cov{1},1);
0164 
0165 % Extra-dipole support
0166 if ~isempty(Model_ext), 
0167   vb_disp('Lead field matrix of extra-dipole ');
0168   Get = vb_load_basis(bayes_parm_abs.extra.basisfile{1});
0169   Ge = Get';
0170 end; 
0171 
0172 %%%%%% Estimated current variance
0173 a_inv = zeros(Nj, Nwindow, Ntask);
0174 sx    = zeros(Nwindow, Ntask);
0175 Covs  = cell(1,Ntask);
0176 
0177 for n=1:Ntask
0178   load(bayesfile{n}, 'Model');
0179   a_inv(:,:,n) = Model.a;
0180   sx(:,n)    = Model.sx;
0181   Covs{n}  = Model.Cov{1};
0182   
0183   % Extra-dipole support
0184   if ~isempty(Model_ext), 
0185     e_inv(:,:,n) = Model_ext.a;
0186   end
0187 end
0188 
0189 Hj = zeros(Nwindow,Ntask);
0190 KW = zeros(Njact_area,Nch,Nwindow,Ntask);
0191 SB_cov = zeros(Nch,Nch,Nwindow,Ntask);
0192 
0193 % Debug code
0194 %S2.a = a_inv(:,1);
0195 %S2.G = G;
0196 %S2.Covs = Covs;
0197 %save('tmp2.mat','S2');
0198 %error('aho');
0199 
0200 if ~isempty(Model_ext), 
0201   %%%%%% Time window loop
0202   for j=1:Nwindow
0203     for n=1:Ntask
0204       Cov = Covs{n};     % Sg
0205     
0206       %%%% Calculate Inverse filter for current reconstruction
0207       A_z = a_inv(:,j,n);
0208       A_e = e_inv(:,j,n);
0209         
0210       % Inverse filter for each session
0211       % SB  = Ga * a_inv *Ga' + Cov
0212       SB      = G * vb_repmultiply(Gt, A_z) + Cov;  % Sigma_B
0213       SB      = SB + Ge * vb_repmultiply(Get, A_e);      
0214       SB_inv  = pinv( SB );        % Sigam_B^{-1}
0215       SB_cov(:,:,j,n) = SB_inv;
0216 
0217       Hj(j,n)  = - 0.5*( vb_log_det(SB) + Nch*log(sx(j,n)) );
0218         
0219       %%%% Calculate Inverse filter for current reconstruction
0220       GSB     = Gt(jx_act,:)*SB_inv;        % Ga' * Sigma_B^{-1}
0221       KW(:,:,j,n)   = vb_repmultiply(GSB, A_z(jx_act));
0222       GSB_ext = Get*SB_inv; 
0223       KW_ext(:,:,j,n) = vb_repmultiply(GSB_ext, A_e);
0224       %        SCVS(:,:,n) = SB_inv*Cov*SB_inv;
0225     end
0226   end % Timewindow loop
0227 else
0228   for j=1:Nwindow
0229     for n=1:Ntask
0230       Cov = Covs{n};     % Sg
0231     
0232       %%%% Calculate Inverse filter for current reconstruction
0233       A_z = a_inv(:,j,n);
0234         
0235       % Inverse filter for each session
0236       % SB  = Ga * a_inv *Ga' + Cov
0237       SB      = G * vb_repmultiply(Gt, A_z) + Cov;  % Sigma_B
0238       SB_inv  = inv( SB );        % Sigam_B^{-1}
0239       SB_cov(:,:,j,n) = SB_inv;
0240         
0241       Hj(j,n)  = - 0.5*( vb_log_det(SB) + Nch*log(sx(j,n)) );
0242         
0243       %%%% Calculate Inverse filter for current reconstruction
0244       GSB     = Gt(jx_act,:)*SB_inv;        % Ga' * Sigma_B^{-1}
0245       KW(:,:,j,n)   = vb_repmultiply(GSB, A_z(jx_act));
0246       %        SCVS(:,:,n) = SB_inv*Cov*SB_inv;
0247     end
0248   end % Timewindow loop
0249   
0250   KW_ext = [];
0251 end
0252 
0253 %VBfilt.a  = a_inv(jx_act,:,:);
0254 %VBfilt.sx = sx;
0255 VBfilt.sx = sx*(Model.bsnorm_meg^2);
0256 VBfilt.Cov = Covs;
0257 
0258 %VBfilt.Hj = Hj;
0259 Nmeg = size(Gact_b{1},1);
0260 Neeg = size(Gact_e{1},1);
0261 ix_meg = 1:Nmeg;
0262 ix_eeg = Nmeg+1:Nmeg+Neeg;
0263 KW(:,ix_meg) = KW(:,ix_meg)./Model.bsnorm_meg;
0264 KW(:,ix_eeg) = KW(:,ix_eeg)./Model.bsnorm_eeg;
0265 
0266 VBfilt.KW = KW;
0267 
0268 if ~isempty(KW_ext)
0269     KW_ext(:,ix_meg) = KW_ext(:,ix_meg)./Model.bsnorm_meg;
0270     KW_ext(:,ix_eeg) = KW_ext(:,ix_eeg)./Model.bsnorm_eeg;    
0271 end
0272 VBfilt.KW_ext = KW_ext;
0273 
0274 %VBfilt.SB_cov = SB_cov;
0275 VBfilt.ix_act = ix_act;
0276 
0277 VBfilt.G = Gact0{1};
0278 
0279 VBfilt.ix_meg = ix_meg;
0280 VBfilt.ix_eeg = ix_eeg;
0281 
0282 %if isfield(fileinfo,'channel_id')
0283 %    VBfilt.channel_id = fileinfo.channel_id;
0284 %end
0285 
0286 % Current Info
0287 Jinfo.version   = vbmeg('version');
0288 Jinfo.curr_type = 'Z-current';
0289 
0290 Jinfo.Lact = Lact;
0291 Jinfo.Wact = Wact;
0292 Jinfo.ix_act = ix_act;
0293 Jinfo.ix_act_ex = ix_act_ex;
0294 
0295 % ix_act : Vertex index corresponding to active current Zact
0296 % ix_act_ex : Vertex index corresponding to active current Jact
0297 % Wact   : Spatial smoothing matrix of focal window
0298 % Jact   = Wact * Zact
0299 
0300 return
0301 
0302 %%%% ---------------
0303 function    [bayes_parm, ix_area] = ...
0304             check_arg(bayes_parm,filt_parm)
0305 
0306 if isfield(filt_parm,'basisfile_meg'), 
0307     bayes_parm.basisfile_meg = filt_parm.basisfile_meg;
0308 end;
0309 if isfield(filt_parm,'basisfile_eeg'), 
0310     bayes_parm.basisfile_eeg = filt_parm.basisfile_eeg;
0311 end;
0312 if ~isfield(filt_parm,'ix_area'),  
0313     ix_area = []; 
0314 else
0315     ix_area = filt_parm.ix_area; 
0316 end;
0317 if isfield(filt_parm,'extra'),
0318   if isfield(filt_parm.extra,'basisfile'), 
0319     bayes_parm.extra.basisfile = filt_parm.extra.basisfile;
0320   end
0321 end
0322 
0323 
0324 return
0325 %%% --- Following fields are used for current calculation
0326 %if isfield(filt_parm,'megfile'),
0327 %    bayes_parm.megfile   = filt_parm.megfile  ;
0328 %end;
0329 %if isfield(filt_parm,'twin_meg'),
0330 %    bayes_parm.twin_meg  = filt_parm.twin_meg ;
0331 %end;
0332 %if isfield(filt_parm,'Tperiod'),
0333 %    bayes_parm.Tperiod   = filt_parm.Tperiod  ;
0334 %end;
0335 %if isfield(filt_parm,'Tnext'),
0336 %    bayes_parm.Tnext     = filt_parm.Tnext    ;
0337 %end;
0338 %
0339 %if ~isfield(filt_parm,'trial_average'),
0340 %    trial_average = ON;
0341 %else
0342 %    trial_average = filt_parm.trial_average;
0343 %end;
0344 %
0345 %bayes_parm.trial_average = trial_average;
0346 %
0347 %if ~isfield(filt_parm,'tsubsamp'),
0348 %    tsubsamp = [];
0349 %else
0350 %    tsubsamp = filt_parm.tsubsamp;
0351 %end;
0352 %if ~isfield(filt_parm,'overlap_mode'),
0353 %    overlap_mode = 0;
0354 %else
0355 %    overlap_mode = filt_parm.overlap_mode;
0356 %end;
vbmeg_current_filter_z_meeg

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
