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)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;