timef.m

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

								% >0: set window length equal to varwin cycles
								%     Bounded above by winsize, which determines
								%     the min. freq. to be computed.
DEFAULT_OVERSMP	= 2;			% Number of times to oversample frequencies 
DEFAULT_MAXFREQ = 50;			% Maximum frequency to display (Hz)
DEFAULT_TITLE	= '';			% Figure title
DEFAULT_ELOC    = 'chan.locs';	% Channel location file
DEFAULT_ALPHA   = NaN;			% Percentile of bins to keep
DEFAULT_MARKTIME= NaN;

% Font sizes:
AXES_FONT       = 10;           % axes text FontSize
TITLE_FONT      = 8;

if (nargin < 1)
	help timef
	return
end

if isstr(X) & strcmp(X,'details')
   more on
   help timefdetails
   more off
   return
end

if (min(size(X))~=1 | length(X)<2)
	error('Data must be a row or column vector.');
end

if (nargin < 2)
	frame = DEFAULT_EPOCH;
elseif (~isnumeric(frame) | length(frame)~=1 | frame~=round(frame))
	error('Value of frames must be an integer.');
elseif (frame <= 0)
	error('Value of frames must be positive.');
elseif (rem(length(X),frame) ~= 0)
	error('Length of data vector must be divisible by frames.');
end

if (nargin < 3)
	tlimits = DEFAULT_TIMLIM;
elseif (~isnumeric(tlimits) | sum(size(tlimits))~=3)
	error('Value of tlimits must be a vector containing two numbers.');
elseif (tlimits(1) >= tlimits(2))
	error('tlimits interval must be ascending.');
end

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

if (nargin < 5)
	varwin = DEFAULT_VARWIN;
elseif (~isnumeric(varwin) | length(varwin)>2)
	error('Value of cycles must be a number.');
elseif (varwin < 0)
	error('Value of cycles must be zero or positive.');
end

% consider structure for these arguments
% --------------------------------------
if ~isempty(varargin)
    try, g = struct(varargin{:}); 
    catch, error('Argument error in the {''param'', value} sequence'); end; 
end;
g.tlimits = tlimits;
g.frame   = frame;
g.srate   = Fs;
g.cycles  = varwin(1);
if length(varwin)>1
	g.cyclesfact = varwin(2);
else 
	g.cyclesfact = 1;
end;

try, g.title;      catch, g.title = DEFAULT_TITLE; end;
try, g.winsize;    catch, g.winsize = max(pow2(nextpow2(g.frame)-3),4); end;
try, g.pad;        catch, g.pad = max(pow2(nextpow2(g.winsize)),4); end;
try, g.timesout;   catch, g.timesout = DEFAULT_NWIN; end;
try, g.padratio;   catch, g.padratio = DEFAULT_OVERSMP; end;
try, g.maxfreq;    catch, g.maxfreq = DEFAULT_MAXFREQ; end;
try, g.topovec;    catch, g.topovec = []; end;
try, g.elocs;      catch, g.elocs = DEFAULT_ELOC; end;
try, g.alpha;      catch, g.alpha = DEFAULT_ALPHA; end;  
try, g.marktimes;  catch, g.marktimes = DEFAULT_MARKTIME; end;
try, g.powbase;    catch, g.powbase = NaN; end;
try, g.pboot;      catch, g.pboot = NaN; end;
try, g.rboot;      catch, g.rboot = NaN; end;
try, g.plotersp;   catch, g.plotersp = 'on'; end;
try, g.plotitc;    catch, g.plotitc  = 'on'; end;
try, g.detrep;     catch, g.detrep = 'off'; end;
try, g.detret;     catch, g.detret = 'off'; end;
try, g.baseline;   catch, g.baseline = 0; end;
try, g.baseboot;   catch, g.baseboot = 1; end;
try, g.linewidth;  catch, g.linewidth = 2; end;
try, g.naccu;      catch, g.naccu = 200; end;
try, g.mtaper;     catch, g.mtaper = []; end;
try, g.vert;       catch, g.vert = []; end;
try, g.type;       catch, g.type = 'phasecoher'; end;
try, g.phsamp;     catch, g.phsamp = 'off'; end;
try, g.plotphase;  catch, g.plotphase = 'on'; end;
try, g.itcmax;     catch, g.itcmax = []; end;
try, g.erspmax;    catch, g.erspmax = []; end;
try, g.verbose;    catch, g.verbose = 'on'; end;

% testing arguments consistency
% -----------------------------
switch lower(g.verbose)
    case { 'on', 'off' }, ;
    otherwise error('verbose must be either on or off');
