crossfold.m

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

% crossf() - Returns estimates and plot of event-related coherence (ERC) changes
%            between data from two input channels. The lower panel gives the
%            coherent phase difference between the processes. In this panel, for Ex.
%               -90 degrees (blue) means xdata leads ydata by a quarter cycle.
%                90 degrees (orange) means ydata leads xdata by a quarter cycle.
%            Click on each subplot to view separately and zoom in/out.
%
% Function description:
%            Uses EITHER fixed-window, zero-padded FFTs (faster) OR constant-Q 
%            0-padded DFTs (better sensitivity), both Hanning-tapered. Output 
%            frequency spacing is the lowest frequency (srate/winsize) divided 
%            by the padratio.
%
%            If number of output arguments > 4, then bootstrap statistics are 
%            computed (from a distribution of 200 (NACCU) surrogate baseline
%            data epochs) for the baseline epoch, and non-significant features 
%            of the output plots are zeroed (e.g., plotted in green). Baseline
%            epoch is all windows with center times < 0 (MAX_BASELN)
%
%            If number of output arguments > 5, coherency angles (lags) at
%            significant coherency (time,frequency) points are plotted as well.
%
% Usage: 
%      >> [coh,mcoh,timesout,freqsout,cohboot,cohangles] = crossf(xdata,ydata,...
%                                              frames,tlimits,titl,          ...
%                                              srate,cycles,winsize,timesout,...
%                                              padratio,maxfreq,alpha,verts);
%
% Inputs:
%       xdata       = first single-channel (1,frames*nepochs) data  {none}
%       ydata       = second single-channel (1,frames*nepochs) data {none}
%       frames      = frames per epoch                        {768}
%       tlimits     = epoch time limits (ms) [mintime maxtime]{-1000 2000}
%       titl        = figure title                            {none}
%       srate       = data sampling rate (Hz)                 {256}
%       cycles      = >0 -> number of cycles in each analysis window (slower)
%                     =0 -> use FFT (constant window length)  {0}
%       winsize     = cycles==0: data subwindow length (2^k<frames)
%                     cycles >0: *longest* window length to use; 
%                     determines the lowest output frequency  {~frames/8}
%       timesout    = number of output times (int<frames-winsize){200}
%       padratio    = FFT-length/winsize (2^k)                {2}
%                     Multiplies the number of output frequencies.
%       maxfreq     = maximum frequency to plot (Hz)          {50}
%       alpha       = Two-tailed bootstrap signif. probability {0.02}
%                     Sets n.s. plotted output values to green (0). 
%                     NOTE that it requires at least FIVE output arguments!
%       verts       = times of vertical lines (other than time 0) {none}
%       caxma       = color axis maximum (magnitude) {default: from data}
%
% Outputs: 
%       coh         = between-channel coherency changes (nfreqs,timesout)
%       mcoh        = vector of mean baseline coherence at each frequency
%       timesout    = vector of output times (subwindow centers) in ms.
%       freqsout    = vector of frequency bin centers in Hz.
%       cohboot     = [2,nfreqs] matrix of [lower;upper] coh significance diffs.
%       cohangle    = coherency angles (nfreqs,timesout) 
%
% Note: when cycles==0, nfreqs is total number of FFT frequencies.
%
% Authors: Sigurd Enghoff & Scott Makeig, SCCN/INC/UCSD, La Jolla, 1998 
%
% See also: timef()

% Copyright (C) 8/1/98 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

% 11-20-98 defined 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 

function [R,mbase,times,freqs,Rboot,Rangle,Rsignif] = crossf(X,Y,epoch,timelim,ftitle,Fs,varwin,winsize,nwin,oversmp,maxfreq,alpha,verts,caxmax)

% Constants set here:
MAX_BASELN      = 0;            % Windows with center times < this are in baseline.
NACCU           = 200;			% Number of sub-windows to accumulate
if nargin>13
  COH_CAXIS_LIMIT = caxmax;
else
  COH_CAXIS_LIMIT = 0;          % 0 -> use data limits; else positive value
end
                                % giving symmetric +/- caxis limits.
AXES_FONT       = 10;
LINEWIDTH       = 2;
TITLE_FONT      = 8;
ANGLEUNITS      = 'deg';        % angle plotting units - 'ms' or 'deg'

% Commandline arg defaults:
DEFAULT_EPOCH	= 768;			% Frames per epoch
DEFAULT_TIMELIM = [-1000 2000];	% Time range of epochs (ms)
DEFAULT_FS		= 256;			% Sampling frequency (Hz)
DEFAULT_NWIN	= 200;			% Number of windows = horizontal resolution
DEFAULT_VARWIN	= 0;			% Fixed window length or base on cycles.
								% =0: fix window length to nwin
								% >0: set window length equal varwin cycles
								%     bounded above by winsize, also determines
								%     the min. freq. to be computed.
DEFAULT_OVERSMP	= 2;			% Number of times to oversample = vertical resolution
DEFAULT_MAXFREQ = 50;			% Maximum frequency to display (Hz)
DEFAULT_TITLE	= '';			% Figure title
DEFAULT_ALPHA   = 0.02;			% Default two-sided significance probability threshold
MARGIN          = 0.12;         % width of marginal plots
DEFAULT_VERTS   = [];           % default no vertical lines

if (nargin < 2)
	help crossf
	return
end

if (min(size(X))~=1 | length(X)<2)
	fprintf('crossf(): xdata must be a row or column vector.\n');
    return
