vbmeg_feature_calc_z

 Current reconstruction using Bayesian inverse filter

 --- Syntax
  [Zact ,Jinfo, bayes_parm, vb_parm, MEGinfo] ...
          = vbmeg_feature_calc_z(proj_root, curr_parm);
 --- Input
 curr_parm.filterfile : VB filter file (.vbfilt.mat)
 curr_parm.currfile   : output file for estimated current feature (.curr.mat)
 curr_parm.megfile    : MEG/EEG data file
 
 curr_parm.freq_range : band frequency for feature extraction [Nband x 2]
          .freq_range(n ,:)  = [freq_low freq_high] Hz (for n-th band)
   For each time window specified by Tperiod & Tnext,
   sum of bandpass frequency component power (square of abs) are calculated.
 The following parameters are independent of bayes_parm used in VB filter
 If there are multiple time windows in VB filter,
    average filter is calculated before current calculation
 curr_parm.twin_meg:  Time range of meg/eeg data used for analysis
   e.g. twin_meg = [1 600]  
 curr_parm.Tperiod :  Length of the time windows. 
 curr_parm.Tnext   :  Moving steps of the time window. 
   e.g. Tperiod = 100, Tnext = 50  
               --> [1 100; 51 150; 101 200;...]
 --- Optional Input
 curr_parm.trial_average = ON : average current over all sessions 
                         = [OFF]  : current for each session

 curr_parm.ix_area : Position index to calculate estimated current
           If 'ix_area' is empty or not given, 
              currents in the active region are calculated
 ---Output
 Zact    : active current feature

 Zact(:, n + Nband*(j-1), :) : n-th bandpass component at j-th time window
 Zact(Nact,Nfeature)         for trial_average = ON 
 Zact(Nact,Nfeature,Ntrial)  for trial_average = OFF
   Nact     : # of active vertex, 
   Nfeature : # of feature = Nwindow * Nband
   Ntrial   : # of trials in all session]

 Jinfo.ix_act  : Vertex index corresponding to active current

 2008/12/2 M. Sato
 2010/04/22 (Sako) substituted BMI_load_eeginfo with vb_load_measurement_info
 ---

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

