Home > vbmeg > functions > job > subdirectory > vb_job_itpc.m

vb_job_itpc

PURPOSE ^

Calculate ITPC for cortical current.

SYNOPSIS ^

function vb_job_itpc(proj_root,tf_parm)

DESCRIPTION ^

 Calculate ITPC for cortical current.
 (VBMEG public function)

 This function calculates inter-trial phase coherence (ITPC; Delorme and
 Makeig, 2004) for cortical current. The value is 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. 
  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 vb_reducepatch). 
  fmin          : <<double>> Minimum frequency.
  fmax          : <<double>> Maximum frequency.
  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. 
  --- 
 ---

 [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 'ITPC'.
  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. 
 ---

 [example]
 The following code is an example to calculate ITPC on template cortical
 surface model. 
 >> tf_parm.curr_file = './subjects/0001/result/0001_LR.curr.mat';
 >> 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.tf_file_out ...
     = './subjects/fsaverage/result/fsaverage_LR_ITPC.tf.mat';
 >> vb_job_itpc(proj_root,tf_parm);

 [history]
 2011-01-04 Taku Yoshioka
 2011-02-28 taku-y
  [enhancement] EEGLAB timefreq function supported. 
 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_itpc(proj_root,tf_parm)
0002 % Calculate ITPC for cortical current.
0003 % (VBMEG public function)
0004 %
0005 % This function calculates inter-trial phase coherence (ITPC; Delorme and
0006 % Makeig, 2004) for cortical current. The value is Analyses can be done
0007 % on the original subject's cortex or on the other subject's cortex
0008 % (commonly template cortical model like fsaverage).
0009 %
0010 % [syntax]
0011 % vb_job_ersp(proj_root,tf_parm)
0012 %
0013 % [input]
0014 % proj_root: <<string>> VBMEG project root directory.
0015 % tf_parm  : <<struct>> Parameters for time-frequency analysis
0016 % --- fields of fs_parm
0017 %  curr_file     : <<string>> Cortical current file.
0018 %  brain_file    : <<string>> Cortical surface model file on which
0019 %                  cortical current was calculated (subject's brain).
0020 %  key           : <<string>> ID of coregistration matrix for each subjects
0021 %                  corresponding to brain files. In common, coregistration
0022 %                  matrix from the subject to a template cortical model is
0023 %                  specified.
0024 %  tbrain_file   : <<string>> Cortical surface model filename of template
0025 %                  brain. This file is used to reduce the number of
0026 %                  vertices.
0027 %  reduce        : <<int or double>> Vertex reduction ratio
0028 %                  (@see vb_reducepatch).
0029 %  fmin          : <<double>> Minimum frequency.
0030 %  fmax          : <<double>> Maximum frequency.
0031 %  elim_time     : <optional> <<double>> In order to eliminate edge
0032 %                  effect, time points with length (elim_time)*(original
0033 %                  time length) are eliminated from start and end of data
0034 %                  after time-frequency decomposition (wavelet). Default
0035 %                  value is 0.05.
0036 %  tf_file_out   : <<string>> Output filename.
0037 %  timefreq_parm : <optional> <<struct>> EEGLAB timefreq parameters. If
0038 %                  this parameter is given, elim_time is ignored.
0039 %  --- fields of timefreq_parm
0040 %   cycles        : <<scalar or vector>> Parameter of timefreq function.
0041 %   winsize       : <<scalar>> Parameter of timefreq function.
0042 %   padratio      : <<scalar>> Parameter of timefreq function.
0043 %   ntimesout     : <<int>> Parameter of timefreq function.
0044 %   nfreqs        : <<int>> Parameter of timefreq function.
0045 %  ---
0046 % ---
0047 %
0048 % [output]
0049 % TimeFrequency file (.tf.mat) is created. This file has the following
0050 % variables:
0051 % data  : <<matrix>> Time-frequency data. F-by-T-by-I matrix, where F, T,
0052 %         I are number of frequency bins (=length(TFinfo.freq), time
0053 %         window length (=length(TFinfo.Tmsec) and the number of vertices
0054 %         (=length(TFinfo.ix_act)), respectively.
0055 % TFinfo: <<struct>> Information of time-frequency data.
0056 % --- fields of TFinfo
0057 %  tf_type   : <<string>> Set to 'ITPC'.
0058 %  Tsample   : <<vector>> Timepoint indices of the original data.
0059 %  Pretrigger: <<int>> Pretrigger period (timepoints).
0060 %  SampleFreq: <<float>> Sampling frequency of the original data.
0061 %  Tmsec     : <<vector>> Actual time of each timepoints in msec
0062 %              (length(Tsample) = length(Tmsec)).
0063 %  Ntrial    : <<vector>> Number of trials of each current data.
0064 %  freq      : <<vector>> Frequency vector [Hz]. Each row of
0065 %              'data(:,f,:)' has frequency TFinfo.freq(f).
0066 %  ix_act    : <<int vector>> Vertex indices on which time-frequency data
0067 %              is calculated.
0068 %  Wact      : <<double matrix>> Smoothing filter matrix. The size of
0069 %              smoothing filter is automatically determined from mean
0070 %              distance between vertices after reduction
0071 %              (tf_parm.reduce). This matrix will be used for displaying
0072 %              purpose.
0073 %  ix_act_ex : <<int vector>> Vertex indices corresponding to the row of
0074 %              the smoothing matrix TFinfo.Wact.
0075 % ---
0076 %
0077 % [example]
0078 % The following code is an example to calculate ITPC on template cortical
0079 % surface model.
0080 % >> tf_parm.curr_file = './subjects/0001/result/0001_LR.curr.mat';
0081 % >> tf_parm.brain_file = './subjects/0001/brain/0001.brain.mat';
0082 % >> tf_parm.key = 'fsaverage.sphere.reg';
0083 % >> tf_parm.tbrain_file ...
0084 %     = './subjects/fsaverage/brain/fsaverage.brain.mat';
0085 % >> tf_parm.reduce = 0.2;
0086 % >> tf_parm.fmin = 2.0;
0087 % >> tf_parm.tf_file_out ...
0088 %     = './subjects/fsaverage/result/fsaverage_LR_ITPC.tf.mat';
0089 % >> vb_job_itpc(proj_root,tf_parm);
0090 %
0091 % [history]
0092 % 2011-01-04 Taku Yoshioka
0093 % 2011-02-28 taku-y
0094 %  [enhancement] EEGLAB timefreq function supported.
0095 % 2011-07-01 taku-y
0096 %  [minor] vb_disp() was replaced by vb_disp_nonl() for displaying
0097 %  progress messages.
0098 %
0099 % Copyright (C) 2011, ATR All Rights Reserved.
0100 % License : New BSD License(see VBMEG_LICENSE.txt)
0101 
0102 vb_disp('--- ITPC calculation');
0103 
0104 %
0105 % Verbose level
0106 %
0107 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0108 const = vb_define_verbose;
0109 vb_struct2vars(const,{'VERBOSE_LEVEL_NONE','VERBOSE_LEVEL_EMERGENCY', ...
0110                     'VERBOSE_LEVEL_WARNING','VERBOSE_LEVEL_NOTICE', ...
0111                     'VERBOSE_LEVEL_INFO','VERBOSE_LEVEL_DEBUG'});
0112 VERBOSE_LEVEL_DEBUG = const.VERBOSE_LEVEL_DEBUG;
0113 
0114 %
0115 % Input arguments
0116 %
0117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0118 vb_struct2vars(tf_parm,{'curr_file','brain_file','key', ...
0119                     'tf_file_out','tbrain_file','reduce','fmin','fmax', ...
0120                     'elim_time','timefreq_parm'});
0121 
0122 if isempty(elim_time), elim_time = 0.05; end
0123 
0124 %
0125 % Preparation of cortical dipole indices
0126 %
0127 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0128 % Number of dipoles on which TF analysis conducted.
0129 if ~isempty(key), % use template brain
0130   % note: consistency check between tbrain_file and key{1} is required,
0131   % but postponed. (7 Dec. 2010, taku-y)
0132   brain_file_tmp = tbrain_file;
0133 else % single subject analysis
0134   brain_file_tmp = brain_file;
0135 end
0136 
0137 % Reduce vertex of the cortical model
0138 V = vb_load_cortex(brain_file_tmp);
0139 ix_act = vb_reduce_cortex(brain_file_tmp,1:size(V,1),reduce);
0140 I = length(ix_act);
0141 clear V;
0142 
0143 % For single subject analysis, find overlap between ix_act and
0144 % Jinfo.ix_act_ex. ix_act is replaced by it and I is overwritten.
0145 if isempty(key), % single subject analysis
0146   Jinfo = vb_load_current_info(curr_file);
0147   [tmp,ixx] = intersect(Jinfo.ix_act_ex,ix_act);
0148   ix_act = Jinfo.ix_act_ex(ixx);
0149   I = length(ix_act);
0150   Wact = Jinfo.Wact(ixx,:);
0151 end
0152 
0153 % Calculate mean distance between reduced vertices to determine smoothing
0154 % filter size
0155 R = 0.0;
0156 N = 0;
0157 [nextDD,nextIX] = vb_load_cortex_neighbour(brain_file_tmp);
0158 for i=1:length(ix_act)
0159   [tmp,ixx] = union(nextIX{ix_act(i)},ix_act);
0160   R = R+sum(nextDD{ix_act(i)}(ixx));
0161   N = N+length(ixx);
0162 end
0163 R = R/N;
0164 
0165 % Calculate smoothing filter. It will be used for displaying purpose, not
0166 % affect further analysis.
0167 V = vb_load_cortex(brain_file_tmp);
0168 [W,ix_act_ex] ...
0169     = vb_spatial_gauss_filter(brain_file_tmp,R,2*R,ix_act,1,-1,1);
0170 W = W./repmat(sum(W,2),[1 size(W,2)]);
0171 TFinfo.Wact = W;
0172 TFinfo.ix_act_ex = ix_act_ex;
0173 clear V;
0174 
0175 %
0176 % Preparation of variables time-frequency data stored
0177 %
0178 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0179 if isempty(timefreq_parm), % wavelet() parameters
0180   % Number of timepoints
0181   Jinfo = vb_load_current_info(curr_file);
0182   T  = length(Jinfo.Tmsec);    % time points
0183   DJ = 0.25;                   % scale of time points spacing
0184   DT = 1/(Jinfo.SampleFreq);   % original time points spacing
0185   S0 = 2*DT;                   % smallest freq. scale
0186 
0187   % Number of frequency bins
0188   [tmp,period,coi] = wavelet(randn(1,T),DT,1,DJ);
0189   TFinfo.freq      = 1./period;
0190   ix_freq          = find(TFinfo.freq>=fmin);
0191   TFinfo.freq      = TFinfo.freq(ix_freq);
0192   F                = length(TFinfo.freq);
0193 
0194   % Eliminate start and end time points
0195   t_elim  = ceil(T*elim_time);
0196   ix_time = (t_elim+1):(T-t_elim);
0197   Tmsec   = Jinfo.Tmsec(ix_time);
0198   T       = length(Tmsec);
0199 
0200   % Whole data
0201   data       = zeros(F,T,I);
0202 else % timefreq() parameters
0203   vb_struct2vars(timefreq_parm,{'cycles','winsize','padratio', ...
0204                       'ntimesout','nfreqs'});
0205   
0206   % Number of timepoints and frequency bins
0207   Jinfo  = vb_load_current_info(curr_file);
0208   frames = length(Jinfo.Tmsec);
0209   tlimits = [min(Jinfo.Tmsec) max(Jinfo.Tmsec)];
0210   srate   = Jinfo.SampleFreq;
0211   if isempty(fmax), 
0212     fmax = srate/2;
0213   end
0214   
0215   [tmp,freqs,Tmsec] ...
0216       = timefreq(randn(frames,1),srate,'cycles',cycles, ...
0217                  'freqs',[fmin fmax],'padratio',padratio, ...
0218                  'winsize',winsize,'tlimits',tlimits, ...
0219                  'ntimesout',ntimesout,'nfreqs',nfreqs, ...
0220                  'verbose','off');
0221   
0222   T                = length(Tmsec);
0223   ix_time          = 1:length(Tmsec);
0224   TFinfo.freq      = freqs;
0225   ix_freq          = find(TFinfo.freq>=fmin);
0226   TFinfo.freq      = TFinfo.freq(ix_freq);
0227   F                = length(TFinfo.freq);
0228   
0229   % Whole data
0230   data = zeros(F,T,I);
0231 end
0232 
0233 vb_disp(['Number of frequency bins       : ' num2str(F)]);
0234 vb_disp(['Number of timepoints           : ' num2str(T)]);
0235 vb_disp(['Number of vertices             : ' num2str(I)]);
0236 
0237 %
0238 % Main routine
0239 %
0240 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0241 Jinfo      = vb_load_current_info(curr_file);
0242 Ntrial     = Jinfo.Ntrial;
0243 
0244 if ~isempty(key), % use template brain
0245   C = vb_load_cortex_pmat(brain_file,key);
0246 
0247   if strcmp(Jinfo.curr_type,'Z-current'), 
0248     CW = sparse(C(ix_act,Jinfo.ix_act_ex)*Jinfo.Wact);
0249   else
0250     CW = sparse(C(ix_act,Jinfo.ix_act_ex));
0251   end
0252 else % single subject analysis
0253   CW = Wact;
0254 end
0255 
0256 prg = 0;
0257 prg_all = Ntrial;
0258 prg_str = sprintf('Current : %d%% processed',0);
0259 vb_disp_nonl(prg_str);
0260 
0261 for j=1:Ntrial % trial loop
0262   [tmp,Z] = vb_load_current(curr_file,Jinfo.curr_type,false,j);
0263     
0264   for k=1:I % vertex loop
0265     ix = find(abs(CW(k,:))>0);
0266     J = CW(k,ix)*Z(ix,:); % J-current at ix_act-th dipole
0267     
0268     if isempty(timefreq_parm), 
0269       wave = wavelet(J,DT,1,DJ);
0270     else
0271       wave = timefreq(J',srate,'cycles',cycles, ...
0272                       'freqs',[fmin fmax],'padratio',padratio, ...
0273                       'winsize',winsize,'tlimits',tlimits, ...
0274                       'ntimesout',ntimesout,'nfreqs',nfreqs, ...
0275                       'verbose','off');
0276     end
0277 
0278     wave = wave(ix_freq,ix_time);
0279     data(:,:,k) = data(:,:,k)+wave./abs(wave);
0280   end
0281   
0282   % progress message
0283   prg = prg+1;
0284   for ii=1:length(prg_str); vb_disp_nonl(sprintf('\b')); end
0285   prg_str = sprintf('Current : %d%% processed',ceil(100*(prg/prg_all)));
0286   vb_disp_nonl(prg_str);
0287 end
0288 
0289 for ii=1:length(prg_str); vb_disp_nonl(sprintf('\b')); end
0290 vb_disp('Current : finished');
0291   
0292 % Normalization
0293 data = abs(data)./Ntrial;
0294 
0295 %
0296 % Make struct 'TFinfo'
0297 %
0298 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0299 Jinfo = vb_load_current_info(curr_file);
0300 
0301 TFinfo.tf_type    = 'ERSP';
0302 TFinfo.Tmsec      = Tmsec;
0303 TFinfo.Ntrial     = Ntrial;
0304 TFinfo.ix_act     = ix_act;
0305 
0306 if isempty(timefreq_parm),
0307   TFinfo.Tsample    = Jinfo.Tsample(ix_time);
0308   TFinfo.Pretrigger = Jinfo.Pretrigger-t_elim;
0309   TFinfo.SampleFreq = Jinfo.SampleFreq;
0310   TFinfo.freq_scale = 'log'; % for wavelet()
0311 else
0312   Tmsec = TFinfo.Tmsec;
0313   [tmp,TFinfo.Pretrigger] = min(abs(Tmsec(:)));
0314   TFinfo.SampleFreq = ((Tmsec(end)-Tmsec(1))./1000)./length(Tmsec);
0315   TFinfo.freq_scale = 'linear';
0316 end
0317 
0318 %
0319 % Save result
0320 %
0321 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0322 vb_fsave(tf_file_out,'tf_parm','data','TFinfo');
0323 
0324 %
0325 % Save project file
0326 %
0327 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0328 proj_file = get_project_filename;
0329 if isempty(proj_file), 
0330   return;
0331 end
0332 
0333 project_file_mgr('load', proj_file);
0334 project_file_mgr('add', 'tf_parm', tf_parm);
0335 
0336 return;

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