vbmeg_current_calc_z

 Current reconstruction using Bayesian inverse filter

 --- Syntax
  [Zact ,Jinfo, bayes_parm, vb_parm, MEGinfo] ...
          = vbmeg_current_calc_z(proj_root, curr_parm);
 --- Input
 curr_parm.resultfile : result file obtained by VB estimation

 --- Optional Input
 curr_parm.trial_average = [ON] : average current over all sessions 
               = OFF  : current for each session

 curr_parm.ix_area : Position index to calculate estimated current
           If 'ix_area' is empty or not given, 
              currents in the active region are calculated
 curr_parm.tsubsamp   : Specify subsampled time index
           If 'tsubsamp' is empty or not given, 
              time subsampling is not done
 curr_parm.overlap_mode 
   = 0 : current is averaged over overlapped time window
   = 1 : current is not averaged for overlapped window
         current time series of each time windows 
         are concatenated sequentially for spectral analysis
 --- If following field is given, these values are used 
     instead of bayes_parm field in result file
 curr_parm.basisfile
 curr_parm.megfile  
 curr_parm.twin_meg 
 curr_parm.Tperiod  
 curr_parm.Tnext    
 ---Output
 Zact    : active current

 Zact(n,t,:) is the current at the vertex 'ix_act(n)' & the time 't'
 Zact(Nact,Nsample)          for trial_average = ON 
 Zact(Nact,Nsample,Ntrials)  for trial_average = OFF
   Nact     : # of active region, 
   Nsample  : # of time sample, 
   Ntrial  : # of trials in all session]
 Jinfo.ix_act  : Vertex index corresponding to active current

 * Inverse filter calculation is done in vb_invfilter_calc
 2006-09-03 M. Sato
 * Non-overlapped concatenation mode is added for spectral analysis
 2008-10-29 M. Sato
 * Inverse filter calculation is done in vb_invfilter_calc
 VBfilt.sx = Observation noise variance (Nwindow,Ntask)
 VBfilt.Hj = Model entropy (Nwindow,Ntask)
 VBfilt.KW = VB inverse filter  (Njact,Nch,Nwindow,Ntask)
 VBfilt.SB_cov = posterior sensor covariance matrix (Nch,Nch,Nwindow,Ntask)

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

