timef.m

含有多种ICA算法的eeglab工具箱

% timef() - Returns estimates and plots of mean event-related spectral
%           perturbation (ERSP) and inter-trial coherence (ITC) changes 
%           across event-related trials (epochs) of a single input time series. 
%
%         * Uses either fixed-window, zero-padded FFTs (fastest), wavelet
%           0-padded DFTs (both Hanning-tapered), OR multitaper spectra ('mtaper').
%         * For the wavelet and FFT methods, output frequency spacing 
%           is the lowest frequency ('srate'/'winsize') divided by 'padratio'.
%           NaN input values (such as returned by eventlock()) are ignored.
%         * If 'alpha' is given, then bootstrap statistics are computed 
%           (from a distribution of 'naccu' surrogate data trials) and 
%           non-significant features of the output plots are zeroed out 
%           (i.e., plotted in green). 
%         * Given a 'topovec' topo vector and 'elocs' electrode location file,
%           the figure also shows a topoplot() of the specified scalp map.
%         * Note: Left-click on subplots to view and zoom in separate windows.
% Usage: 
%        >> [ersp,itc,powbase,times,freqs,erspboot,itcboot] = ...
%                timef(data,frames,tlimits,srate,cycles,...
%                        'key1',value1,'key2',value2, ... );        
% NOTE:                                        
%        >> timef details  % scrolls more detailed information about timef
%
% Required inputs:     
%       data        = Single-channel data vector (1,frames*ntrials) (required)
%       frames      = Frames per trial                               {750}
%       tlimits     = [mintime maxtime] (ms) Epoch time limits    {[-1000 2000]}
%       srate       = data sampling rate (Hz)                        {250}
%       cycles      = If 0 -> Use FFTs (with constant window length) {0}
%                     If >0 -> Number of cycles in each analysis wavelet 
%                     If [wavecycles factor] -> wavelet cycles increase with frequency 
%                     beginning at wavecyles (0<factor<1; factor=1 -> no increase,
%                     standard wavelets; factor=0 -> fixed epoch length, as in FFT.
%                     Else, 'mtaper' -> multitaper decomposition 
%
%    Optional Inter-Irial Coherence (ITC) type:
%       'type'      = ['coher'|'phasecoher'] Compute either linear coherence 
%                      ('coher') or phase coherence ('phasecoher') also known
%                      as the phase coupling factor           {'phasecoher'}.
%
%    Optional detrending:
%       'detret'    = ['on'|'off'], Detrend data in time.       {'off'}
%       'detrep'    = ['on'|'off'], Detrend data across trials  {'off'}
%
%    Optional FFT/DFT parameters:
%       'winsize'   = If cycles==0: data subwindow length (fastest, 2^n<frames);
%                     If cycles >0: *longest* window length to use. This
%                      determines the lowest output frequency       {~frames/8}
%       'timesout'  = Number of output times (int<frames-winframes)       {200}
%       'padratio'  = FFT-length/winframes (2^k)                          {2}
%                      Multiplies the number of output frequencies by
%                      dividing their spacing. When cycles==0, frequency
%                      spacing is (low_freq/padratio).
%       'maxfreq'   = Maximum frequency (Hz) to plot (& to output, if cycles>0) 
%                      If cycles==0, all FFT frequencies are output.      {50}
%       'baseline'  = Spectral baseline end-time (in ms).                 {0}
%       'powbase'   = Baseline spectrum to log-subtract. {def|NaN->from data}
%
%    Optional multitaper parameters:
%       'mtaper'    = If [N W], performs multitaper decomposition. 
%                      (N is the time resolution and W the frequency resolution; 
%                      maximum taper number is 2NW-1). Overwrites 'winsize' and 'padratio'. 
%                     If [N W K], forces the use of K Slepian tapers (if possible).
%                      Phase is calculated using standard methods.
%                      The use of mutitaper with wavelets (cycles>0) is not 
%                      recommended (as multiwavelets are not implemented). 
%                      Uses Matlab functions DPSS, PMTM.      {no multitaper}
%
%    Optional bootstrap parameters:
%       'alpha'     = If non-0, compute two-tailed bootstrap significance prob. 
%                     level. Show non-signif. output values in green       {0}
%       'naccu'     = Number of bootstrap replications to accumulate       {200}
%       'baseboot'  = Bootstrap baseline to subtract (1 -> use 'baseline'(see above)
%                                                     0 -> use whole trial) {1}
%    Optional scalp map:
%       'topovec'   = Scalp topography (map) to plot                       {none}
%       'elocs'     = Electrode location file for scalp map   
%                     File should be ascii in format of  >> topoplot example   
%                     May also be an EEG.chanlocs struct. {file named in icadefs.m}
%    Optional plotting parameters:
%       'ploterps'  = ['on'|'off'] Plot power spectral perturbations       {'on'} 
%       'plotitc'   = ['on'|'off'] Plot inter trial coherence              {'on'}
%       'plotphase' = ['on'|'off'] Plot sign of the phase in the ITC panel, i.e.
%                     green->red, pos.-phase ITC, green->blue, neg.-phase ITC {'on'}
%       'erspmax'   = [real dB] set the ERSP max. for the scale (min= -max){auto}
%       'itcmax'    = [real<=1] set the ITC maximum for the scale          {auto}
%       'title'     = Optional figure title                                {none}
%       'marktimes' = Non-0 times to mark with a dotted vertical line (ms) {none}
%       'linewidth' = Line width for 'marktimes' traces (thick=2, thin=1)  {2}
%       'pboot'     = Bootstrap power limits (e.g., from timef())    {from data}
%       'rboot'     = Bootstrap ITC limits (e.g., from timef())      {from data}
%       'axesfont'  = Axes text font size                                  {10}
%       'titlefont' = Title text font size                                 {8}
%       'vert'      = [times_vector] -> plot vertical dashed lines at given times in ms.
%       'verbose'   = ['on'|'off'] print text                              {'on'}
% Outputs: 
%            ersp   = Matrix (nfreqs,timesout) of log spectral diffs. from baseline (dB) 
%            itc    = Matrix of inter-trial coherencies (nfreqs,timesout) (range: [0 1])
%          powbase  = Baseline power spectrum (subtracted from each window to compute 
%                     the ERSP).
%            times  = Vector of output times (sub-window centers) (in ms).
%            freqs  = Vector of frequency bin centers (in Hz).
%         erspboot  = Matrix (2,nfreqs) of [lower;upper] ERSP significance diffs.
%          itcboot  = Matrix (2,nfreqs) of [lower;upper] ITC thresholds (not diffs).
%
% Plot description:
%   Assuming both 'plotersp' and 'plotitc' options are 'on' (= default). The upper panel
%   presents the data ERSP (Event-Related Spectral Perturbation) in dB, with mean baseline 
%   spectral activity (in dB) subtracted. Use "'baseline', NaN" to prevent timef() from 
%   removing the baseline. The lower panel presents the data ITC (Inter-Trial Coherence). 
%   Click on any plot axes to pop up a new window (using 'axcopy()')
%   -- Upper left marginal panel presents the mean spectrum during the baseline period
%      (blue), and when significance is set, the significance threshold at each frequency
%      (dotted green-black trace).
%   -- The marginal panel under the ERSP image shows the maximum (green) and minimum 
%      (blue) ERSP values relative to baseline power at each frequency.
%   -- The lower left marginal panel shows mean ITC across the imaged time range (blue),
%      and when significance is set, the significance threshold (dotted green-black).
%   -- The marginal panel under the ITC image shows the ERP (which is produced by ITC 
%      across the data spectral pass band).
%
% Author: Sigurd Enghoff, Arnaud Delorme & Scott Makeig
%          CNL / Salk Institute 1998- | SCCN/INC, UCSD 2002-
%
% Known problems:
%   Significance masking currently fails for linear coherence.
%
% See also: crossf()
 
%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) 1998 Sigurd Enghoff, Scott Makeig, Arnaud Delorme, 
% CNL / Salk Institute 8/1/98-8/28/01
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% $Log: timef.m,v $
% Revision 1.73  2004/03/25 15:12:08  scott
% help message re EEG.chanlocs struct possibility
%
% Revision 1.72  2004/03/09 17:35:33  arno
% msg
%
% Revision 1.71  2004/03/01 16:17:49  arno
% default baseboot 1
%
% Revision 1.70  2004/03/01 02:25:06  arno
% change help message
%
% Revision 1.69  2004/02/09 18:01:27  arno
% debug multitaper
%
% Revision 1.68  2004/01/23 02:55:15  scott
% documenting maxersp
%
% Revision 1.67  2004/01/23 02:23:35  scott
% add erspmax
%
% Revision 1.66  2004/01/06 17:00:37  arno
% header typo
%
% Revision 1.65  2003/12/03 02:32:41  arno
% verbose option
%
% Revision 1.64  2003/10/22 17:37:03  arno
% slight imprecision on frequency for FFT when padratio > 1
%
% Revision 1.63  2003/08/27 21:49:11  arno
% header typo
%
% Revision 1.62  2003/08/25 21:58:50  arno
% header typo
%
% Revision 1.61  2003/08/05 22:01:52  scott
% header msg edits
%
% Revision 1.60  2003/08/04 16:38:36  arno
% plot description and curve color
%
% Revision 1.59  2003/08/04 14:40:51  arno
% description of plot
%
% Revision 1.58  2003/08/02 23:07:10  arno
% debug h(9)
%
% Revision 1.57  2003/07/23 00:28:22  scott
% help msg
%
% Revision 1.56  2003/05/12 23:41:54  arno
% text typo
%
% Revision 1.55  2003/05/05 16:27:11  arno
% debug FFT scale
%
% Revision 1.54  2003/03/14 01:04:12  arno
% typo header
%
% Revision 1.53  2003/03/12 20:09:28  scott
% header edits -sm
%
% Revision 1.52  2003/01/08 23:37:52  arno
% add discalimer about linear coherence bug
%
% Revision 1.51  2002/11/15 03:04:07  arno
% header for web
%
% Revision 1.50  2002/10/18 16:49:53  arno
% win' -> transpose(win)
%
% Revision 1.49  2002/08/14 21:07:05  arno
% hanning debug
%
% Revision 1.48  2002/08/14 21:02:19  arno
% implementing hanning funciton
%
% Revision 1.47  2002/08/12 01:47:29  arno
% color
%
% Revision 1.46  2002/08/11 22:30:20  arno
% color
%
% Revision 1.45  2002/08/09 22:29:44  arno
% implementing wavelet factor
%
% Revision 1.44  2002/07/22 14:28:56  arno
% debugging input baseline spectrum
%
% Revision 1.43  2002/07/11 18:19:30  arno
% header typo
%
% Revision 1.42  2002/07/11 15:27:28  arno
% debuging linear coherence
%
% Revision 1.41  2002/07/11 15:18:40  arno
% programing phase coherence 2
%
% Revision 1.40  2002/07/11 00:04:57  arno
% same
%
% Revision 1.39  2002/07/11 00:03:09  arno
% debugging itcmax
%
% Revision 1.38  2002/07/10 23:59:57  arno
% implement itcmax
%
% Revision 1.37  2002/07/10 21:46:03  arno
% adding phase argument
%
% Revision 1.36  2002/05/03 03:17:45  arno
% diffdlt 0.5 -> 1
%
% Revision 1.35  2002/05/01 22:25:21  arno
% no modif
%
% Revision 1.34  2002/05/01 21:23:24  arno
% no changes
%
% Revision 1.33  2002/04/30 04:27:54  arno
% trying to debug phsamp, still unsucessfull
%
% Revision 1.32  2002/04/29 15:13:50  scott
% debugging PA -sm
%
% Revision 1.31  2002/04/29 15:06:51  scott
% same -sm
%
% Revision 1.30  2002/04/29 15:03:31  scott
% same -sm
%
% Revision 1.29  2002/04/29 15:01:58  scott
% debugging -sm
%
% Revision 1.28  2002/04/29 14:58:56  scott
% debugging -sm
%
% Revision 1.27  2002/04/29 14:32:56  scott
% removing debugging statements -sm
%
% Revision 1.26  2002/04/29 14:32:31  scott
% adding ; -sm
%
% Revision 1.25  2002/04/29 14:29:53  scott
% debugging cumulX/PA -sm
%
% Revision 1.24  2002/04/29 14:24:13  scott
% debugging -sm
%
% Revision 1.23  2002/04/29 14:21:23  scott
% same -sm
%
% Revision 1.22  2002/04/29 14:15:58  scott
% insured cumulX sized -sm
%
% Revision 1.21  2002/04/29 14:12:33  scott
% used switches to test g.phsamp -sm
%
% Revision 1.20  2002/04/29 14:07:08  scott
% fixed typo -sm
%
% Revision 1.19  2002/04/29 14:02:34  scott
% made sure cumulX is computed -sm
%
% Revision 1.18  2002/04/29 13:57:49  scott
% modified PC->PA, 'phasecouple'->'phsamp', made output format (phs,amp,time) -sm
%
% Revision 1.17  2002/04/27 21:26:24  scott
% debugging PC -sm
%
% Revision 1.16  2002/04/27 21:19:02  scott
% debugging PC -sm
%
% Revision 1.15  2002/04/27 21:17:27  scott
% debugging PC -sm
%
% Revision 1.14  2002/04/27 21:13:57  scott
% added undocumented arg 'phasecouple',{'on'|'off'} -sm
%
% Revision 1.13  2002/04/25 02:56:03  arno
% redebugging topovec
%
% Revision 1.12  2002/04/25 02:54:33  arno
% improved topovec check
%
% Revision 1.11  2002/04/23 18:34:29  arno
% modified baseline way of computation
%
% Revision 1.10  2002/04/23 18:28:02  arno
% correcting coher computation
%
% Revision 1.9  2002/04/11 19:56:12  arno
% debuging baseboot -ad & lf
%
% Revision 1.8  2002/04/11 02:10:27  arno
% correcting typo
%
% Revision 1.7  2002/04/09 00:24:55  arno
% editing latest modifications
%
% Revision 1.6  2002/04/09 00:21:17  arno
% adding vertical bars
%
% Revision 1.5  2002/04/08 19:57:19  arno
% Editing, maxfreq-> Niquist
%
% Revision 1.4  2002/04/07 02:49:35  scott
% clarified hlpe message, changed default srate and winsize -sm
%
% Revision 1.3  2002/04/06 03:48:07  arno
% changing input for 1 channel for topoplot
%
% Revision 1.2  2002/04/06 03:40:58  arno
% modifying location file check
%
% Revision 1.1  2002/04/05 17:36:45  jorn
% Initial revision
%

