crossf.m

来自「含有多种ICA算法的eeglab工具箱」· M 代码 · 共 1,560 行 · 第 1/4 页

% crossf() - Returns estimates and plots event-related coherence (ERCOH) 
%        between two input data time series. A lower panel (optionally) shows 
%        the coherence phase difference between the processes. In this panel:
%           -90 degrees (blue)   means x leads y by a quarter cycle.
%            90 degrees (orange) means y leads x by a quarter cycle.
%        Click on any subplot to view separately and zoom in/out.
%
% Function description:
%        Uses EITHER fixed-window, zero-padded FFTs (fastest) OR constant-Q 
%        0-padded wavelet DFTs (more even sensitivity across frequencies), 
%        both Hanning-tapered.  Output frequency spacing is the lowest 
%        frequency ('srate'/'winsize') divided by the 'padratio'.
%
%        If an 'alpha' value is given, then bootstrap statistics are 
%        computed (from a distribution of 'naccu' (200) surrogate baseline
%        data epochs) for the baseline epoch, and non-significant features 
%        of the output plots are zeroed (and shown in green). The baseline
%        epoch is all windows with center latencies < the given 'baseline' value 
%        or, if 'baseboot' is 1, the whole epoch. 
% Usage: 
%        >> [coh,mcoh,timesout,freqsout,cohboot,cohangles] ...
%                       = crossf(x,y,frames,tlimits,srate,cycles, ...
%                                        'key1', 'val1', 'key2', val2' ...);
%
% Required inputs:
%       x       = first single-channel data set (1,frames*nepochs)      
%       y       = second single-channel data set (1,frames*nepochs)     
%       frames  = frames per epoch                                   {750}
%       tlimits = [mintime maxtime] (ms) epoch latency limits  {[-1000 2000]}
%       srate   = data sampling rate (Hz)                            {250}
%       cycles  = 0  -> Use FFTs (with constant window length) 
%               = >0 -> Number of cycles in each analysis wavelet 
%               = [cycles expfactor] -> if 0 < expfactor < 1,  the number 
%                 of wavelet cycles expands with frequency from cycles
%                 If expfactor = 1, no expansion; if = 0, constant
%                 window length (as in FFT)            {default cycles: 0}
%
%    Optional Coherence Type:
%       'type'  = ['coher'|'phasecoher'] Compute either linear coherence
%                 ('coher') or phase coherence ('phasecoher') also known
%                 as phase coupling factor' {default: 'phasecoher'}.
%       'subitc' = ['on'|'off'] subtract stimulus locked Inter-Trial Coherence 
%                 from x and y. This computes the  'intrinsic' coherence
%                 x and y not arising from common synchronization to 
%                 experimental events. See notes. {default: 'off'}
%       'shuffle' = integer indicating the number of estimates to compute
%                 bootstrap coherence based on shuffled trials. This estimates
%                 the coherence arising only from time locking of x and y
%                 to experimental events (opposite of 'subitc'). {default: 0}. 
%
%    Optional Detrend:
%       'detret' = ['on'|'off'], Linearly detrend data within epochs. {'off'}
%       'detrep' = ['on'|'off'], Linearly detrend data across trials  {'off'}
%
%    Optional FFT/DFT:
%       '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 latencies (int<frames-winsize) {def: 200}
%       'padratio' = FFTlength/winsize (2^k)                     {def: 2}
%                    Multiplies the number of output frequencies by
%                    dividing their spacing. When cycles==0, frequency
%                    spacing is (low_frequency/padratio).
%       'maxfreq'  = Maximum frequency (Hz) to plot (& output if cycles>0) 
%                    If cycles==0, all FFT frequencies are output.{def: 50}
%       'baseline' = Coherence baseline end latency (ms). NaN=no baseline  {NaN}
%       'powbase'  = Baseline spectrum to log-subtract.  {default: from data}
%
%    Optional Bootstrap:
%       'alpha'    = If non-0, compute two-tailed bootstrap significance prob.
%                    level. Show non-signif output values as green. {0}
%       'naccu'    = Number of bootstrap replications to compute {200}
%       'boottype' = ['times'|'timestrials'] Bootstrap type: Either shuffle
%                    windows ('times') or windows and trials ('timestrials')
%                    Option 'timestrials' requires more memory {default 'times'}
%       'memory'   = ['low'|'high'] 'low' -> decrease memory use {default: 'high'}
%       'baseboot' = Extent of bootstrap shuffling (0=to 'baseline'; 1=whole epoch). 
%                    If no baseline is given (NaN), extent of bootstrap shuffling 
%                    is the whole epoch                         {default: 0}
%       'rboot'    = Input bootstrap coherence limits (e.g., from crossf()) 
%                    The bootstrap type should be identical to that used
%                    to obtain the input limits. {default: compute from data}
% Optional Scalp Map:
%       'topovec'  = (2,nchans) matrix, plot scalp maps to plot {[]}
%                    ELSE (c1,c2), plot two cartoons showing channel locations.
%       'elocs'    = Electrode location file for scalp map       {none}
%                    File should be ascii in format of >> topoplot example   
%
% Optional Plot and Compute Features:
%       'compute'   = ['matlab'|'C'] Use C sub-routine to speed up the
%                     computation                      {default 'matlab'}
%       'plotamp'   = ['on'|'off'], Plot coherence magnitude       {'on'}
%       'maxamp'    = [real] Set the maximum for the amp. scale    {auto}
%       'plotphase' = ['on'|'off'], Plot coherence phase angle     {'on'}
%       'angleunit' = Phase units: 'ms' for msec or 'deg' for degrees {'deg'}
%       'title'     = Optional figure title {none}
%       'vert'      = Latencies to mark with a dotted vertical line   {none}
%       'linewidth' = Line width for marktimes traces (thick=2, thin=1) {2}
%       'cmax'      = Maximum amplitude for color scale  { use data limits }
%       'axesfont'  = Axes font size                               {10}
%       'titlefont' = Title font size                              {8}
%
% Outputs: 
%       coh         = Matrix (nfreqs,timesout) of coherence magnitudes 
%       mcoh        = Vector of mean baseline coherence at each frequency
%       timesout    = Vector of output latencies (window centers) (ms).
%       freqsout    = Vector of frequency bin centers (Hz).
%       cohboot     = Matrix (nfreqs,2) of [lower;upper] coher signif. limits
%                     if 'boottype' is 'trials',  (nfreqs,timesout, 2)
%       cohangle    = (nfreqs,timesout) matrix of coherence angles 
%
% Plot description:
%   Assuming both 'plotamp' and 'plotphase' options are 'on' (=default), the upper panel
%   presents the magnitude of either phase coherence or linear coherence, depending on 
%   the 'type' parameter (above). The lower panel presents the coherence phase difference 
%   (in degrees). Click on any plot to pop up a new window (using 'axcopy()').
%   -- The upper left marginal panel shows mean coherence during the baseline period
%      (blue), and when significance is set, the significance threshold (dotted black-green).
%   -- The horizontal panel under the coherence magnitude image indicates the maximum 
%      (green) and minimum (blue) coherence values across all frequencies. When significance 
%      is set (using option 'trials' for 'boottype'), an additional curve indicates the 
%      significance threshold (dotted black-green).
%
% Notes: 1) When cycles==0, nfreqs is total number of FFT frequencies.
%        2) 'blue' coherence lag -> x leads y; 'red' -> y leads x
%        3) The 'boottype' should be ideally 'timesframes', but this creates high 
%           memory demands, so the 'times' method must be used in many cases.
%        4) If 'boottype' is 'trials', the average of the complex bootstrap
%           is subtracted from the coherence to compensate for phase differences 
%           (the average is also subtracted from the bootstrap distribution). 
%           For other bootstraps, this is not necessary since the phase is random.
%        5) If baseline is non-NaN, the baseline is subtracted from
%           the complex coherence. On the left hand side of the coherence
%           amplitude image, the baseline is displayed as a magenta line
%           (if no baseline is selected, this curve represents the average
%           coherence at every given frequency).
%        6) If a out-of-memory error occurs, set the 'memory' option to 'low'
%           (Makes computation time slower; Only the 'times' bootstrap method 
%           can be used in this mode).
%
% Authors: Arnaud Delorme, Sigurd Enghoff & Scott Makeig
%          CNL/Salk Institute 1998-2001; SCCN/INC/UCSD, La Jolla, 2002-
%
% See also: timef()

