Current reconstruction using Bayesian inverse filter
--- Syntax
[Jact ,ix_act, Jbck, ix_bck, Var] =...
vbmeg_current_reconstruct(proj_root, bayesfile)
[Jact ,ix_act, Jbck, ix_bck, Var] =...
vbmeg_current_reconstruct(proj_root, bayesfile, ix_area)
[Jact ,ix_act, Jbck, ix_bck, Var] =...
vbmeg_current_reconstruct(proj_root, bayesfile, ix_area, trial_average)
[Jact ,ix_act, Jbck, ix_bck, Var] =...
vbmeg_current_reconstruct(proj_root, bayesfile, ix_area, ...
trial_average, tsubsamp)
--- Input
bayesfile : result file obtained by VB estimation
MEG data : B{Nsession}(Nsensors, Tsample, Ntrials)
--- Optional Input
trial_average = [ON] : average current over all sessions
= OFF : current for each session
ix_area : Position index to calculate estimated current
If 'ix_area' is empty or not given,
currents in the active region are calculated
tsubsamp : Specify subsampled time index
If 'tsubsamp' is empty or not given,
time subsampling is not done
---Output
Jact : active current
ix_act : Vertex index corresponding to active current
Jact(n,t,:) is the current at the vertex 'ix_act(n)' & the time 't'
Jact(Nact,Nsample) for trial_average = ON
Jact(Nact,Nsample,Ntrials) for trial_average = OFF
Nact : # of active region,
Nsample : # of time sample,
Ntrials : # of trials in all session]
Jbck : background current
ix_bck : Vertex index corresponding to background current
same for 'Jact'
Var : posterior variance of act_area
When active area is empty,
posterior variance for background area is calculated
2004-9-25 M. Sato
2004-12-5 M. Sato : Modified
2004-12-22 M. Sato : Modified
2005-1-5 M. Sato Modified
2005-1-15 M. Sato Modified
2005-08-19 O. Yamashita
* Posterior variance is calculated simultaneously
* Variables are renamed
* Corresponding to new bayesfile format
2005-08-19 O. Yamashita ver.30b
2005-09-30 O. Yamashita
* Minor bug fix when Lact >= 2 & Lact_var = 1
2006-08-31 M. Sato
* Remove focal area from global area
* VBMEG result & liedfield is loaded by vb_load_result
* Inverse filter calculation is done in vb_invfilter_calc
2006-09-03 M. Sato
* Non-overlapped concatenation mode is added for spectral analysis
2009-04-02 Taku Yoshioka
Parameter name changed within this code for readability
(just replacing 'resultfile' to bayesfile)
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [Jact ,ix_act, Jbck, ix_bck, Var] ... 0002 = vbmeg_current_reconstruct(proj_root, bayesfile, ... 0003 ix_area, trial_average, tsubsamp, mode) 0004 % Current reconstruction using Bayesian inverse filter 0005 % 0006 % --- Syntax 0007 % [Jact ,ix_act, Jbck, ix_bck, Var] =... 0008 % vbmeg_current_reconstruct(proj_root, bayesfile) 0009 % [Jact ,ix_act, Jbck, ix_bck, Var] =... 0010 % vbmeg_current_reconstruct(proj_root, bayesfile, ix_area) 0011 % [Jact ,ix_act, Jbck, ix_bck, Var] =... 0012 % vbmeg_current_reconstruct(proj_root, bayesfile, ix_area, trial_average) 0013 % [Jact ,ix_act, Jbck, ix_bck, Var] =... 0014 % vbmeg_current_reconstruct(proj_root, bayesfile, ix_area, ... 0015 % trial_average, tsubsamp) 0016 % 0017 % --- Input 0018 % bayesfile : result file obtained by VB estimation 0019 % 0020 % MEG data : B{Nsession}(Nsensors, Tsample, Ntrials) 0021 % 0022 % --- Optional Input 0023 % trial_average = [ON] : average current over all sessions 0024 % = OFF : current for each session 0025 % 0026 % ix_area : Position index to calculate estimated current 0027 % If 'ix_area' is empty or not given, 0028 % currents in the active region are calculated 0029 % tsubsamp : Specify subsampled time index 0030 % If 'tsubsamp' is empty or not given, 0031 % time subsampling is not done 0032 % 0033 % ---Output 0034 % Jact : active current 0035 % ix_act : Vertex index corresponding to active current 0036 % 0037 % Jact(n,t,:) is the current at the vertex 'ix_act(n)' & the time 't' 0038 % Jact(Nact,Nsample) for trial_average = ON 0039 % Jact(Nact,Nsample,Ntrials) for trial_average = OFF 0040 % Nact : # of active region, 0041 % Nsample : # of time sample, 0042 % Ntrials : # of trials in all session] 0043 % 0044 % Jbck : background current 0045 % ix_bck : Vertex index corresponding to background current 0046 % same for 'Jact' 0047 % Var : posterior variance of act_area 0048 % When active area is empty, 0049 % posterior variance for background area is calculated 0050 % 0051 % 2004-9-25 M. Sato 0052 % 2004-12-5 M. Sato : Modified 0053 % 2004-12-22 M. Sato : Modified 0054 % 2005-1-5 M. Sato Modified 0055 % 2005-1-15 M. Sato Modified 0056 % 2005-08-19 O. Yamashita 0057 % * Posterior variance is calculated simultaneously 0058 % * Variables are renamed 0059 % * Corresponding to new bayesfile format 0060 % 2005-08-19 O. Yamashita ver.30b 0061 % 2005-09-30 O. Yamashita 0062 % * Minor bug fix when Lact >= 2 & Lact_var = 1 0063 % 2006-08-31 M. Sato 0064 % * Remove focal area from global area 0065 % * VBMEG result & liedfield is loaded by vb_load_result 0066 % * Inverse filter calculation is done in vb_invfilter_calc 0067 % 2006-09-03 M. Sato 0068 % * Non-overlapped concatenation mode is added for spectral analysis 0069 % 2009-04-02 Taku Yoshioka 0070 % Parameter name changed within this code for readability 0071 % (just replacing 'resultfile' to bayesfile) 0072 % 0073 % Copyright (C) 2011, ATR All Rights Reserved. 0074 % License : New BSD License(see VBMEG_LICENSE.txt) 0075 0076 0077 if ~exist('trial_average','var'), trial_average = ON; end; 0078 0079 if ~exist('ix_area','var'), ix_area = []; end; 0080 if ~exist('tsubsamp','var'), tsubsamp = []; end; 0081 if ~exist('mode','var'), mode = 0; end; 0082 0083 error('This function is a canditate to be removed, but invoked.'); 0084 0085 %%%%%% load result and MEG data 0086 [B, Model, nGact, nGall, nGGall, Cov, Wact, Wbck, ... 0087 bsnorm, ix_act, ix_bck, vb_parm] ... 0088 = vb_load_result(proj_root, bayesfile, ix_area); 0089 0090 % B : MEG data 0091 % Model : Estimated model 0092 % nGact : Normalized leadfield of focal window 0093 % nGall : Normalized leadfield of global window 0094 % Cov : Sensor noise covariance matrix 0095 % bsnorm : Normalizing constant 0096 % Wact : Spatial smoothing matrix of focal window 0097 % Wbck : Spatial smoothing matrix of global window 0098 % ix_act : Vertex index corresponding to active current 0099 % ix_bck : Vertex index corresponding to background current 0100 0101 % # of vertex 0102 Nvact = vb_parm.Nvact; % # of active vertex 0103 Njact_area = vb_parm.Njact_area; 0104 Njall_area = vb_parm.Njall_area; 0105 0106 %% Constant 0107 Nsession = vb_parm.Nsession; % Number of sessions 0108 Ntrials = zeros(1,Nsession); 0109 for n = 1:Nsession 0110 Ntrials(n) = size(B{n},3); 0111 end 0112 Ntotal = sum(Ntrials); % Total number of trials in all sessions 0113 0114 Tsample = vb_parm.Tsample; % Number of total time sample 0115 Twindow = vb_parm.Twindow; 0116 0117 % Temporal subsampling index 0118 if isempty(tsubsamp), 0119 tsubsamp = 1:Tsample; 0120 end 0121 0122 % Temporal smoothing window weight 0123 [Tweight ,Tindex, Nindex, Nsample] = ... 0124 vb_overlapped_timewindow_weight(Twindow, Tsample, tsubsamp, mode); 0125 0126 Nwindow = length(Nindex); % # of time window 0127 0128 % Nsample : # of subsampled data 0129 0130 if mode == 1,fprintf('Non-overlapped concatenation mode\n'); end; 0131 0132 %% Initialization 0133 if trial_average == ON 0134 % Current averaged over trials 0135 Jact = zeros(Njact_area,Nsample); 0136 Jbck = zeros(Njall_area,Nsample); 0137 else 0138 % Current for each trial 0139 Jact = zeros(Njact_area,Nsample,Ntotal); 0140 Jbck = zeros(Njall_area,Nsample,Ntotal); 0141 end 0142 0143 % Posterior variance 0144 if Njact_area > 0 0145 Var = zeros(Njact_area,Nsample); 0146 elseif Njall_area > 0 0147 Var = zeros(Njall_area,Nsample); 0148 end 0149 0150 %% Estimated current variance 0151 if isfield(Model,'a') 0152 a_inv = Model.a; 0153 else 0154 a_inv = zeros(Nvact,Nwindow); 0155 end 0156 if isfield(Model,'v') 0157 v_inv = Model.v; 0158 else 0159 v_inv = zeros(Nwindow,1); 0160 end 0161 0162 % Time window loop 0163 for j=1:Nwindow 0164 % Subsampling time index 0165 Tid = Tindex{j}; % subsampled time index 0166 Nid = Nindex{j}; % index in the total subsampled data 0167 Nt = length(Tid); 0168 0169 if Nt == 0, continue; end; 0170 0171 if Njact_area > 0 0172 Var_sum = zeros(Njact_area,1); 0173 elseif Njall_area > 0 0174 Var_sum = zeros(Njall_area,1); 0175 end 0176 0177 % Session loop 0178 Stry = 1; 0179 0180 for n=1:Nsession 0181 Ntry = Ntrials(n); 0182 0183 G = nGact{n}; % Ga 0184 GbW = nGall{n}; % Gb*Wb' 0185 GG = nGGall{n}; % Gb*Gb' 0186 Covs = Cov{n}; % Sg 0187 Nch = size(G,1); 0188 0189 %%%% Calculate Inverse filter for current reconstruction 0190 [KW, KW0, var_tmp] = vb_invfilter_calc(a_inv(:,j), v_inv(j), ... 0191 G, GbW, GG, Covs, Wact, Wbck); 0192 0193 %%%% Time window smoothing 0194 weight = repmat(Tweight{j}, [Nch 1]); 0195 0196 if trial_average == ON 0197 % MEG trial_average 0198 Bt = sum( B{n}(:,Tid,:), 3).* weight; 0199 % Current reconstruction 0200 Jact(:,Nid) = Jact(:,Nid) + (KW * Bt)/(bsnorm*Ntotal); 0201 Jbck(:,Nid) = Jbck(:,Nid) + (KW0 * Bt)/(bsnorm*Ntotal); 0202 else 0203 % MEG for each trial 0204 Bt = B{n}(:,Tid,:); 0205 for m=1:Ntry 0206 Bt(:,:,m) = Bt(:,:,m) .* weight; 0207 end 0208 Bt = reshape( Bt, [Nch, Nt*Ntry]); 0209 0210 % Trial index 0211 ixTry = [Stry:Stry+Ntry-1]; 0212 Stry = Stry+Ntry; 0213 0214 % Current reconstruction 0215 Jact(:,Nid,ixTry) = Jact(:,Nid,ixTry) ... 0216 + reshape( (KW * Bt)/bsnorm, [Njact_area, Nt, Ntry]); 0217 Jbck(:,Nid,ixTry) = Jbck(:,Nid,ixTry) ... 0218 + reshape( (KW0 * Bt)/bsnorm, [Njall_area, Nt, Ntry]); 0219 end 0220 0221 %%%% Posterior variance calculation 0222 0223 % trial average of inverse variance 0224 var_tmp = max(var_tmp,eps); 0225 Var_sum = Var_sum + (Ntry./var_tmp); 0226 0227 end; % Session loop END 0228 0229 % Var_sum = Var_sum/Ntotal; 0230 Var_sum = Model.sx(j)./(Var_sum * bsnorm^2); 0231 Var(:,Nid) = Var(:,Nid) + Var_sum * Tweight{j}; 0232 0233 end % Timewindow loop END 0234