0001 function [Zact_ave ,Jinfo, bayes_parm, vb_parm, MEGinfo] ...
0002           = vbmeg_current_calc_z(proj_root, curr_parm);
0003 % Current reconstruction using Bayesian inverse filter
0004 %
0005 % --- Syntax
0006 %  [Zact ,Jinfo, bayes_parm, vb_parm, MEGinfo] ...
0007 %          = vbmeg_current_calc_z(proj_root, curr_parm);
0008 % --- Input
0009 % curr_parm.resultfile : result file obtained by VB estimation
0010 %
0011 % --- Optional Input
0012 % curr_parm.trial_average = [ON] : average current over all sessions
0013 %               = OFF  : current for each session
0014 %
0015 % curr_parm.ix_area : Position index to calculate estimated current
0016 %           If 'ix_area' is empty or not given,
0017 %              currents in the active region are calculated
0018 % curr_parm.tsubsamp   : Specify subsampled time index
0019 %           If 'tsubsamp' is empty or not given,
0020 %              time subsampling is not done
0021 % curr_parm.overlap_mode
0022 %   = 0 : current is averaged over overlapped time window
0023 %   = 1 : current is not averaged for overlapped window
0024 %         current time series of each time windows
0025 %         are concatenated sequentially for spectral analysis
0026 % --- If following field is given, these values are used
0027 %     instead of bayes_parm field in result file
0028 % curr_parm.basisfile
0029 % curr_parm.megfile
0030 % curr_parm.twin_meg
0031 % curr_parm.Tperiod
0032 % curr_parm.Tnext
0033 % ---Output
0034 % Zact    : active current
0035 %
0036 % Zact(n,t,:) is the current at the vertex 'ix_act(n)' & the time 't'
0037 % Zact(Nact,Nsample)          for trial_average = ON
0038 % Zact(Nact,Nsample,Ntrials)  for trial_average = OFF
0039 %   Nact     : # of active region,
0040 %   Nsample  : # of time sample,
0041 %   Ntrial  : # of trials in all session]
0042 % Jinfo.ix_act  : Vertex index corresponding to active current
0043 %
0044 % * Inverse filter calculation is done in vb_invfilter_calc
0045 % 2006-09-03 M. Sato
0046 % * Non-overlapped concatenation mode is added for spectral analysis
0047 % 2008-10-29 M. Sato
0048 % * Inverse filter calculation is done in vb_invfilter_calc
0049 % VBfilt.sx = Observation noise variance (Nwindow,Ntask)
0050 % VBfilt.Hj = Model entropy (Nwindow,Ntask)
0051 % VBfilt.KW = VB inverse filter  (Njact,Nch,Nwindow,Ntask)
0052 % VBfilt.SB_cov = posterior sensor covariance matrix (Nch,Nch,Nwindow,Ntask)
0053 %
0054 % Copyright (C) 2011, ATR All Rights Reserved.
0055 % License : New BSD License(see VBMEG_LICENSE.txt)
0056 
0057 MinExp = -50;
0058 
0059 % multiple bayse file for multiple conditions
0060 
0061 if ~isempty(proj_root)
0062     filterfile = [proj_root filesep curr_parm.filterfile];
0063 else
0064     filterfile = curr_parm.filterfile;
0065 end  
0066 
0067 %%%%%% load VBMEG estimated result
0068 load(filterfile, 'VBfilt','Jinfo','bayes_parm','vb_parm');
0069 
0070 %%%%%% check parameter of 'curr_parm'
0071 % Values of 'curr_parm' fields dominates over
0072 %   those of 'bayes_parm' in bayesfile
0073 [bayes_parm, ix_area, trial_average, tsubsamp, overlap_mode, verbose] = ...
0074             check_arg(bayes_parm, curr_parm);
0075             
0076 if ~isempty(proj_root);
0077   bayes_parm_abs = vb_parm_absolute_path(proj_root, bayes_parm);
0078 else
0079   bayes_parm_abs = bayes_parm;
0080 end
0081 
0082 % current file directory
0083 if ~isempty(proj_root)
0084     currfile = fullfile(proj_root, curr_parm.currfile);
0085 else
0086     currfile = [curr_parm.currfile];
0087 end  
0088 
0089 curr_root  = fileparts(currfile);
0090 
0091 % jactdir is relative path from current file
0092 if isfield(curr_parm,'jactdir') && ~isempty(curr_parm.jactdir)
0093     jactdir    = curr_parm.jactdir;
0094     jactdir_ab = fullfile(curr_root, jactdir);
0095     
0096     %check directory
0097     if exist(jactdir_ab,'dir') ~= 0
0098         msg = [ 'jact_dir : ' jactdir_ab ' exists. Do you over write?'];
0099         btn = questdlg(msg, 'confirm', 'Yes', 'No', 'Yes');
0100         if strcmp(btn, 'No')
0101             error('processing aborted.');
0102             return;
0103         end
0104     else
0105         res = mkdir(curr_root, jactdir);
0106         if  res~= 1,  error('mkdir failed.'); end
0107     end
0108 else
0109     jactdir_ab = [];
0110 end
0111 
0112 %%%%%% MEG data preparation
0113 % B      : MEG data
0114 [B, Ntrial, Nch, Tsample, Twindow] = ...
0115     vb_megdata_preparation(bayes_parm_abs);
0116 
0117 MEGinfo = vb_load_meg_info(bayes_parm_abs.megfile{1});
0118 
0119 Nsession = length(B);     % Number of files
0120 
0121 fprintf('Number of session: %d\n', Nsession);
0122 
0123 %% Constant
0124 
0125 % Temporal subsampling index
0126 if isempty(tsubsamp),
0127     tsubsamp = 1:Tsample;
0128 end
0129 
0130 %%%%%% Temporal smoothing window weight
0131 [Tweight ,Tindex, Nindex, Nsample] = ...
0132     vb_overlapped_timewindow_weight(Twindow, Tsample, tsubsamp, overlap_mode);
0133 
0134 Nwindow = size(Twindow,1);       % # of time window
0135 
0136 sx = VBfilt.sx;
0137 Hj = VBfilt.Hj ;
0138 KW = VBfilt.KW ;
0139 SB_cov = VBfilt.SB_cov ;
0140 
0141 % # of active vertex
0142 Njact = size(KW,1);
0143 Ntask = size(KW,4);
0144 Nch = size(SB_cov,1);
0145 
0146 fprintf('Active vertex number = %d\n',Njact)
0147 fprintf('Active sensor number = %d\n',Nch)
0148 
0149 if overlap_mode == 1,fprintf('Non-overlapped concatenation mode\n'); end;
0150 
0151 %%%%%% Initialization
0152 Zact = zeros(Njact,Nsample);
0153 Zact_ave = zeros(Njact,Nsample);
0154 LP   = zeros(Ntask,1);
0155 
0156 % Current Info
0157 Jinfo.Tsample = Tsample;
0158 Jinfo.tindex  = tsubsamp;
0159 Jinfo.Nwindow = Nwindow;
0160 Jinfo.Ntask   = Ntask;
0161 
0162 % MEG time window index
0163 Tstart  = bayes_parm.twin_meg(1);
0164 Tend    = bayes_parm.twin_meg(2);
0165 if isempty(tsubsamp)
0166     Jinfo.Tsample = Tstart:Tend;
0167 else
0168     Jinfo.Tsample = tsubsamp + Tstart - 1;
0169 end
0170 
0171 %%%%%% Time window loop
0172 % MEG for each trial
0173 filename = [];
0174 
0175 Ntrial_all = 0;
0176 %nall = 0;
0177 
0178 for n=1:Nsession
0179     Ntry = Ntrial(n);
0180     
0181     if Ntry == 0, continue; end;
0182 
0183     Bs = B{n};
0184     
0185     for m=1:Ntry
0186         
0187         Zact = zeros(Njact,Nsample);
0188         Post = zeros(Ntask,Nwindow);
0189         
0190         for j=1:Nwindow
0191             % Subsampling time index
0192             Tid = Tindex{j};    % subsampled time index
0193             Nid = Nindex{j};    % index in the total subsampled data
0194             Nt  = length(Tid);
0195             
0196             if Nt == 0, continue; end;
0197     
0198             %%%% Current reconstruction
0199             [Zp ,post]= vb_current_calc(Bs(:,Tid,m), ...
0200                         KW(:,:,j,:),SB_cov(:,:,j,:),Hj(j,:),sx(j,:));
0201             
0202             %%%% Time window smoothing
0203             weight = Tweight{j};
0204             Zact(:,Nid) = Zact(:,Nid) ...
0205                         + vb_repmultiply(Zp , weight);
0206             Post(:,j)  = post;
0207             
0208         end % Timewindow loop
0209         
0210         % Trial average current
0211         Zact_ave   = Zact_ave + Zact;
0212         Ntrial_all = Ntrial_all + 1;
0213         
0214         if ~isempty(jactdir_ab)
0215             fname = sprintf('data_s%04dt%04d',n,m);
0216             filename{m} = fname;
0217             vb_fsave([jactdir_ab '/' fname], 'Zact','Jinfo','Post');
0218             if verbose==1,
0219                 fprintf('.')
0220                 if rem(m, 20) == 0 % linefeed per 20
0221                     fprintf('\n');
0222                 end
0223             end
0224         end
0225     end; % Trial loop
0226 end
0227 
0228 Zact_ave = Zact_ave/Ntrial_all;
0229 
0230 Jinfo.Ntrial   = Ntrial;
0231 Jinfo.Nsession = Nsession;
0232 Jinfo.jactdir  = jactdir;
0233 Jinfo.filename = filename;
0234 
0235 if verbose==1,
0236     fprintf('\n')
0237 end
0238 
0239 % ix_act : Vertex index corresponding to active current Zact
0240 % ix_act_ex : Vertex index corresponding to active current Jact
0241 % Wact   : Spatial smoothing matrix of focal window
0242 % Jact   = Wact * Zact
0243 
0244 return
0245 
0246 %%%% ---------------
0247 function [bayes_parm,ix_area,trial_average,tsubsamp,overlap_mode,verbose]= ...
0248             check_arg(bayes_parm,curr_parm)
0249 
0250 if isfield(curr_parm,'basisfile'), 
0251     bayes_parm.basisfile = curr_parm.basisfile;
0252 end;
0253 if isfield(curr_parm,'megfile'), 
0254     bayes_parm.megfile   = curr_parm.megfile  ;
0255 end;
0256 if isfield(curr_parm,'twin_meg'), 
0257     bayes_parm.twin_meg  = curr_parm.twin_meg ;
0258 end;
0259 if isfield(curr_parm,'Tperiod'), 
0260     bayes_parm.Tperiod   = curr_parm.Tperiod  ;
0261 end;
0262 if isfield(curr_parm,'Tnext'), 
0263     bayes_parm.Tnext     = curr_parm.Tnext    ;
0264 end;
0265 
0266 if ~isfield(curr_parm,'trial_average'), 
0267     trial_average = ON; 
0268 else
0269     trial_average = curr_parm.trial_average; 
0270 end;
0271 
0272 bayes_parm.trial_average = trial_average;
0273 
0274 if ~isfield(curr_parm,'ix_area'),  
0275     ix_area = []; 
0276 else
0277     ix_area = curr_parm.ix_area; 
0278 end;
0279 if ~isfield(curr_parm,'tsubsamp'), 
0280     tsubsamp = []; 
0281 else
0282     tsubsamp = curr_parm.tsubsamp; 
0283 end;
0284 if ~isfield(curr_parm,'overlap_mode'),     
0285     overlap_mode = 0; 
0286 else
0287     overlap_mode = curr_parm.overlap_mode; 
0288 end;
0289 if ~isfield(curr_parm,'verbose'),     
0290     verbose = 1; 
0291 else
0292     verbose = curr_parm.verbose; 
0293 end;