Home > functions > job > vb_job_ersp.m

vb_job_ersp

PURPOSE ^

Calculate ERSP for cortical current.

SYNOPSIS ^

function vb_job_ersp(proj_root,tf_parm)

DESCRIPTION ^

 Calculate ERSP for cortical current.
 (VBMEG public function)

 This function calculates event-related spectral perturbation (ERSP;
 Delorme and Makeig, 2004) for cortical current. The value is
 baseline-normalized and its unit is dB (10log10). Analyses can be done
 on the original subject's cortex or on the other subject's cortex
 (commonly template cortical model like fsaverage). 

 [syntax]
 vb_job_ersp(proj_root,tf_parm)

 [input]
 proj_root: <<string>> VBMEG project root directory. 
 tf_parm  : <<struct>> Parameters for time-frequency analysis
 --- fields of fs_parm
  curr_file     : <<string>> Cortical current file. 
  curr_file_base: <optional> <<string>> Cortical current files from
                  which baseline current amplitude is calculated. If not
                  given, curr_file is used to calculate baseline. 
  brain_file    : <<string>> Cortical surface model file on which
                  cortical current was calculated (subject's brain). 
  key           : <<string>> ID of coregistration matrix for each subjects
                  corresponding to brain files. In common, coregistration
                  matrix from the subject to a template cortical model is
                  specified. 
  tbrain_file   : <<string>> Cortical surface model filename of template
                  brain. This file is used to reduce the number of
                  vertices. 
  reduce        : <<int or double>> Vertex reduction ratio 
                  (@see reducepatch). 
  fmin          : <<double>> Minimum frequency.
  fmax          : <<double>> Maximum frequency.
  baseline      : <<double vector>> Start and end of baseline time window
                  in millisecond. 
  elim_time     : <optional> <<double>> In order to eliminate edge
                  effect, time points with length (elim_time)*(original
                  time length) are eliminated from start and end of data
                  after time-frequency decomposition (wavelet). Default
                  value is 0.05. 
  tf_file_out   : <<string>> Output filename. 
  timefreq_parm : <optional> <<struct>> EEGLAB timefreq parameters. If
                  this parameter is given, elim_time is ignored. 
  --- fields of timefreq_parm
   cycles        : <<scalar or vector>> Parameter of timefreq function. 
   winsize       : <<scalar>> Parameter of timefreq function. 
   padratio      : <<scalar>> Parameter of timefreq function. 
   ntimesout     : <<int>> Parameter of timefreq function.
   nfreqs        : <<int>> Parameter of timefreq function. 
   ntimesout_base: <<int>> Parameter of timefreq function for baseline. 
  --- 
 ---

 [output]
 TimeFrequency file (.tf.mat) is created. This file has the following
 variables: 
 data  : <<matrix>> Time-frequency data. F-by-T-by-I matrix, where F, T, 
         I are number of frequency bins (=length(TFinfo.freq), time
         window length (=length(TFinfo.Tmsec) and the number of vertices
         (=length(TFinfo.ix_act)), respectively.
 TFinfo: <<struct>> Information of time-frequency data. 
 --- fields of TFinfo
  tf_type   : <<string>> Set to 'ERSP'.
  Tsample   : <<vector>> Timepoint indices of the original data. 
  Pretrigger: <<int>> Pretrigger period (timepoints). 
  SampleFreq: <<float>> Sampling frequency of the original data. 
  Tmsec     : <<vector>> Actual time of each timepoints in msec
              (length(Tsample) = length(Tmsec)). 
  Ntrial    : <<vector>> Number of trials of each current data. 
  freq      : <<vector>> Frequency vector [Hz]. Each row of
              'data(:,f,:)' has frequency TFinfo.freq(f). 
  ix_act    : <<int vector>> Vertex indices on which time-frequency data
              is calculated.
  Wact      : <<double matrix>> Smoothing filter matrix. The size of
              smoothing filter is automatically determined from mean
              distance between vertices after reduction
              (tf_parm.reduce). This matrix will be used for displaying
              purpose. 
  ix_act_ex : <<int vector>> Vertex indices corresponding to the row of
              the smoothing matrix TFinfo.Wact. 
  base      : <<double matrix>> Baseline at each frequency and
              vertex. It is temporal average of the baseline ERSP. 
 ---

 [example]
 The following code is an example to calculate ERSP on template cortical
 surface model. 
 >> tf_parm.curr_file = './subjects/0001/result/0001_LR.curr.mat';
 >> tf_parm.curr_file_base = tf_parm.curr_file;
 >> tf_parm.brain_file = './subjects/0001/brain/0001.brain.mat';
 >> tf_parm.key = 'fsaverage.sphere.reg';
 >> tf_parm.tbrain_file ...
     = './subjects/fsaverage/brain/fsaverage.brain.mat';
 >> tf_parm.reduce = 0.2;
 >> tf_parm.fmin = 2.0;
 >> tf_parm.baseline = [-300 0];
 >> tf_parm.tf_file_out ...
     = './subjects/fsaverage/result/fsaverage_LR_ERSP.tf.mat';
 >> vb_job_ersp(proj_root,tf_parm);

 [history]
 2011-01-02 Taku Yoshioka
 2011-01-25 taku-y
  [enhancement] EEGLAB timefreq function supported. 
 2011-02-03 taku-y
  [trivial] Baseline ERSP 'base' added to TFinfo. 
 2011-02-28 taku-y
  [debug] [0 fmax] -> [fmin fmax] in timefreq().
 2011-07-01 taku-y
  [minor] vb_disp() was replaced by vb_disp_nonl() for displaying
  progress messages. 

 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 vb_job_ersp(proj_root,tf_parm)
0002 % Calculate ERSP for cortical current.
0003 % (VBMEG public function)
0004 %
0005 % This function calculates event-related spectral perturbation (ERSP;
0006 % Delorme and Makeig, 2004) for cortical current. The value is
0007 % baseline-normalized and its unit is dB (10log10). Analyses can be done
0008 % on the original subject's cortex or on the other subject's cortex
0009 % (commonly template cortical model like fsaverage).
0010 %
0011 % [syntax]
0012 % vb_job_ersp(proj_root,tf_parm)
0013 %
0014 % [input]
0015 % proj_root: <<string>> VBMEG project root directory.
0016 % tf_parm  : <<struct>> Parameters for time-frequency analysis
0017 % --- fields of fs_parm
0018 %  curr_file     : <<string>> Cortical current file.
0019 %  curr_file_base: <optional> <<string>> Cortical current files from
0020 %                  which baseline current amplitude is calculated. If not
0021 %                  given, curr_file is used to calculate baseline.
0022 %  brain_file    : <<string>> Cortical surface model file on which
0023 %                  cortical current was calculated (subject's brain).
0024 %  key           : <<string>> ID of coregistration matrix for each subjects
0025 %                  corresponding to brain files. In common, coregistration
0026 %                  matrix from the subject to a template cortical model is
0027 %                  specified.
0028 %  tbrain_file   : <<string>> Cortical surface model filename of template
0029 %                  brain. This file is used to reduce the number of
0030 %                  vertices.
0031 %  reduce        : <<int or double>> Vertex reduction ratio
0032 %                  (@see reducepatch).
0033 %  fmin          : <<double>> Minimum frequency.
0034 %  fmax          : <<double>> Maximum frequency.
0035 %  baseline      : <<double vector>> Start and end of baseline time window
0036 %                  in millisecond.
0037 %  elim_time     : <optional> <<double>> In order to eliminate edge
0038 %                  effect, time points with length (elim_time)*(original
0039 %                  time length) are eliminated from start and end of data
0040 %                  after time-frequency decomposition (wavelet). Default
0041 %                  value is 0.05.
0042 %  tf_file_out   : <<string>> Output filename.
0043 %  timefreq_parm : <optional> <<struct>> EEGLAB timefreq parameters. If
0044 %                  this parameter is given, elim_time is ignored.
0045 %  --- fields of timefreq_parm
0046 %   cycles        : <<scalar or vector>> Parameter of timefreq function.
0047 %   winsize       : <<scalar>> Parameter of timefreq function.
0048 %   padratio      : <<scalar>> Parameter of timefreq function.
0049 %   ntimesout     : <<int>> Parameter of timefreq function.
0050 %   nfreqs        : <<int>> Parameter of timefreq function.
0051 %   ntimesout_base: <<int>> Parameter of timefreq function for baseline.
0052 %  ---
0053 % ---
0054 %
0055 % [output]
0056 % TimeFrequency file (.tf.mat) is created. This file has the following
0057 % variables:
0058 % data  : <<matrix>> Time-frequency data. F-by-T-by-I matrix, where F, T,
0059 %         I are number of frequency bins (=length(TFinfo.freq), time
0060 %         window length (=length(TFinfo.Tmsec) and the number of vertices
0061 %         (=length(TFinfo.ix_act)), respectively.
0062 % TFinfo: <<struct>> Information of time-frequency data.
0063 % --- fields of TFinfo
0064 %  tf_type   : <<string>> Set to 'ERSP'.
0065 %  Tsample   : <<vector>> Timepoint indices of the original data.
0066 %  Pretrigger: <<int>> Pretrigger period (timepoints).
0067 %  SampleFreq: <<float>> Sampling frequency of the original data.
0068 %  Tmsec     : <<vector>> Actual time of each timepoints in msec
0069 %              (length(Tsample) = length(Tmsec)).
0070 %  Ntrial    : <<vector>> Number of trials of each current data.
0071 %  freq      : <<vector>> Frequency vector [Hz]. Each row of
0072 %              'data(:,f,:)' has frequency TFinfo.freq(f).
0073 %  ix_act    : <<int vector>> Vertex indices on which time-frequency data
0074 %              is calculated.
0075 %  Wact      : <<double matrix>> Smoothing filter matrix. The size of
0076 %              smoothing filter is automatically determined from mean
0077 %              distance between vertices after reduction
0078 %              (tf_parm.reduce). This matrix will be used for displaying
0079 %              purpose.
0080 %  ix_act_ex : <<int vector>> Vertex indices corresponding to the row of
0081 %              the smoothing matrix TFinfo.Wact.
0082 %  base      : <<double matrix>> Baseline at each frequency and
0083 %              vertex. It is temporal average of the baseline ERSP.
0084 % ---
0085 %
0086 % [example]
0087 % The following code is an example to calculate ERSP on template cortical
0088 % surface model.
0089 % >> tf_parm.curr_file = './subjects/0001/result/0001_LR.curr.mat';
0090 % >> tf_parm.curr_file_base = tf_parm.curr_file;
0091 % >> tf_parm.brain_file = './subjects/0001/brain/0001.brain.mat';
0092 % >> tf_parm.key = 'fsaverage.sphere.reg';
0093 % >> tf_parm.tbrain_file ...
0094 %     = './subjects/fsaverage/brain/fsaverage.brain.mat';
0095 % >> tf_parm.reduce = 0.2;
0096 % >> tf_parm.fmin = 2.0;
0097 % >> tf_parm.baseline = [-300 0];
0098 % >> tf_parm.tf_file_out ...
0099 %     = './subjects/fsaverage/result/fsaverage_LR_ERSP.tf.mat';
0100 % >> vb_job_ersp(proj_root,tf_parm);
0101 %
0102 % [history]
0103 % 2011-01-02 Taku Yoshioka
0104 % 2011-01-25 taku-y
0105 %  [enhancement] EEGLAB timefreq function supported.
0106 % 2011-02-03 taku-y
0107 %  [trivial] Baseline ERSP 'base' added to TFinfo.
0108 % 2011-02-28 taku-y
0109 %  [debug] [0 fmax] -> [fmin fmax] in timefreq().
0110 % 2011-07-01 taku-y
0111 %  [minor] vb_disp() was replaced by vb_disp_nonl() for displaying
0112 %  progress messages.
0113 %
0114 % Copyright (C) 2011, ATR All Rights Reserved.
0115 % License : New BSD License(see VBMEG_LICENSE.txt)
0116 
0117 vb_disp('--- ERSP calculation');
0118 
0119 %
0120 % Verbose level
0121 %
0122 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0123 const = vb_define_verbose;
0124 vb_struct2vars(const,{'VERBOSE_LEVEL_NONE','VERBOSE_LEVEL_EMERGENCY', ...
0125                     'VERBOSE_LEVEL_WARNING','VERBOSE_LEVEL_NOTICE', ...
0126                     'VERBOSE_LEVEL_INFO','VERBOSE_LEVEL_DEBUG'});
0127 VERBOSE_LEVEL_DEBUG = const.VERBOSE_LEVEL_DEBUG;
0128 
0129 %
0130 % Input arguments
0131 %
0132 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0133 vb_struct2vars(tf_parm,{'curr_file','brain_file','key', ...
0134                     'tf_file_out','tbrain_file','reduce','fmin','fmax', ...
0135                     'baseline','elim_time','curr_file_base', ...
0136                     'timefreq_parm'});
0137 
0138 if isempty(curr_file_base), curr_file_base = curr_file; end
0139 if isempty(elim_time), elim_time = 0.05; end
0140 
0141 %
0142 % Check arguments
0143 %
0144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0145 % Sampling frequency
0146 Jinfo1 = vb_load_current_info(curr_file);
0147 Jinfo2 = vb_load_current_info(curr_file_base);
0148 if Jinfo1.SampleFreq~=Jinfo2.SampleFreq, 
0149   error(['Sampling frequencies of current file and baseline ' ...
0150          'current file must be the same.']);
0151 end
0152 
0153 %
0154 % Preparation of cortical dipole indices
0155 %
0156 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0157 % Number of dipoles on which TF analysis conducted.
0158 if ~isempty(key), % use template brain
0159   % note: consistency check between tbrain_file and key{1} is required,
0160   % but postponed. (7 Dec. 2010, taku-y)
0161   brain_file_tmp = tbrain_file;
0162 else % single subject analysis
0163   brain_file_tmp = brain_file;
0164 end
0165 
0166 % Reduce vertex of the cortical model
0167 V = vb_load_cortex(brain_file_tmp);
0168 ix_act = vb_reduce_cortex(brain_file_tmp,1:size(V,1),reduce);
0169 I = length(ix_act);
0170 clear V;
0171 
0172 % For single subject analysis, find overlap between ix_act and
0173 % Jinfo.ix_act_ex. ix_act is replaced by it and I is overwritten.
0174 if isempty(key), % single subject analysis
0175   Jinfo = vb_load_current_info(curr_file);
0176   [tmp,ixx] = intersect(Jinfo.ix_act_ex,ix_act);
0177   ix_act = Jinfo.ix_act_ex(ixx);
0178   I = length(ix_act);
0179   Wact = Jinfo.Wact(ixx,:);
0180 end
0181 
0182 % Calculate mean distance between reduced vertices to determine smoothing
0183 % filter size
0184 R = 0.0;
0185 N = 0;
0186 [nextDD,nextIX] = vb_load_cortex_neighbour(brain_file_tmp);
0187 for i=1:length(ix_act)
0188   [tmp,ixx] = union(nextIX{ix_act(i)},ix_act);
0189   R = R+sum(nextDD{ix_act(i)}(ixx));
0190   N = N+length(ixx);
0191 end
0192 R = R/N;
0193 
0194 % Calculate smoothing filter. It will be used for displaying purpose, not
0195 % affect further analysis.
0196 V = vb_load_cortex(brain_file_tmp);
0197 [W,ix_act_ex] ...
0198     = vb_spatial_gauss_filter(brain_file_tmp,R,2*R,ix_act,1,-1,1);
0199 W = W./repmat(sum(W,2),[1 size(W,2)]);
0200 TFinfo.Wact = W;
0201 TFinfo.ix_act_ex = ix_act_ex;
0202 clear V;
0203 
0204 %
0205 % Preparation of variables time-frequency data stored
0206 %
0207 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0208 if isempty(timefreq_parm), % wavelet() parameters
0209   % Number of timepoints
0210   Jinfo = vb_load_current_info(curr_file);
0211   T  = length(Jinfo.Tmsec);    % time points
0212   DJ = 0.25;                   % scale of time points spacing
0213   DT = 1/(Jinfo.SampleFreq);   % original time points spacing
0214   S0 = 2*DT;                   % smallest freq. scale
0215 
0216   % Number of frequency bins
0217   [tmp,period,coi] = wavelet(randn(1,T),DT,1,DJ);
0218   TFinfo.freq      = 1./period;
0219   ix_freq          = find(TFinfo.freq>=fmin);
0220   TFinfo.freq      = TFinfo.freq(ix_freq);
0221   F                = length(TFinfo.freq);
0222 
0223   % Eliminate start and end time points
0224   t_elim  = ceil(T*elim_time);
0225   ix_time = (t_elim+1):(T-t_elim);
0226   Tmsec   = Jinfo.Tmsec(ix_time);
0227   T       = length(Tmsec);
0228 
0229   % Whole data
0230   data       = zeros(F,T,I);
0231 
0232   % Baseline
0233   if ~isempty(baseline), 
0234     Jinfo_base      = vb_load_current_info(curr_file_base);
0235     Tbase           = length(Jinfo_base.Tmsec);
0236     t_elim_base     = ceil(Tbase*elim_time);
0237     ix_time_base    = (t_elim_base+1):(Tbase-t_elim_base);
0238     Tmsec_base      = Jinfo_base.Tmsec(ix_time_base);
0239     Tbase           = length(Tmsec_base);
0240     base            = zeros(F,Tbase,I);
0241     Ntrial_base     = Jinfo_base.Ntrial;
0242   else
0243     Tbase           = [];
0244   end
0245 else % timefreq() parameters
0246   vb_struct2vars(timefreq_parm,{'cycles','winsize','padratio', ...
0247                       'ntimesout','nfreqs','ntimesout_base'});
0248   
0249   % Number of timepoints and frequency bins
0250   Jinfo  = vb_load_current_info(curr_file);
0251   frames = length(Jinfo.Tmsec);
0252   tlimits = [min(Jinfo.Tmsec) max(Jinfo.Tmsec)];
0253   srate   = Jinfo.SampleFreq;
0254   if isempty(fmax), 
0255     fmax = srate/2;
0256   end
0257   
0258   [tmp,freqs,Tmsec] ...
0259       = timefreq(randn(frames,1),srate,'cycles',cycles, ...
0260                  'freqs',[fmin fmax],'padratio',padratio, ...
0261                  'winsize',winsize,'tlimits',tlimits, ...
0262                  'ntimesout',ntimesout,'nfreqs',nfreqs, ...
0263                  'verbose','off');
0264   
0265   T                = length(Tmsec);
0266   ix_time          = 1:length(Tmsec);
0267   TFinfo.freq      = freqs;
0268   ix_freq          = find(TFinfo.freq>=fmin);
0269   TFinfo.freq      = TFinfo.freq(ix_freq);
0270   F                = length(TFinfo.freq);
0271   
0272   % Whole data
0273   data = zeros(F,T,I);
0274   
0275   % Baseline
0276   if ~isempty(baseline),
0277     Jinfo_base   = vb_load_current_info(curr_file_base);
0278     frames_base  = length(Jinfo_base.Tmsec);
0279     tlimits_base = [min(Jinfo.Tmsec) max(Jinfo.Tmsec)];
0280     
0281     [tmp,freqs_base,Tmsec_base] ...
0282       = timefreq(randn(frames_base,1),srate,'cycles',cycles, ...
0283                  'freqs',[fmin fmax],'padratio',padratio, ...
0284                  'winsize',winsize,'tlimits',tlimits_base, ...
0285                  'ntimesout',ntimesout_base,'nfreqs',nfreqs, ...
0286                  'verbose','off');
0287     
0288     Tbase           = length(Tmsec_base);
0289     ix_time_base    = 1:Tbase;
0290     base            = zeros(F,Tbase,I);
0291     Ntrial_base     = Jinfo_base.Ntrial;
0292   else
0293     Tbase           = [];
0294   end
0295 end
0296 
0297 vb_disp(['Number of frequency bins       : ' num2str(F)]);
0298 vb_disp(['Number of timepoints           : ' num2str(T)]);
0299 vb_disp(['Number of vertices             : ' num2str(I)]);
0300 
0301 if ~isempty(Tbase),
0302   vb_disp(['Number of timepoints (baseline): ' num2str(Tbase)]);
0303 end
0304 
0305 %
0306 % Main routine
0307 %
0308 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0309 Jinfo      = vb_load_current_info(curr_file);
0310 Ntrial     = Jinfo.Ntrial;
0311 
0312 if ~isempty(key), % use template brain
0313   C = vb_load_cortex_pmat(brain_file,key);
0314 
0315   if strcmp(Jinfo.curr_type,'Z-current'), 
0316     CW = sparse(C(ix_act,Jinfo.ix_act_ex)*Jinfo.Wact);
0317   else
0318     CW = sparse(C(ix_act,Jinfo.ix_act_ex));
0319   end
0320 else % single subject analysis
0321   CW = Wact;
0322 end
0323 
0324 prg = 0;
0325 prg_all = Ntrial;
0326 prg_str = sprintf('Current : %d%% processed',0);
0327 vb_disp_nonl(prg_str);
0328 
0329 for j=1:Ntrial % trial loop
0330   [tmp,Z] = vb_load_current(curr_file,Jinfo.curr_type,false,j);
0331     
0332   for k=1:I % vertex loop
0333     ix = find(abs(CW(k,:))>0);
0334     J = CW(k,ix)*Z(ix,:); % J-current at ix_act-th dipole
0335     
0336     if isempty(timefreq_parm), 
0337       wave = wavelet(J,DT,1,DJ);
0338     else
0339       wave = timefreq(J',srate,'cycles',cycles, ...
0340                       'freqs',[fmin fmax],'padratio',padratio, ...
0341                       'winsize',winsize,'tlimits',tlimits, ...
0342                       'ntimesout',ntimesout,'nfreqs',nfreqs, ...
0343                       'verbose','off');
0344     end
0345 
0346     wave = wave(ix_freq,ix_time);
0347     data(:,:,k) = data(:,:,k)+abs(wave).^2;
0348   end
0349   
0350   % progress message
0351   prg = prg+1;
0352   for ii=1:length(prg_str); vb_disp_nonl(sprintf('\b')); end
0353   prg_str = sprintf('Current : %d%% processed',ceil(100*(prg/prg_all)));
0354   vb_disp_nonl(prg_str);
0355 end
0356 
0357 for ii=1:length(prg_str); vb_disp_nonl(sprintf('\b')); end
0358 vb_disp('Current : finished');
0359 
0360 % Baseline
0361 Jinfo_base = vb_load_current_info(curr_file_base);
0362 Ntrial     = Jinfo_base.Ntrial;
0363 
0364 if ~isempty(key), % use template brain
0365   C = vb_load_cortex_pmat(brain_file,key);
0366 
0367   if strcmp(Jinfo_base.curr_type,'Z-current'), 
0368     CW = sparse(C(ix_act,Jinfo_base.ix_act_ex)*Jinfo_base.Wact);
0369   else
0370     CW = sparse(C(ix_act,Jinfo_base.ix_act_ex));
0371   end
0372 else % single subject analysis
0373   CW = Wact;
0374 end
0375 
0376 prg = 0;
0377 prg_all = Ntrial_base;
0378 prg_str = sprintf('Baseline: %d%% processed',0);
0379 vb_disp_nonl(prg_str);
0380 
0381 for j=1:Ntrial_base % trial loop
0382   [tmp,Z] = vb_load_current(curr_file_base, ...
0383                             Jinfo_base.curr_type,false,j);
0384 
0385   for k=1:I % vertex loop
0386     ix = find(abs(CW(k,:))>0);
0387     J = CW(k,ix)*Z(ix,:); % J-current at ix_act-th dipole
0388     if isempty(timefreq_parm), 
0389       wave = wavelet(J,DT,1,DJ);
0390     else
0391       wave = timefreq(J',srate,'cycles',cycles, ...
0392                       'freqs',[fmin fmax],'padratio',padratio, ...
0393                       'winsize',winsize,'tlimits',tlimits_base, ...
0394                       'ntimesout',ntimesout_base,'nfreqs',nfreqs, ...
0395                       'verbose','off');
0396     end
0397 
0398     wave = wave(ix_freq,ix_time_base);        
0399     base(:,:,k) = base(:,:,k)+abs(wave).^2;
0400   end
0401   
0402   % progress message
0403   prg = prg+1;
0404   for ii=1:length(prg_str); vb_disp_nonl(sprintf('\b')); end
0405   prg_str = sprintf('Baseline: %d%% processed',ceil(100*(prg/prg_all)));
0406   vb_disp_nonl(prg_str);
0407 end
0408 
0409 for ii=1:length(prg_str); vb_disp_nonl(sprintf('\b')); end
0410 vb_disp('Baseline: finished');
0411 
0412 % Normalization
0413 [tmp,ix_base_start] = min(abs(Tmsec_base-tf_parm.baseline(1)));
0414 [tmp,ix_base_end]   = min(abs(Tmsec_base-tf_parm.baseline(2)));
0415 ix_base = ix_base_start:ix_base_end;
0416 data = data./Ntrial;
0417 base = base./Ntrial_base;
0418 data = 10*log10(data./repmat(mean(base(:,ix_base,:),2),[1 T 1]));
0419 
0420 %
0421 % Make struct 'TFinfo'
0422 %
0423 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0424 Jinfo = vb_load_current_info(curr_file);
0425 
0426 TFinfo.tf_type    = 'ERSP';
0427 TFinfo.Tmsec      = Tmsec;
0428 TFinfo.Ntrial     = Ntrial;
0429 TFinfo.ix_act     = ix_act;
0430 TFinfo.base       = mean(base(:,ix_base,:),2);
0431 
0432 if isempty(timefreq_parm),
0433   TFinfo.Tsample    = Jinfo.Tsample(ix_time);
0434   TFinfo.Pretrigger = Jinfo.Pretrigger-t_elim;
0435   TFinfo.SampleFreq = Jinfo.SampleFreq;
0436   TFinfo.freq_scale = 'log'; % for wavelet()
0437 else
0438   Tmsec = TFinfo.Tmsec;
0439   [tmp,TFinfo.Pretrigger] = min(abs(Tmsec(:)));
0440   TFinfo.SampleFreq = ((Tmsec(end)-Tmsec(1))./1000)./length(Tmsec);
0441   TFinfo.freq_scale = 'linear';
0442 end
0443 
0444 %
0445 % Save result
0446 %
0447 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0448 vb_fsave(tf_file_out,'tf_parm','data','TFinfo');
0449 
0450 %
0451 % Save project file
0452 %
0453 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0454 proj_file = get_project_filename;
0455 if isempty(proj_file), 
0456   return;
0457 end
0458 
0459 project_file_mgr('load', proj_file);
0460 project_file_mgr('add', 'tf_parm', tf_parm);
0461 
0462 return;

Generated on Tue 27-Aug-2013 11:46:04 by m2html © 2005