elseif (min(size(Y))~=1 | length(Y)<2)
	fprintf('crossf(): ydata must be a row or column vector.\n');
    return
elseif (length(X) ~= length(Y))
	fprintf('crossf(): xdata and ydata must have same length.\n');
    return
end

if (nargin < 3)
	epoch = DEFAULT_EPOCH;
elseif (~isnumeric(epoch) | length(epoch)~=1 | epoch~=round(epoch))
	fprintf('crossf(): Value of frames must be an integer.\n');
    return
elseif (epoch <= 0)
	fprintf('crossf(): Value of frames must be positive.\n');
    return
elseif (rem(length(X),epoch) ~= 0)
	fprintf('crossf(): Length of data vectors must be divisible by frames.\n');
    return
end

if (nargin < 4)
	timelim = DEFAULT_TIMELIM;
elseif (~isnumeric(timelim) | sum(size(timelim))~=3)
	error('crossf(): Value of tlimits must be a vector containing two numbers.');
elseif (timelim(1) >= timelim(2))
	error('crossf(): tlimits interval must be [min,max].');
end

if (nargin < 5)
	ftitle = DEFAULT_TITLE;
elseif (~ischar(ftitle))
	error('crossf(): Plot title argument must be a quoted string.');
end

if (nargin < 6)
	Fs = DEFAULT_FS;
elseif (~isnumeric(Fs) | length(Fs)~=1)
	error('crossf(): Value of srate must be a number.');
elseif (Fs <= 0)
	error('crossf(): Value of srate must be positive.');
end

if (nargin < 7)
	varwin = DEFAULT_VARWIN;
elseif (~isnumeric(varwin) | length(varwin)~=1)
	error('crossf(): Value of cycles must be a number.');
elseif (varwin < MAX_BASELN)
	error('crossf(): Value of cycles must be either zero or positive.');
end

if (nargin < 8)
	winsize = max(pow2(nextpow2(epoch)-3),4);
elseif (~isnumeric(winsize) | length(winsize)~=1 | winsize~=round(winsize))
	error('crossf(): Value of winsize must be an integer number.');
elseif (winsize <= 0)
	error('crossf(): Value of winsize must be positive.');
elseif (varwin == 0 & pow2(nextpow2(winsize)) ~= winsize)
	error('crossf(): Value of winsize must be an integer power of two [1,2,4,8,16,...]');
elseif (winsize > epoch)
	error('crossf(): Value of winsize must be less than epoch length.');
end

if (nargin < 9)
	nwin = DEFAULT_NWIN;
elseif (~isnumeric(nwin) | length(nwin)~=1 | nwin~=round(nwin))
	error('crossf(): Value of nwin must be an integer number.');
elseif (nwin <= 0)
	error('crossf(): Value of nwin must be positive.');
end
if (nwin > epoch-winsize)
	error('crossf(): Value of nwin must be <= epoch-winsize.');
end

if (nargin < 10)
	oversmp = DEFAULT_OVERSMP;
elseif (~isnumeric(oversmp) | length(oversmp)~=1 | oversmp~=round(oversmp))
	error('crossf(): Value of oversmp must be an integer number.');
elseif (oversmp <= 0)
	error('crossf(): Value of oversmp must be positive.');
elseif (pow2(nextpow2(oversmp)) ~= oversmp)
	error('crossf(): Value of oversmp must be an integer power of two [1,2,4,8,16,...]');
end

if (nargin < 11)
	maxfreq = DEFAULT_MAXFREQ;
elseif (~isnumeric(maxfreq) | length(maxfreq)~=1)
	error('crossf(): Value of maxfreq must be a number.');
elseif (maxfreq <= 0)
	error('crossf(): Value of maxfreq must be positive.');
end

if (nargin < 12)
	alpha = DEFAULT_ALPHA;
elseif (~isnumeric(alpha) | length(alpha)~=1)
	error('crossf(): Value of alpha must be a number.');
elseif (round(NACCU*alpha) < 2 | alpha > .5)
	fprintf('crossf(): Value of alpha must be in the range (~0,0.5]');
    return 
else
    if round(NACCU*alpha)<1,
      alpha = 1/NACCU;
	  fprintf(...
        'Using alpha = %0.3f. To decrease, must raise NACCU in source code.\n',...
                       alpha);
    end
end
if (nargin < 13)
   verts = DEFAULT_VERTS;
end

if (varwin == 0) % FFT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	freqs = Fs/winsize*[1:2/oversmp:winsize]/2;
	win = hanning(winsize);

	R = zeros(oversmp*winsize/2,nwin);
	RR = zeros(oversmp*winsize/2,nwin);
	Rboot = zeros(oversmp*winsize/2,NACCU);

else % wavelet DFT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	freqs = Fs*varwin/winsize*[2:2/oversmp:winsize]/2;
	win = dftfilt(winsize,maxfreq/Fs,varwin,oversmp,1);

	R  = zeros(size(win,2),nwin);
	RR = zeros(size(win,2),nwin);
	Rboot = zeros(size(win,2),NACCU);
end

wintime = 500*winsize/Fs;
times = timelim(1)+wintime:(timelim(2)-timelim(1)-2*wintime)/(nwin-1):timelim(2)-wintime;

baseln = find(times < 0);
dispf = find(freqs <= maxfreq);
stp = (epoch-winsize)/(nwin-1);
trials = length(X)/epoch;