end;
if (~ischar(g.title))
	error('Title must be a string.');
end

if (~isnumeric(g.winsize) | length(g.winsize)~=1 | g.winsize~=round(g.winsize))
	error('Value of winsize must be an integer number.');
elseif (g.winsize <= 0)
	error('Value of winsize must be positive.');
elseif (g.cycles == 0 & pow2(nextpow2(g.winsize)) ~= g.winsize)
	error('Value of winsize must be an integer power of two [1,2,4,8,16,...]');
elseif (g.winsize > g.frame)
	error('Value of winsize must be less than frame length.');
end

if (~isnumeric(g.timesout) | length(g.timesout)~=1 | g.timesout~=round(g.timesout))
	error('Value of timesout must be an integer number.');
elseif (g.timesout <= 0)
	error('Value of timesout must be positive.');
end
if (g.timesout > g.frame-g.winsize)
	g.timesout = g.frame-g.winsize;
	disp(['Value of timesout must be <= frame-winsize, timeout adjusted to ' int2str(g.timesout) ]);
end

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

if (~isnumeric(g.maxfreq) | length(g.maxfreq)~=1)
	error('Value of maxfreq must be a number.');
elseif (g.maxfreq <= 0)
	error('Value of maxfreq must be positive.');
elseif (g.maxfreq > Fs/2)
	myprintf(g.verbose,['Warning: value of maxfreq reduced to Nyquist rate' ...
		 ' (%3.2f)\n\n'], Fs/2);
	g.maxfreq = Fs/2;
end

if isempty(g.topovec)
	g.topovec = [];
	if isempty(g.elocs)
		error('Channel location file must be specified.');
	end;
end
if isempty(g.elocs)
	g.elocs = DEFAULT_ELOC;
elseif (~ischar(g.elocs)) & ~isstruct(g.elocs)
	error('Channel location file must be a valid text file.');
end

if (~isnumeric(g.alpha) | length(g.alpha)~=1)
	error('timef(): Value of g.alpha must be a number.\n');
elseif (round(g.naccu*g.alpha) < 2)
	myprintf(g.verbose,'Value of g.alpha is out of the normal range [%g,0.5]\n',2/g.naccu);
    g.naccu = round(2/g.alpha);
	myprintf(g.verbose,'  Increasing the number of bootstrap iterations to %d\n',g.naccu);
end
if g.alpha>0.5 | g.alpha<=0
    error('Value of g.alpha is out of the allowed range (0.00,0.5).');
end
if ~isnan(g.alpha)
   if g.baseboot > 0
     myprintf(g.verbose,'Bootstrap analysis will use data in baseline (pre-0) subwindows only.\n')
   else
     myprintf(g.verbose,'Bootstrap analysis will use data in all subwindows.\n')
   end
end
if ~isnumeric(g.vert)
    error('vertical line(s) option must be a vector');
else
	if min(g.vert) < g.tlimits(1) | max(g.vert) > g.tlimits(2)
        error('vertical line(s) time out-of-bound');
	end;
end;

if ~isnan (g.rboot)
  if size(g.rboot) == [1,1]
    if g.cycles == 0
        g.rboot = g.rboot*ones(g.winsize*g.padratio/2);
    end
  end
end;

if ~isempty(g.mtaper) % mutitaper, inspired from Bijan Pesaran matlab function
  if length(g.mtaper) < 3
        %error('mtaper arguement must be [N W] or [N W K]');
    
    if g.mtaper(1) * g.mtaper(2) < 1
        error('mtaper 2 first arguments'' product must be higher than 1');
    end;
    if length(g.mtaper) == 2
        g.mtaper(3) = floor( 2*g.mtaper(2)*g.mtaper(1) - 1);
    end
    if length(g.mtaper) == 3
        if g.mtaper(3) > 2 * g.mtaper(1) * g.mtaper(2) -1
            error('mtaper number too high (maximum (2*N*W-1))');
        end;
    end
    disp(['Using ' num2str(g.mtaper(3)) ' tapers.']);
    NW = g.mtaper(1)*g.mtaper(2);   % product NW
    N  = g.mtaper(1)*g.srate;     
    [e,v] = dpss(N, NW, 'calc');
    e=e(:,1:g.mtaper(3));
    g.alltapers = e;
  else    
    g.alltapers = g.mtaper;
    disp('mtaper argument not [N W] or [N W K]; considering raw taper matrix');
  end;
  g.winsize = size(g.alltapers, 1);
  g.pad = max(pow2(nextpow2(g.winsize)),256); % pad*nextpow
  %nfk = floor([0 g.maxfreq]./g.srate.*g.pad); % not used any more
  %g.padratio = 2*nfk(2)/g.winsize;
  g.padratio = g.pad/g.winsize;
 
  %compute number of frequencies
  %nf = max(256, g.pad*2^nextpow2(g.winsize+1)); 
  %nfk = floor([0 g.maxfreq]./g.srate.*nf);
  
  %freqs = linspace( 0, g.maxfreq, diff(nfk)); % this also work in the case of a FFT
  
