Home > vbmeg > functions > estimation > bayes > vb_calc_current.m

vb_calc_current

PURPOSE ^

Current reconstruction using Bayesian inverse filter

SYNOPSIS ^

function [Z,W,Jinfo] = vb_calc_current(varargin)

DESCRIPTION ^

 Current reconstruction using Bayesian inverse filter
 (VBMEG public function)

 [syntax]
 [Z,W,Jinfo] = vb_calc_current(VBfilt,bayes_parm,curr_parm)
 [Z,W,Jinfo] = vb_calc_current(VBfilt,bexp,curr_parm)
 [Z,W,Jinfo] = vb_calc_current(proj_root,VBfilt,bayes_parm,curr_parm)
 [Z,W,Jinfo] = vb_calc_current(proj_root,VBfilt,bexp,curr_parm)

 [input]
 proj_root : <optional> <<string>> VBMEG project root directory. 
 VBfilt    : <<struct>> Inverse filter and related variables.
 bayes_parm: <<struct>> Parameters for variational Bayes method. Fields
             of this parameter is the same with that of job_vb.m. 
 bexp      : <<cell>> Cell of MEG data. If 'bexp' is given instead of
             struct 'bayes_parm', the following variables, required for
             current reconstruction will be set to ones in struct
             'VBfilt' (see example): 
                 Ntrial
                 Tsample
                 Twindow
 curr_parm : <<struct>> Current reconstruction parameters. 
 --- fields of curr_parm
  overlap_mode : <<int>> Specify how cortical current is estimated. 
   = 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
  ix_area      : <optional> <<vector>> Vertex indices to calculate
                 cortical current. If 'ix_area' is empty or not given,
                 currents in the active region are calculated. 
  trial_average: <<bool>> If true, average current over all sessions and
                 trials is calculated. 
  ix_trial     : <<vector>> Trial indices for which currents are
                 reconstructed.
  tsubsamp     : <optional> <<bool>> Specify subsampled time index
                 If 'tsubsamp' is empty or not given, time subsampling
                 is not done. 
 ---

 [output]
 Z    : <<matrix>> Z-current ([I,Nsample,Ntrial]). 
 W    : <<matrix>> Spatial smoothing matrix. 
 Jinfo: <<struct>> Information of estimated current. 
 --- fields of Jinfo
  ix_act   : <<vector>> Cortical vertex indices of Z-current. 
  ix_act_ex: <<vector>> Cortical vertex indices of J-current.
 ---

 [example]
 >> bexp{1} = vb_load_meg_data(meg_file); % assume single session
 >> load(filter_file,'VBfilt','bayes_parm'); % .vbfilt.mat
 >> vb_struct2vars(bayes_parm,{'twin_meg','Tperiod','Tnext'});
 >> bexp{1} = bexp{1}(:,twin_meg(1):twin_meg(2),:);
 >> VBfilt.Ntrial(1) = size(bexp{1},3);
 >> VBfilt.Tsample = twin_meg(2)-twin_meg(1)+1;
 >> VBfilt.Twindow = vb_calc_timewindow(VBfilt.Tsample,Tperiod,Tnext);
 >> curr_parm.overlap_mode  = 0;
 >> curr_parm.ix_area       = [];
 >> curr_parm.trial_average = true;
 >> curr_parm.ix_trial      = [];
 >> curr_parm.tsubsamp      = [];
 >> [Z,W,Jinfo] = vb_calc_current(proj_root,VBfilt,bexp,curr_parm);

 [history]
 2010-06-28 Taku Yoshioka
 2010-07-12 Taku Yoshioka
  'trial_average' support.
 2010-09-29 Taku Yoshioka
  New usage of function implemented ('bexp' given instead of
  'bayes_parm'). 

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Z,W,Jinfo] = vb_calc_current(varargin)
0002 % Current reconstruction using Bayesian inverse filter
0003 % (VBMEG public function)
0004 %
0005 % [syntax]
0006 % [Z,W,Jinfo] = vb_calc_current(VBfilt,bayes_parm,curr_parm)
0007 % [Z,W,Jinfo] = vb_calc_current(VBfilt,bexp,curr_parm)
0008 % [Z,W,Jinfo] = vb_calc_current(proj_root,VBfilt,bayes_parm,curr_parm)
0009 % [Z,W,Jinfo] = vb_calc_current(proj_root,VBfilt,bexp,curr_parm)
0010 %
0011 % [input]
0012 % proj_root : <optional> <<string>> VBMEG project root directory.
0013 % VBfilt    : <<struct>> Inverse filter and related variables.
0014 % bayes_parm: <<struct>> Parameters for variational Bayes method. Fields
0015 %             of this parameter is the same with that of job_vb.m.
0016 % bexp      : <<cell>> Cell of MEG data. If 'bexp' is given instead of
0017 %             struct 'bayes_parm', the following variables, required for
0018 %             current reconstruction will be set to ones in struct
0019 %             'VBfilt' (see example):
0020 %                 Ntrial
0021 %                 Tsample
0022 %                 Twindow
0023 % curr_parm : <<struct>> Current reconstruction parameters.
0024 % --- fields of curr_parm
0025 %  overlap_mode : <<int>> Specify how cortical current is estimated.
0026 %   = 0 : current is averaged over overlapped time window.
0027 %   = 1 : current is not averaged for overlapped window.
0028 %         current time series of each time windows
0029 %         are concatenated sequentially for spectral analysis
0030 %  ix_area      : <optional> <<vector>> Vertex indices to calculate
0031 %                 cortical current. If 'ix_area' is empty or not given,
0032 %                 currents in the active region are calculated.
0033 %  trial_average: <<bool>> If true, average current over all sessions and
0034 %                 trials is calculated.
0035 %  ix_trial     : <<vector>> Trial indices for which currents are
0036 %                 reconstructed.
0037 %  tsubsamp     : <optional> <<bool>> Specify subsampled time index
0038 %                 If 'tsubsamp' is empty or not given, time subsampling
0039 %                 is not done.
0040 % ---
0041 %
0042 % [output]
0043 % Z    : <<matrix>> Z-current ([I,Nsample,Ntrial]).
0044 % W    : <<matrix>> Spatial smoothing matrix.
0045 % Jinfo: <<struct>> Information of estimated current.
0046 % --- fields of Jinfo
0047 %  ix_act   : <<vector>> Cortical vertex indices of Z-current.
0048 %  ix_act_ex: <<vector>> Cortical vertex indices of J-current.
0049 % ---
0050 %
0051 % [example]
0052 % >> bexp{1} = vb_load_meg_data(meg_file); % assume single session
0053 % >> load(filter_file,'VBfilt','bayes_parm'); % .vbfilt.mat
0054 % >> vb_struct2vars(bayes_parm,{'twin_meg','Tperiod','Tnext'});
0055 % >> bexp{1} = bexp{1}(:,twin_meg(1):twin_meg(2),:);
0056 % >> VBfilt.Ntrial(1) = size(bexp{1},3);
0057 % >> VBfilt.Tsample = twin_meg(2)-twin_meg(1)+1;
0058 % >> VBfilt.Twindow = vb_calc_timewindow(VBfilt.Tsample,Tperiod,Tnext);
0059 % >> curr_parm.overlap_mode  = 0;
0060 % >> curr_parm.ix_area       = [];
0061 % >> curr_parm.trial_average = true;
0062 % >> curr_parm.ix_trial      = [];
0063 % >> curr_parm.tsubsamp      = [];
0064 % >> [Z,W,Jinfo] = vb_calc_current(proj_root,VBfilt,bexp,curr_parm);
0065 %
0066 % [history]
0067 % 2010-06-28 Taku Yoshioka
0068 % 2010-07-12 Taku Yoshioka
0069 %  'trial_average' support.
0070 % 2010-09-29 Taku Yoshioka
0071 %  New usage of function implemented ('bexp' given instead of
0072 %  'bayes_parm').
0073 %
0074 % Copyright (C) 2011, ATR All Rights Reserved.
0075 % License : New BSD License(see VBMEG_LICENSE.txt)
0076 
0077 %
0078 % Verbose level
0079 %
0080 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0081 global vbmeg_inst;
0082 const = vb_define_verbose;
0083 if isfield(vbmeg_inst,'verbose_level'), 
0084   verbose_level = vbmeg_inst.verbose_level;
0085 else
0086   verbose_level = const.VERBOSE_LEVEL_NOTICE;
0087 end;
0088 
0089 %
0090 % Check input arguments
0091 %
0092 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0093 if length(varargin) == 3
0094   proj_root = [];
0095   VBfilt = varargin{1};
0096   curr_parm = varargin{3};
0097   if isstruct(varargin{2}), 
0098     bayes_parm = varargin{2};
0099   else
0100     bexp = varargin{2}; % 'bexp' is a cell array.
0101   end;
0102 elseif length(varargin) == 4
0103   proj_root = varargin{1};
0104   VBfilt = varargin{2};
0105   if isstruct(varargin{3}), 
0106     bayes_parm = varargin{3};
0107   else
0108     bexp = varargin{3}; % 'bexp' is a cell array.
0109   end;
0110   curr_parm = varargin{4};
0111 else
0112   error('Error: invalid usage of vb_calc_current.');
0113 end;
0114 
0115 vb_struct2vars(curr_parm,{'overlap_mode','ix_area','trial_average', ...
0116                     'ix_trial','tsubsamp'});
0117 
0118 if exist('bayes_parm'), 
0119   if ~isempty(proj_root);
0120     bayes_parm_abs = vb_parm_absolute_path(proj_root,bayes_parm);
0121   else
0122     bayes_parm_abs = bayes_parm;
0123   end;
0124 end;
0125 
0126 %
0127 % MEG data preparation
0128 %
0129 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0130 if exist('bayes_parm_abs'), 
0131   [bexp,Ntrial,Nch,Tsample,Twindow] ...
0132       = vb_megdata_preparation(bayes_parm_abs);
0133   Nsession = length(bexp);    % Number of sessions
0134 else
0135   Tsample = VBfilt.Tsample;
0136   Twindow = VBfilt.Twindow;
0137   Ntrial = VBfilt.Ntrial;
0138   Nsession = length(bexp);    % Number of sessions
0139 end;
0140 
0141 %
0142 % Area index in which current is calculated
0143 %
0144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0145 ix_act = VBfilt.ix_act;
0146 ix_act_ex = VBfilt.ix_act_ex;
0147 
0148 if ~isempty(ix_area),
0149   % Select vertex index 'ix_area' within the active current region
0150   [jx_area_ex, ix_area_ex] = vb_index2indexsub(ix_area, ix_act_ex);
0151 else
0152   jx_area_ex = 1:length(ix_act_ex);
0153 end
0154 
0155 W = VBfilt.W(jx_area_ex,:);
0156 jx_act = find(sum(W,1)>0);
0157 W = W(:,jx_act);
0158 
0159 ix_act = ix_act(jx_act); % active index of Z-current
0160 ix_act_ex = ix_act_ex(jx_area_ex); % active index of J-current
0161 
0162 % # of active vertex
0163 Njact_area = length(jx_act);
0164 
0165 % # of extra dipoles
0166 if isempty(VBfilt.KW_ext), Njext = 0;
0167 else Njext = size(VBfilt.KW_ext{1},1); end;
0168 
0169 Jinfo.ix_act = ix_act;
0170 Jinfo.ix_act_ex = ix_act_ex;
0171 
0172 %
0173 % Temporal smoothing window weight
0174 %
0175 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0176 if isempty(tsubsamp),
0177   tsubsamp = 1:Tsample;
0178 end;
0179 
0180 [Tweight ,Tindex, Nindex, Nsample] = ...
0181     vb_overlapped_timewindow_weight(Twindow, Tsample, tsubsamp, overlap_mode);
0182 
0183 Nwindow = length(Nindex);       % # of time window
0184 vb_disp(['Number of time windows : ' num2str(Nwindow)]);
0185 
0186 if overlap_mode == 1,
0187   vb_disp('Non-overlapped concatenation mode'); 
0188 end;
0189 
0190 %
0191 % Calculate Z-current
0192 %
0193 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0194 Ntrial_all = 0;
0195 
0196 % Session loop
0197 for n=1:Nsession
0198   % Trial average
0199   if trial_average, 
0200     bexp{n} = mean(bexp{n},3); 
0201     ix_trial(n) = 1;
0202   elseif isempty(curr_parm.ix_trial), 
0203     ix_trial(n) = VBfilt.Ntrial(n);
0204   end;
0205   
0206   % Number of trials
0207   Ntrial(n) = length(ix_trial(n));
0208 
0209   % Get inverse filter for all time window
0210   for j=1:Nwindow
0211     Nid = Nindex{j};    % index in the total subsampled data
0212     %KW{j} = VBfilt.KW(:,:,j,n);
0213     KW{j} = VBfilt.KW(jx_act,:,j,n);
0214   end;
0215 
0216   % Trial loop
0217   for m=ix_trial(n)
0218     Zact = zeros(Njact_area,Nsample);
0219         
0220     for j=1:Nwindow
0221       % Subsampling time index
0222       Tid = Tindex{j};    % subsampled time index
0223       Nid = Nindex{j};    % index in the total subsampled data
0224 
0225       if isempty(Nid), continue; end;
0226             
0227       % Time window smoothing
0228       Bt = vb_repmultiply(bexp{n}(:,Tid,m), Tweight{j});
0229       Zact(:,Nid) = Zact(:,Nid) + (KW{j} * Bt);
0230     end; % Timewindow loop END
0231     
0232     Ntrial_all = Ntrial_all+1;
0233     Z(:,:,Ntrial_all) = Zact;
0234   end; % Trial loop END
0235   
0236   vb_disp_nonl('.');
0237 end; % Session loop END
0238 vb_disp_nonl(sprintf('\n'));
0239 
0240 Jinfo.Ntrial   = Ntrial;
0241 Jinfo.Nsession = Nsession;
0242 
0243 %
0244 % Calculate extra-dipole current (not implemented: 2010-06-28)
0245 %
0246 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0247 
0248 return;

Generated on Mon 22-May-2023 06:53:56 by m2html © 2005