% NOTE: one hidden parameter 'savecoher', 0 or 1

% Copyright (C) 8/1/98  Arnaud Delorme, Sigurd Enghoff & Scott Makeig, SCCN/INC/UCSD
%
% 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: crossf.m,v $
% Revision 1.63  2003/11/26 18:18:52  scott
% help msg
%
% Revision 1.62  2003/08/05 22:06:49  scott
% header edits
%
% Revision 1.61  2003/08/04 16:22:13  arno
% marginal plot color
%
% Revision 1.60  2003/05/05 16:31:18  arno
% debug last
%
% Revision 1.59  2003/05/05 16:25:24  arno
% debuging FFT freq scale
%
% Revision 1.58  2003/01/31 23:57:08  arno
% debuging plotting title
%
% Revision 1.57  2003/01/06 18:12:15  arno
% shuffle along the 3rd dimension
%
% Revision 1.56  2003/01/02 22:40:43  arno
% updateing linear coherence formulas
%
% Revision 1.55  2003/01/02 19:33:55  arno
% removing reverse angle for wavelet
%
% Revision 1.54  2002/11/15 01:54:52  arno
% header for web
%
% Revision 1.53  2002/11/15 01:45:06  scott
% can not -> cannot
%
% Revision 1.52  2002/10/24 15:08:12  arno
% debugging plot scalp map for channel
%
% Revision 1.51  2002/10/18 16:50:40  arno
% win' -> transpose(win)
%
% Revision 1.50  2002/09/26 00:40:48  arno
% updating header
%
% Revision 1.49  2002/09/12 03:22:44  scott
% help msg -sm
%
% Revision 1.48  2002/08/14 21:06:33  arno
% hanning debug
%
% Revision 1.47  2002/08/14 21:02:54  arno
% implementing hanning funciton
%
% Revision 1.46  2002/08/12 01:48:01  arno
% color
%
% Revision 1.45  2002/08/11 22:30:19  arno
% color
%
% Revision 1.44  2002/08/11 22:28:36  arno
% updating cycles
%
% Revision 1.43  2002/08/09 22:37:32  arno
% debugging cyclefact
%
% Revision 1.42  2002/08/09 22:29:42  arno
% implementing wavelet factor
%
% Revision 1.41  2002/07/26 00:27:56  arno
% significance for all plots
%
% Revision 1.40  2002/07/25 23:05:15  arno
% implementing condition comparison
%
% Revision 1.39  2002/07/16 23:08:25  arno
% correcting freqs for wavelet output
%
% Revision 1.38  2002/07/16 23:00:03  arno
% testing
%
% Revision 1.37  2002/07/16 04:51:07  arno
% *** empty log message ***
%
% Revision 1.36  2002/07/12 23:11:03  arno
% dftfilt->1
%
% Revision 1.35  2002/07/11 20:56:28  arno
% debugging maxamp
%
% Revision 1.34  2002/07/11 18:22:07  arno
% same
%
% Revision 1.33  2002/07/11 18:20:21  arno
% maxamp debug
%
% Revision 1.32  2002/07/11 18:16:00  arno
% implmenting maxamp
%
% Revision 1.31  2002/07/11 17:01:30  arno
% implementing phase coherence 2, speeding up 'coher' and 'phasecoher'
%
% Revision 1.30  2002/07/10 19:12:56  arno
% removing warning messages
%
% Revision 1.29  2002/07/03 23:45:03  arno
% correcting rboot output
%
% Revision 1.28  2002/07/02 21:51:41  scott
% edited help message -sm
%
% Revision 1.27  2002/06/28 17:57:10  arno
% returning coherence magnitudes
%
% Revision 1.26  2002/06/25 22:23:37  arno
% NEW MODULAR VERSION WITH SUBITC
%
% Revision 1.24  2002/05/23 23:33:35  scott
% remove subtraction of mean coherence -sm
%
% Revision 1.23  2002/04/25 02:54:07  arno
% debugging topovec
%
% Revision 1.22  2002/04/25 02:38:12  arno
% change significance (only one sided
%
% Revision 1.21  2002/04/25 02:18:19  arno
% debugging topovec
%
% Revision 1.20  2002/04/24 22:04:49  arno
% debugging error check
%
% Revision 1.19  2002/04/24 21:59:11  arno
% debugging
%
% Revision 1.18  2002/04/24 21:45:18  scott
% topovec bug -sm
%
% Revision 1.17  2002/04/24 21:43:00  scott
% editing topovec code -sm
%
% Revision 1.16  2002/04/24 21:02:28  scott
% added topoplots of two heads -sm
%
% Revision 1.15  2002/04/24 02:43:18  arno
% debugging amplitude coherence
%
% Revision 1.14  2002/04/20 00:53:14  arno
% restorings some outputs options
%
% Revision 1.13  2002/04/19 23:20:23  arno
% changing trial bootstrap, not optimal, waiting for further inputs
%
% Revision 1.12  2002/04/19 19:46:28  arno
% crossf with new trial coherence bootstrap (minus mean)
%
% Revision 1.11  2002/04/12 18:10:55  scott
% added note
%
% Revision 1.10  2002/04/12 01:30:43  arno
% compatibility for returning frequencies with timef
%
% Revision 1.9  2002/04/12 01:13:40  arno
% debuging no ploting option
%
% Revision 1.8  2002/04/12 01:08:13  arno
% change plotamps to plotamp in help message
%
% Revision 1.7  2002/04/12 00:41:37  arno
% programming baseboot
%
% Revision 1.6  2002/04/11 02:39:34  arno
% updated header message
%
% Revision 1.5  2002/04/10 01:29:45  arno
% adding vert optional input
%
% Revision 1.4  2002/04/09 19:36:38  arno
% corrected bootstrap optional input
%
% Revision 1.3  2002/04/09 18:59:06  arno
% corrected typo in header that made the function to crash
%
% Revision 1.2  2002/04/07 02:24:36  scott
% worked on hlpe message, changed some defaults -sm
%
% Revision 1.1  2002/04/05 17:36:45  jorn
% Initial revision
%

% 11-20-98 defined g.linewidth constant -sm
% 04-01-99 made number of frequencies consistent -se
% 06-29-99 fixed constant-Q freq indexing -se
% 08-13-99 added cohangle plotting -sm
% 08-20-99 made bootstrap more efficient -sm
% 08-24-99 allow nan values introduced by possible eventlock() preproc. -sm
% 03-16-00 added lead/lag interpretation to help msg - sm & eric visser
% 03-16-00 added axcopy() feature -sm & tpj
% 04-20-00 fixed Rangle sign for wavelets, added verts array -sm
% 01-22-01 corrected help msg when nargin<2 -sm & arno delorme
% 01-25-02 reformated help & license, added links -ad 
% 03-09-02 function restructuration -ad
%  add 'key', val arguments (+ external baseboot, baseline, color axis, angleunit...)
%  add detrending (across time and trials) + 'coher' option for amplitude coherence
%  significance only if alpha is given, ploting options in 'plotamp' and 'plotphase'
% 03-16-02 timeout automatically adjusted if too high -ad 
% 04-03-02 added new options for bootstrap -ad 

% There are 3 "objects" Tf, Coher and Boot which are handled
% - by specific functions under Matlab
%    (Tf) function Tf = tfinit(...) - create object Time Frequency (Tf) associated with some data
%    (Tf) function [Tf, itcvals] = tfitc(...) - compute itc for the selected data
%    (Tf) function [Tf, itcvals] = tfitcpost(Tf, trials) - itc normlisation 
%    (Tf) function [Tf, tmpX] = tfcomp(Tf, trials, times) - compute time freq. decomposition
%    (Coher) function Coher = coherinit(...) - initialize coherence object
%    (Coher) function Coher = cohercomp(Coher, tmpX, tmpY, trial, time) - compute coherence
%    (Coher) function Coher = cohercomppost(Coher, trials) - coherence normalization
%    (Boot) function Boot = bootinit(...) - intialize bootstrap object
%    (Boot) function Boot = bootcomp(...) - compute bootstrap
%    (Boot) function [Boot, Rbootout] = bootcomppost(...) - bootstrap normalization
% - by real objects under C++ (see C++ code)

function [R,mbase,times,freqs,Rbootout,Rangle, trialcoher, Tfx, Tfy] = crossf(X, Y, frame, tlimits, Fs, varwin, varargin)

%varwin,winsize,nwin,oversmp,maxfreq,alpha,verts,caxmax)

% Commandline arg defaults:
DEFAULT_ANGLEUNITS = 'deg';     % angle plotting units - 'ms' or 'deg'
DEFAULT_EPOCH	= 750;			% Frames per epoch