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

vbmeg_current_reconstruct

PURPOSE ^

Current reconstruction using Bayesian inverse filter

SYNOPSIS ^

function [Jact ,ix_act, Jbck, ix_bck, Var]= vbmeg_current_reconstruct(proj_root, bayesfile,ix_area, trial_average, tsubsamp, mode)

DESCRIPTION ^

 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)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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