% 10-19-98 avoided division by zero (using MIN_ABS) -sm
% 10-19-98 improved usage message and commandline info printing -sm
% 10-19-98 made valid [] values for tvec and g.elocs -sm
% 04-01-99 added missing freq in freqs and plots, fixed log scaling bug -se & -tpj
% 06-29-99 fixed frequency indexing for constant-Q -se
% 08-24-99 reworked to handle NaN input values -sm
% 12-07-99 adjusted ERPtimes to plot ERP under ITC -sm
% 12-22-99 debugged ERPtimes, added BASE_BOOT -sm 
% 01-10-00 debugged BASE_BOOT=0 -sm
% 02-28-00 added NOTE on formula derivation below -sm
% 03-16-00 added axcopy() feature -sm & tpj
% 04-16-00 added multiple marktimes loop -sm
% 04-20-00 fixed ITC cbar limits when spcified in input -sm
% 07-29-00 changed frequencies displayed msg -sm
% 10-12-00 fixed bug in freqs when cycles>0 -sm
% 02-07-01 fixed inconsistency in BASE_BOOT use -sm
% 08-28-01 matlab 'key' value arguments -ad
% 08-28-01 multitaper decomposition -ad
% 01-25-02 reformated help & license -ad 
% 03-08-02 debug & compare to old timef function -ad 
% 03-16-02 timeout automatically adjusted if too high -ad 
% 04-02-02 added 'coher' option -ad 

