Current reconstruction for each trial using Bayesian inverse filter
--- Syntax
[Jact ,ix_act, Jbck, ix_bck, Varact, JactInfo] ...
= vbmeg_current_reconstruct_tr(proj_root, bayesfile, ...
ix_area, tsubsamp, jactdir)
This program output current files for each trial to specified directory
If trial average current is needed, use 'vbmeg_current_reconstruct'
--- Input
bayesfile : result file obtained by VB estimation
--- Optional Input
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
jactdir : output directory of current files for each trial
relative path from proj_root
do not include '/' or '\' in the end of string
If 'jactdir' is not given, jactdir = 'temp'
---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,Ntrial)
Nact : # of active region,
Nsample : # of time sample,
Ntrial : # of trials
Jbck : background current
ix_bck : Vertex index corresponding to background current
same for 'Jact'
Varact : posterior variance of act_area
When active area is empty,
posterior variance for background area is calculated
----
MEG data : B{Nsession}(Nsensors, Tsample, Ntrials)
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
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, JactInfo] = ... 0002 vbmeg_current_reconstruct_tr(proj_root, bayesfile, ix_area, ... 0003 tsubsamp, jactdir, mode) 0004 % Current reconstruction for each trial using Bayesian inverse filter 0005 % --- Syntax 0006 % [Jact ,ix_act, Jbck, ix_bck, Varact, JactInfo] ... 0007 % = vbmeg_current_reconstruct_tr(proj_root, bayesfile, ... 0008 % ix_area, tsubsamp, jactdir) 0009 % This program output current files for each trial to specified directory 0010 % If trial average current is needed, use 'vbmeg_current_reconstruct' 0011 % 0012 % --- Input 0013 % bayesfile : result file obtained by VB estimation 0014 % 0015 % --- Optional Input 0016 % ix_area : Position index to calculate estimated current 0017 % If 'ix_area' is empty or not given, 0018 % currents in the active region are calculated 0019 % tsubsamp : Specify subsampled time index 0020 % If 'tsubsamp' is empty or not given, 0021 % time subsampling is not done 0022 % jactdir : output directory of current files for each trial 0023 % relative path from proj_root 0024 % do not include '/' or '\' in the end of string 0025 % If 'jactdir' is not given, jactdir = 'temp' 0026 % 0027 % ---Output 0028 % Jact : active current 0029 % ix_act : Vertex index corresponding to active current 0030 % 0031 % Jact(n,t,:) is the current at the vertex 'ix_act(n)' & the time 't' 0032 % Jact(Nact,Nsample,Ntrial) 0033 % Nact : # of active region, 0034 % Nsample : # of time sample, 0035 % Ntrial : # of trials 0036 % 0037 % Jbck : background current 0038 % ix_bck : Vertex index corresponding to background current 0039 % same for 'Jact' 0040 % Varact : posterior variance of act_area 0041 % When active area is empty, 0042 % posterior variance for background area is calculated 0043 % ---- 0044 % MEG data : B{Nsession}(Nsensors, Tsample, Ntrials) 0045 % 0046 % 2004-9-25 M. Sato 0047 % 2004-12-5 M. Sato : Modified 0048 % 2004-12-22 M. Sato : Modified 0049 % 2005-1-5 M. Sato Modified 0050 % 2005-1-15 M. Sato Modified 0051 % 2005-08-19 O. Yamashita 0052 % * Posterior variance is calculated simultaneously 0053 % * Variables are renamed 0054 % * Corresponding to new bayesfile format 0055 % 2005-08-19 O. Yamashita ver.30b 0056 % 2005-09-30 O. Yamashita 0057 % * Minor bug fix when Lact >= 2 & Lact_var = 1 0058 % 2006-08-31 M. Sato 0059 % * Remove focal area from global area 0060 % * VBMEG result & liedfield is loaded by vb_load_result 0061 % * Inverse filter calculation is done in vb_invfilter_calc 0062 % 2009-04-02 Taku Yoshioka 0063 % Parameter name changed within this code for readability 0064 % (just replacing 'resultfile' to bayesfile) 0065 % 0066 % Copyright (C) 2011, ATR All Rights Reserved. 0067 % License : New BSD License(see VBMEG_LICENSE.txt) 0068 0069 0070 if ~exist('ix_area','var'), ix_area = [];end; 0071 if ~exist('tsubsamp','var'), tsubsamp = [];end; 0072 if ~exist('jactdir','var'), jactdir = 'temp'; end; 0073 if ~exist('mode','var'), mode = 0; end; 0074 0075 0076 %%%%%% load result 0077 if ~isempty(proj_root) 0078 jactdir_ab = fullfile(proj_root, jactdir); 0079 else 0080 jactdir_ab = [jactdir]; 0081 end 0082 0083 %check directory 0084 if exist(jactdir_ab) ~= 0 0085 msg = [ 'jact_dir : ' jactdir_ab ' exists. Do you over write?']; 0086 btn = questdlg(msg, 'confirm', 'Yes', 'No', 'Yes'); 0087 if strcmp(btn, 'No') 0088 error('processing aborted.'); 0089 return; 0090 end 0091 end 0092 0093 if mkdir(proj_root, jactdir) ~= 1, error('mkdir failed.'); end 0094 0095 %%%%%% load result and MEG data 0096 [B, Model, nGact, nGall, nGGall, Cov, Wact, Wbck, ... 0097 bsnorm, ix_act, ix_bck, vb_parm] ... 0098 = vb_load_result(proj_root, bayesfile, ix_area); 0099 0100 % B : MEG data 0101 % Model : Estimated model 0102 % nGact : Normalized leadfield of focal window 0103 % nGall : Normalized leadfield of global window 0104 % Cov : Sensor noise covariance matrix 0105 % bsnorm : Normalizing constant 0106 % Wact : Spatial smoothing matrix of focal window 0107 % Wbck : Spatial smoothing matrix of global window 0108 % ix_act : Vertex index corresponding to active current 0109 % ix_bck : Vertex index corresponding to background current 0110 0111 % # of vertex 0112 Nvact = vb_parm.Nvact; % # of active vertex 0113 Njact_area = vb_parm.Njact_area; 0114 Njall_area = vb_parm.Njall_area; 0115 0116 %% Constant 0117 Nsession = vb_parm.Nsession; % Number of sessions 0118 Ntrials = zeros(Nsession,1); 0119 Tsample = vb_parm.Tsample; % Number of total time sample 0120 Twindow = vb_parm.Twindow; 0121 0122 % Temporal subsampling index 0123 if isempty(tsubsamp), 0124 tsubsamp = 1:Tsample; 0125 end 0126 0127 % Temporal smoothing window weight 0128 [Tweight ,Tindex, Nindex, Nsample] = ... 0129 vb_overlapped_timewindow_weight(Twindow, Tsample, tsubsamp, mode); 0130 0131 Nwindow = length(Nindex); % # of time window 0132 0133 % Nsample : # of subsampled data 0134 0135 if mode == 1,fprintf('Non-overlapped concatenation mode\n'); end; 0136 0137 %% Posterior variance 0138 if Njact_area > 0 0139 Var = zeros(Njact_area,Nsample); 0140 Varact = zeros(Njact_area,Nsample); 0141 elseif Njall_area > 0 0142 Var = zeros(Njall_area,Nsample); 0143 Varact = zeros(Njall_area,Nsample); 0144 end 0145 0146 %% Current 0147 Jact_ave = zeros(Njact_area,Nsample); 0148 Jbck_ave = zeros(Njall_area,Nsample); 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 %% Inverse filter 0163 KW = cell(Nwindow,1); 0164 KW0 = cell(Nwindow,1); 0165 0166 Ntrial_all = 0; 0167 0168 % Session loop 0169 for n=1:Nsession 0170 Ntry = size(B{n},3); 0171 Ntrials(n) = Ntry; 0172 0173 % Following notation is used in comments 0174 % G0 : original leadfield matrix (non-smoothed) : [] 0175 % A : variance of active current : [A_inv] 0176 % tau : variance of active current : [v_inv] 0177 % Wa : spatial smoothing matrix of active current : [Wact] 0178 % Wb : spatial smoothing matrix of background current : [Wback] 0179 % Ga : leadfield of active current (smoothed) G0*Wa : [G] 0180 % Gb : leadfield of background current (smoothed) G0*Wb : [] 0181 0182 % Lead field for each session 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 for j=1:Nwindow 0191 % Time window loop 0192 % Subsampling time index 0193 Nid = Nindex{j}; % index in the total subsampled data 0194 Nt = length(Nid); 0195 0196 if isempty(Nid), continue; end; 0197 0198 [KWt, KW0t, var_tmp] = vb_invfilter_calc(a_inv(:,j), v_inv(j), ... 0199 G, GbW, GG, Covs, Wact, Wbck); 0200 KW{j} = KWt; 0201 KW0{j} = KW0t; 0202 var_tmp = Model.sx(j) * var_tmp./(bsnorm^2); 0203 Varact(:,Nid) = Varact(:,Nid) + var_tmp * Tweight{j}; 0204 end % Timewindow loop END 0205 0206 %%%% Current reconstruction 0207 for m=1:Ntry 0208 Jact = zeros(Njact_area,Nsample); 0209 Jbck = zeros(Njall_area,Nsample); 0210 0211 for j=1:Nwindow 0212 % Subsampling time index 0213 Tid = Tindex{j}; % subsampled time index 0214 Nid = Nindex{j}; % index in the total subsampled data 0215 0216 if isempty(Nid), continue; end; 0217 0218 weight = repmat(Tweight{j}, [Nch 1]); 0219 0220 Bt = B{n}(:,Tid,m) .* weight; % Time window smoothing 0221 0222 Jact(:,Nid) = Jact(:,Nid) + (KW{j} * Bt)/bsnorm; 0223 Jbck(:,Nid) = Jbck(:,Nid) + (KW0{j}* Bt)/bsnorm; 0224 end 0225 0226 % Trial average current 0227 Jact_ave = Jact_ave + Jact; 0228 Jbck_ave = Jbck_ave + Jbck; 0229 Ntrial_all = Ntrial_all + 1; 0230 0231 fname = sprintf('data_s%04dt%04d',n,m); 0232 vb_save([jactdir_ab '/' fname], 'Jact','Jbck', ... 0233 'Jbck','ix_act','ix_bck','Varact'); 0234 fprintf('progress... session:%04d/%04d, trial:%04d/%04d',... 0235 n,Nsession,m,Ntry); 0236 %fprintf('%s',repmat(char(sprintf('\b')),1,46)); 0237 fprintf('\n'); 0238 end; % Trial loop END 0239 0240 Var = Var + Ntry./max(Varact,eps); 0241 end; % Session loop END 0242 0243 JactInfo.Ntrial = Ntrials; 0244 JactInfo.Nsession = Nsession; 0245 JactInfo.jactdir = jactdir; 0246 JactInfo.NJact = Njact_area; 0247 JactInfo.Tsample = Tsample; 0248 0249 % Trial average current 0250 Jact = Jact_ave/Ntrial_all; 0251 Jbck = Jbck_ave/Ntrial_all; 0252 0253 %Var = Var/Ntrial_all; 0254 Var = 1./Var;