0001 function [Zact ,Jinfo, bayes_parm, vb_parm, MEGinfo] ...
0002           = vbmeg_feature_calc_z(proj_root, curr_parm);
0003 % Current reconstruction using Bayesian inverse filter
0004 %
0005 % --- Syntax
0006 %  [Zact ,Jinfo, bayes_parm, vb_parm, MEGinfo] ...
0007 %          = vbmeg_feature_calc_z(proj_root, curr_parm);
0008 % --- Input
0009 % curr_parm.filterfile : VB filter file (.vbfilt.mat)
0010 % curr_parm.currfile   : output file for estimated current feature (.curr.mat)
0011 % curr_parm.megfile    : MEG/EEG data file
0012 %
0013 % curr_parm.freq_range : band frequency for feature extraction [Nband x 2]
0014 %          .freq_range(n ,:)  = [freq_low freq_high] Hz (for n-th band)
0015 %   For each time window specified by Tperiod & Tnext,
0016 %   sum of bandpass frequency component power (square of abs) are calculated.
0017 % The following parameters are independent of bayes_parm used in VB filter
0018 % If there are multiple time windows in VB filter,
0019 %    average filter is calculated before current calculation
0020 % curr_parm.twin_meg:  Time range of meg/eeg data used for analysis
0021 %   e.g. twin_meg = [1 600]
0022 % curr_parm.Tperiod :  Length of the time windows.
0023 % curr_parm.Tnext   :  Moving steps of the time window.
0024 %   e.g. Tperiod = 100, Tnext = 50
0025 %               --> [1 100; 51 150; 101 200;...]
0026 % --- Optional Input
0027 % curr_parm.trial_average = ON : average current over all sessions
0028 %                         = [OFF]  : current for each session
0029 %
0030 % curr_parm.ix_area : Position index to calculate estimated current
0031 %           If 'ix_area' is empty or not given,
0032 %              currents in the active region are calculated
0033 % ---Output
0034 % Zact    : active current feature
0035 %
0036 % Zact(:, n + Nband*(j-1), :) : n-th bandpass component at j-th time window
0037 % Zact(Nact,Nfeature)         for trial_average = ON
0038 % Zact(Nact,Nfeature,Ntrial)  for trial_average = OFF
0039 %   Nact     : # of active vertex,
0040 %   Nfeature : # of feature = Nwindow * Nband
0041 %   Ntrial   : # of trials in all session]
0042 %
0043 % Jinfo.ix_act  : Vertex index corresponding to active current
0044 %
0045 % 2008/12/2 M. Sato
0046 % 2010/04/22 (Sako) substituted BMI_load_eeginfo with vb_load_measurement_info
0047 % ---
0048 %
0049 % Copyright (C) 2011, ATR All Rights Reserved.
0050 % License : New BSD License(see VBMEG_LICENSE.txt)
0051 
0052 if ~isempty(proj_root)
0053     filterfile = fullfile(proj_root, curr_parm.filterfile);
0054 else
0055     filterfile = curr_parm.filterfile;
0056 end  
0057 
0058 %%%%%% load VBMEG estimated result
0059 load(filterfile, 'VBfilt','Jinfo','bayes_parm','vb_parm');
0060 
0061 %%%%%% check parameter of 'curr_parm'
0062 % Values of 'curr_parm' fields dominates over
0063 %   those of 'bayes_parm' in bayesfile
0064 [bayes_parm, ix_area, trial_average, tsubsamp, overlap_mode, verbose] = ...
0065             check_arg(bayes_parm, curr_parm);
0066             
0067 if ~isempty(proj_root);
0068   bayes_parm_abs = vb_parm_absolute_path(proj_root, bayes_parm);
0069 else
0070   bayes_parm_abs = bayes_parm;
0071 end
0072 
0073 %%%%%% MEG data preparation
0074 % B      : MEG data
0075 [B, Ntrial, Nch, Tsample, Twindow] = ...
0076     vb_megdata_preparation(bayes_parm_abs);
0077 
0078 MEGinfo = vb_load_measurement_info(bayes_parm_abs.megfile{1});
0079 
0080 % sampling frequency
0081 Fs = MEGinfo.SampleFrequency;
0082 
0083 if ~isfield(curr_parm,'freq_range'), 
0084     freq_range = [5 9; 11 16; 21 26]
0085 else
0086     freq_range = curr_parm.freq_range; 
0087 end;
0088 if ~isfield(curr_parm,'fft_window'), 
0089     fft_window = 0;
0090 else
0091     fft_window = curr_parm.fft_window;
0092 end
0093 
0094 % Number of session
0095 Nsession = length(B);     
0096 fprintf('Number of session: %d\n', Nsession);
0097 
0098 % Temporal window
0099 Nwindow = size(Twindow,1);
0100 Tindex  = cell(Nwindow,1);
0101 
0102 for j=1:Nwindow
0103     Tindex{j} = Twindow(j,1):Twindow(j,2);
0104 end
0105 
0106 % FFT point number
0107 NFFT = 0;
0108 for j=1:Nwindow
0109     Nt   = length(Tindex{j});
0110     nfft = 2^nextpow2(Nt);
0111     if NFFT < nfft, NFFT=nfft; end;
0112 end
0113 
0114 % frequency for FFT
0115 df=Fs/NFFT;
0116 freq=0:df:df*(NFFT/2);
0117 
0118 if max(freq_range(:)) > Fs/2, 
0119     error('max freq is larger than Nyquist freq')
0120 end
0121 
0122 % band frequency
0123 [Nband ,D] = size(freq_range);
0124 if D ~=2, error('freq_range should be [Nband x 2]'); end;
0125 
0126 % index list for band frequency
0127 flist  = zeros(Nband ,2);
0128 for n=1:Nband*2
0129     [tmp, ix] = min(abs(freq - freq_range(n)));
0130     flist(n) = ix;
0131 end
0132 
0133 fxlist = [];
0134 ixlist = [];
0135 nxlist = [];
0136 Nlast  = 0;
0137 %Ilist  = zeros(Nband,2);
0138 
0139 for n=1:Nband
0140     Nlist  = flist(n,2) - flist(n,1) + 1;
0141     Nstart = Nlast + 1;
0142     Nlast  = Nlast + Nlist;
0143     
0144     fxlist = [fxlist, flist(n,1):flist(n,2)];
0145     ixlist = [ixlist, Nstart:Nlast];
0146     nxlist = [nxlist, repmat(n, [1, Nlist])];
0147     
0148 %    Ilist(n,:) = [Nstart Nlast];
0149 end
0150 
0151 % Sparse matrix to get sum of frequency component in each band
0152 Nall = length(ixlist);
0153 Wfreq = sparse( ixlist , nxlist , 1 , Nall, Nband) ;
0154 % Bf = fft( vb_repmultiply(B, hanning(Nt)') ,NFFT, 2);
0155 % Zf = Bf(:,fxlist) ).^2 * Wfreq
0156 % Zf(:,n) : n-th band power summed over [freq_range(n,1) freq_range(n,2)]
0157 
0158 % VBfilt.KW = VB inverse filter  (Njact,Nch,Nwindow,Ntask)
0159 [Njact,Nch,Nw,Ntask] = size(VBfilt.KW);
0160 
0161 if Nw==1 && Ntask==1, 
0162     KW = VBfilt.KW(:,:,1,1) ;
0163 else
0164     fprintf('(Nwindow=%d ,Ntask=%d)\n',Nw,Ntask)
0165     fprintf('Take average of VB filters\n')
0166     KW = squeeze(mean( mean( VBfilt.KW,3) ,4));
0167 end
0168 
0169 fprintf('Active vertex number = %d\n',Njact)
0170 fprintf('Active sensor number = %d\n',Nch)
0171 fprintf('Number of time window= %d\n',Nwindow)
0172 fprintf('Number of bandpass components= %d\n',Nband)
0173 fprintf('Number of FFT points= %d\n',NFFT)
0174 
0175 %
0176 % --- Set Current Info
0177 %
0178 
0179 %Jinfo.Pretrigger = MEGinfo.Pretrigger;
0180 %Jinfo.Wact = Wact;
0181 
0182 %Lact      = Jinfo.Lact;
0183 %ix_act    = Jinfo.ix_act;
0184 %ix_act_ex = Jinfo.ix_act_ex;
0185 %
0186 %Jinfo = [];
0187 %Jinfo.Lact       = Lact;
0188 %Jinfo.ix_act     = ix_act;
0189 %Jinfo.ix_act_ex  = ix_act_ex;
0190 
0191 VBfreq.KW         = KW;
0192 VBfreq.fxlist     = fxlist;
0193 VBfreq.Wfreq      = Wfreq;
0194 VBfreq.NFFT       = NFFT;
0195 VBfreq.Twindow    = Twindow;
0196 VBfreq.freq_range = freq_range;
0197 
0198 Jinfo.VBfreq     = VBfreq;
0199 
0200 Jinfo.Nwindow    = Nwindow;
0201 Jinfo.SampleFreq = Fs;
0202 
0203 % MEG time window index
0204 Tstart  = bayes_parm.twin_meg(1);
0205 Tend    = bayes_parm.twin_meg(2);
0206 if isempty(tsubsamp)
0207     Jinfo.Tsample = Tstart:Tend;
0208 else
0209     Jinfo.Tsample = tsubsamp + Tstart - 1;
0210 end
0211 
0212 %%%%%% Time window loop
0213 % MEG for each trial
0214 err  = 0;
0215 nall = 0;
0216 Zact = zeros( Njact, Nwindow*Nband, sum(Ntrial));
0217 
0218 for n=1:Nsession
0219     Ntry = Ntrial(n);
0220     
0221     if Ntry == 0, continue; end;
0222 
0223     Bs = B{n};
0224     
0225     for m=1:Ntry
0226         nall = nall + 1;
0227         for j=1:Nwindow
0228             % Subsampling time index
0229             Tid = Tindex{j};    % subsampled time index
0230             Nt  = length(Tid);
0231             
0232             if Nt == 0, continue; end;
0233             
0234             % FFT of MEG/EEG data with Hanning window
0235             Bt = Bs(:,Tid,m);
0236             
0237             if fft_window==1,
0238                 Bt = vb_repmultiply(Bt, hanning(Nt)');
0239             end
0240             
0241             Bf = fft( Bt ,NFFT, 2);
0242             
0243             % Current reconstruction in frequency space
0244             Zf = abs( KW * Bf(:,fxlist) ).^2 ;
0245             
0246             % Frequency band component
0247             Zact(:, (1:Nband) + (j-1)*Nband, nall) = Zf * Wfreq;
0248             
0249 %            % Debug
0250 %            Zb  = Zf * Wfreq;
0251 %            Zbb = zeros(Njact,Nband);
0252 %
0253 %            for k=1:Nband
0254 %                Zbb(:,k) = sum(Zf(:,Ilist(k,1):Ilist(k,2)) ,2);
0255 %            end
0256 %
0257 %            err = err + sum(sum(abs(Zb - Zbb)))/sum(sum(abs(Zbb)));
0258         end % Timewindow loop
0259     end; % Trial loop
0260 end
0261 
0262 if trial_average == ON
0263     Zact = mean(Zact,3);
0264 end
0265 
0266 Jinfo.Ntrial   = Ntrial;
0267 Jinfo.Nsession = Nsession;
0268 
0269 if verbose==1,
0270     fprintf('\n')
0271 end
0272 
0273 %err
0274 
0275 % ix_act : Vertex index corresponding to active current Zact
0276 % ix_act_ex : Vertex index corresponding to active current Jact
0277 % Wact   : Spatial smoothing matrix of focal window
0278 % Jact   = Wact * Zact
0279 
0280 return
0281 
0282 %%%% ---------------
0283 function [bayes_parm,ix_area,trial_average,tsubsamp,overlap_mode,verbose]= ...
0284             check_arg(bayes_parm,curr_parm)
0285 
0286 if isfield(curr_parm,'basisfile'), 
0287     bayes_parm.basisfile = curr_parm.basisfile;
0288 end;
0289 if isfield(curr_parm,'megfile'), 
0290     bayes_parm.megfile   = curr_parm.megfile  ;
0291 end;
0292 if isfield(curr_parm,'twin_meg'), 
0293     bayes_parm.twin_meg  = curr_parm.twin_meg ;
0294 end;
0295 if isfield(curr_parm,'Tperiod'), 
0296     bayes_parm.Tperiod   = curr_parm.Tperiod  ;
0297 end;
0298 if isfield(curr_parm,'Tnext'), 
0299     bayes_parm.Tnext     = curr_parm.Tnext    ;
0300 end;
0301 
0302 if ~isfield(curr_parm,'trial_average'), 
0303     trial_average = ON; 
0304 else
0305     trial_average = curr_parm.trial_average; 
0306 end;
0307 
0308 bayes_parm.trial_average = trial_average;
0309 
0310 if ~isfield(curr_parm,'ix_area'),  
0311     ix_area = []; 
0312 else
0313     ix_area = curr_parm.ix_area; 
0314 end;
0315 if ~isfield(curr_parm,'tsubsamp'), 
0316     tsubsamp = []; 
0317 else
0318     tsubsamp = curr_parm.tsubsamp; 
0319 end;
0320 if ~isfield(curr_parm,'overlap_mode'),     
0321     overlap_mode = 0; 
0322 else
0323     overlap_mode = curr_parm.overlap_mode; 
0324 end;
0325 if ~isfield(curr_parm,'verbose'),     
0326     verbose = 1; 
0327 else
0328     verbose = curr_parm.verbose; 
0329 end;