function [P,R,mbase,times,freqs,Pboot,Rboot,PA] = timef( X, frame, tlimits, Fs, varwin, varargin);

% Note: PA is output of 'phsamp','on' 

%varwin,winsize,g.timesout,g.padratio,g.maxfreq,g.topovec,g.elocs,g.alpha,g.marktimes,g.powbase,g.pboot,g.rboot)

% ITC:   Normally, R = |Sum(Pxy)| / (Sum(|Pxx|)*Sum(|Pyy|)) is coherence.
%        But here, we consider    Phase(Pyy) = 0 and |Pyy| = 1 -> Pxy = Pxx
%        Giving, R = |Sum(Pxx)|/Sum(|Pxx|), the inter-trial coherence (ITC)
%        Also called 'phase-locking factor' by Tallon-Baudry et al. (1996)

% Constants set here:
ERSP_CAXIS_LIMIT = 0;           % 0 -> use data limits; else positive value
                                % giving symmetric +/- caxis limits.
ITC_CAXIS_LIMIT  = 0;           % 0 -> use data limits; else positive value
                                % giving symmetric +/- caxis limits.
MIN_ABS          = 1e-8;        % avoid division by ~zero

% Commandline arg defaults:
DEFAULT_EPOCH	= 750;			% Frames per trial
DEFAULT_TIMLIM = [-1000 2000];	% Time range of g.frames (ms)
DEFAULT_FS	= 250;			% Sampling frequency (Hz)
DEFAULT_NWIN	= 200;			% Number of windows = horizontal resolution
DEFAULT_VARWIN	= 0;			% Fixed window length or fixed number of cycles.
								% =0: fix window length to that determined by nwin