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)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;