end;           

switch lower(g.plotphase)
    case { 'on', 'off' }, ;
    otherwise error('plotphase must be either on or off');
end;
switch lower(g.plotersp)
    case { 'on', 'off' }, ;
    otherwise error('plotersp must be either on or off');
end;
switch lower(g.plotitc)
    case { 'on', 'off' }, ;
    otherwise error('plotitc must be either on or off');
end;
switch lower(g.detrep)
    case { 'on', 'off' }, ;
    otherwise error('detrep must be either on or off');
end;
switch lower(g.detret)
    case { 'on', 'off' }, ;
    otherwise error('detret must be either on or off');
end;
switch lower(g.phsamp)
    case { 'on', 'off' }, ;
    otherwise error('phsamp must be either on or off');
end;
if ~isnumeric(g.linewidth)
    error('linewidth must be numeric');
end;
if ~isnumeric(g.naccu)
    error('naccu must be numeric');
end;
if ~isnumeric(g.baseline)
    error('baseline must be numeric');
end;
switch g.baseboot
    case {0,1}, ;
    otherwise, error('baseboot must be 0 or 1');
end;
switch g.type
    case { 'coher', 'phasecoher', 'phasecoher2' },;
    otherwise error('Type must be either ''coher'' or ''phasecoher''');
end;    

if (g.cycles == 0) %%%%%%%%%%%%%% constant window-length FFTs %%%%%%%%%%%%%%%%
    %freqs = g.srate/g.winsize*[1:2/g.padratio:g.winsize]/2 % incorect for padratio > 2
    freqs = linspace(0, g.srate/2, g.padratio*g.winsize/2+1);
    freqs = freqs(2:end);
    win = hanning(g.winsize);

    P  = zeros(g.padratio*g.winsize/2,g.timesout); % summed power
	PP = zeros(g.padratio*g.winsize/2,g.timesout); % power
	R  = zeros(g.padratio*g.winsize/2,g.timesout); % mean coherence
	RR = zeros(g.padratio*g.winsize/2,g.timesout); % (coherence)
	Pboot = zeros(g.padratio*g.winsize/2,g.naccu); % summed bootstrap power
	Rboot = zeros(g.padratio*g.winsize/2,g.naccu); % summed bootstrap coher
    Rn = zeros(1,g.timesout);
    Rbn = 0;
	switch g.type
	    case { 'coher' 'phasecoher2' },
           cumulX = zeros(g.padratio*g.winsize/2,g.timesout);
           cumulXboot = zeros(g.padratio*g.winsize/2,g.naccu);
        case 'phasecoher'
           switch g.phsamp
             case 'on'
               cumulX = zeros(g.padratio*g.winsize/2,g.timesout);
           end
    end;        

else % %%%%%%%%%%%%%%%%%% cycles>0, Constant-Q (wavelet) DFTs %%%%%%%%%%%%%%%%%%%%

   	freqs = g.srate*g.cycles/g.winsize*[2:2/g.padratio:g.winsize]/2;
    dispf = find(freqs <= g.maxfreq);
    freqs = freqs(dispf);

	%win = dftfilt(g.winsize,g.maxfreq/g.srate,g.cycles,g.padratio,0.5);
	win = dftfilt(g.winsize,g.maxfreq/g.srate,g.cycles,g.padratio,g.cyclesfact);
	P = zeros(size(win,2),g.timesout);       % summed power
	R = zeros(size(win,2),g.timesout);       % mean coherence
	PP = repmat(NaN,size(win,2),g.timesout); % initialize with NaN
	RR = repmat(NaN,size(win,2),g.timesout); % initialize with NaN
	Pboot = zeros(size(win,2),g.naccu);  % summed bootstrap power
	Rboot = zeros(size(win,2),g.naccu);  % summed bootstrap coher
    Rn = zeros(1,g.timesout);
    Rbn = 0;
	switch g.type
	  case { 'coher' 'phasecoher2' },
           cumulX = zeros(size(win,2),g.timesout);
           cumulXboot = zeros(size(win,2),g.naccu);
      case 'phasecoher'
           switch g.phsamp
             case 'on'
               cumulX = zeros(size(win,2),g.timesout);
           end
   end;        
end

switch g.phsamp
  case 'on'
    PA = zeros(size(P,1),size(P,1),g.timesout); % NB: (freqs,freqs,times)
end                                             %       phs   amp

wintime = 1000/g.srate*(g.winsize/2); % (1000/g.srate)*(g.winsize/2);
times = [g.tlimits(1)+wintime:(g.tlimits(2)-g.tlimits(1)-2*wintime)/(g.timesout-1):g.tlimits(2)-wintime];
ERPtimes = [g.tlimits(1):(g.tlimits(2)-g.tlimits(1))/(g.frame-1):g.tlimits(2)+0.000001];
ERPindices = [];
for ti=times
 [tmp indx] = min(abs(ERPtimes-ti));
 ERPindices  = [ERPindices indx];
end
ERPtimes = ERPtimes(ERPindices); % subset of ERP frames on t/f window centers

if ~isempty(find(times < g.baseline))
   baseln = find(times < g.baseline); % subtract means of pre-0 (centered) windows
else
   baseln = 1:length(times); % use all times as baseline
end
if ~isnan(g.alpha) & length(baseln)==0
  myprintf(g.verbose,'timef(): no window centers in baseline (times<%g) - shorten (max) window length.\n', g.baseline)
  return
elseif ~isnan(g.alpha) & g.baseboot
  myprintf(g.verbose,'   %d bootstrap windows in baseline (times<%g).\n',...
          length(baseln), g.baseline)
end
dispf = find(freqs <= g.maxfreq);
stp = (g.frame-g.winsize)/(g.timesout-1);

myprintf(g.verbose,'Computing Event-Related Spectral Perturbation (ERSP) and\n');
switch g.type
    case 'phasecoher',  myprintf(g.verbose,'  Inter-Trial Phase Coherence (ITC) images based on %d trials\n',length(X)/g.frame);
    case 'phasecoher2', myprintf(g.verbose,'  Inter-Trial Phase Coherence 2 (ITC) images based on %d trials\n',length(X)/g.frame);
    case 'coher',       myprintf(g.verbose,'  Linear Inter-Trial Coherence (ITC) images based on %d trials\n',length(X)/g.frame);
end;
myprintf(g.verbose,'  of %d frames sampled at %g Hz.\n',g.frame,g.srate);
myprintf(g.verbose,'Each trial contains samples from %d ms before to\n',g.tlimits(1));
myprintf(g.verbose,'  %d ms after the timelocking event.\n',g.tlimits(2));
myprintf(g.verbose,'The window size used is %d samples (%g ms) wide.\n',g.winsize,2*wintime);
myprintf(g.verbose,'The window is applied %d times at an average step\n',g.timesout);
myprintf(g.verbose,'  size of %g samples (%gms).\n',stp,1000*stp/g.srate);
myprintf(g.verbose,'Results are oversampled %d times; the %d frequencies\n',g.padratio,length(dispf));
myprintf(g.verbose,'  displayed are from %2.1f Hz to %3.1f Hz.\n',freqs(dispf(1)),freqs(dispf(end)));
if ~isnan(g.alpha)
  myprintf(g.verbose,'Only significant values (bootstrap p<%g) will be colored;\n',g.alpha) 
  myprintf(g.verbose,'  non-significant values will be plotted in green\n');
end

trials = length(X)/g.frame;
baselength = length(baseln);
myprintf(g.verbose,'\nOf %d trials total, processing trial:',trials);

% detrend over epochs (trials) if requested
% -----------------------------------------
switch g.detrep
    case 'on'
        X = reshape(X, g.frame, length(X)/g.frame);
        X = X - mean(X,2)*ones(1, length(X(:))/g.frame);
        X = X(:)';
end;        

for i=1:trials
    if (rem(i,100)==0)
        myprintf(g.verbose,'\n');
    end
	if (rem(i,10) == 0)
		myprintf(g.verbose,'%d',i);
	elseif (rem(i,2) == 0)
		myprintf(g.verbose,'.');
	end

    ERP = blockave(X,g.frame); % compute the ERP trial average

    Wn = zeros(1,g.timesout);
	for j=1:g.timesout,
		tmpX = X([1:g.winsize]+floor((j-1)*stp)+(i-1)*g.frame); 
                                                      % pull out data g.frame
		tmpX = tmpX - mean(tmpX); % remove the mean for that window
        switch g.detret, case 'on', tmpX = detrend(tmpX); end;