timef.m

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

		if ~any(isnan(tmpX))
		  if (g.cycles == 0) % FFT
            if ~isempty(g.mtaper)   % apply multitaper (no hanning window)
                tmpXMT = fft(g.alltapers .* ...
                             (tmpX(:) * ones(1,size(g.alltapers,2))), g.pad);
			    %tmpXMT = tmpXMT(nfk(1)+1:nfk(2),:);
			    tmpXMT = tmpXMT(2:g.padratio*g.winsize/2+1,:);
                PP(:,j) = mean(abs(tmpXMT).^2, 2); 
                  % power; can also ponderate multitaper by their eigenvalues v
		        tmpX = win .* tmpX(:);
               	tmpX = fft(tmpX, g.pad);
    		    tmpX = tmpX(2:g.padratio*g.winsize/2+1);
            else
		        tmpX = win .* tmpX(:);
               	tmpX = fft(tmpX,g.padratio*g.winsize);
    		    tmpX = tmpX(2:g.padratio*g.winsize/2+1);
                PP(:,j) = abs(tmpX).^2; % power
            end;    
          else % wavelet
            if ~isempty(g.mtaper)  % apply multitaper
			    tmpXMT = g.alltapers .* (tmpX(:) * ones(1,size(g.alltapers,2)));
			    tmpXMT = transpose(win) * tmpXMT;
                PP(:,j) = mean(abs(tmpXMT).^2, 2); % power
		        tmpX = transpose(win) * tmpX(:);
            else
		        tmpX = transpose(win) * tmpX(:); 
                PP(:,j) = abs(tmpX).^2; % power
            end    
          end
		
          if abs(tmpX) < MIN_ABS
		        RR(:,j) = zeros(size(RR(:,j)));
          else
		      switch g.type
		        case { 'coher' },
		          RR(:,j) = tmpX; 
                  cumulX(:,j) = cumulX(:,j)+abs(tmpX).^2;
		        case { 'phasecoher2' },
		          RR(:,j) = tmpX; 
                  cumulX(:,j) = cumulX(:,j)+abs(tmpX);
		        case 'phasecoher',
		          RR(:,j) = tmpX ./ abs(tmpX); % normalized cross-spectral vector
         		  switch g.phsamp
             		  case 'on'
                	    cumulX(:,j) = cumulX(:,j)+abs(tmpX); % accumulate for PA
          		  end
              end;
         end
          Wn(j) = 1;
        end

        switch g.phsamp
         case 'on' %PA (freq x freq x time)
          PA(:,:,j) = PA(:,:,j)  + (tmpX ./ abs(tmpX)) * ((PP(:,j)))';
                                           % x-product: unit phase column
                                           % times amplitude row
        end
	end % window

	%figure; imagesc(times, freqs, angle(RR));
	%figure; imagesc(times, freqs, log(PP));
	
	if ~isnan(g.alpha) % save surrogate data for bootstrap analysis
        j = 1;
        goodbasewins = find(Wn==1);
        if g.baseboot % use baseline windows only
          goodbasewins = find(goodbasewins<=baselength); 
        end
        ngdbasewins = length(goodbasewins);
        if ngdbasewins>1
		  while j <= g.naccu
            i=ceil(rand*ngdbasewins);
            i=goodbasewins(i);
			Pboot(:,j) = Pboot(:,j) + PP(:,i);
		    Rboot(:,j) = Rboot(:,j) + RR(:,i);
		    switch g.type
		        case 'coher',       cumulXboot(:,j) = cumulXboot(:,j)+abs(tmpX).^2;
		        case 'phasecoher2', cumulXboot(:,j) = cumulXboot(:,j)+abs(tmpX);
            end;
            j = j+1;
          end
          Rbn = Rbn + 1;
	    end
	end % bootstrap
	
    Wn = find(Wn>0);
    if length(Wn)>0
	  P(:,Wn) = P(:,Wn) + PP(:,Wn); % add non-NaN windows
	  R(:,Wn) = R(:,Wn) + RR(:,Wn);
	  Rn(Wn) = Rn(Wn) + ones(1,length(Wn)); % count number of addends
    end
end % trial

% if coherence, perform the division
% ----------------------------------
switch g.type
 case 'coher',
  R = R ./ ( sqrt( trials*cumulX ) );
  if ~isnan(g.alpha)
	  Rboot = Rboot ./ ( sqrt( trials*cumulXboot ) );
  end;
 case 'phasecoher2',
  R = R ./ ( cumulX );
  if ~isnan(g.alpha)
	  Rboot = Rboot ./ cumulXboot;
  end;   
 case 'phasecoher',
  R = R ./ (ones(size(R,1),1)*Rn);
end;        

switch g.phsamp
 case 'on'
  tmpcx(1,:,:) = cumulX; % allow ./ below
  for j=1:g.timesout
    PA(:,:,j) = PA(:,:,j) ./ repmat(PP(:,j)', [size(PP,1) 1]);
  end
end

if min(Rn) < 1
  myprintf(g.verbose,'timef(): No valid timef estimates for windows %s of %d.\n',...
                         int2str(find(Rn==0)),length(Rn));
  Rn(find(Rn<1))==1;
  return
end
P = P ./ (ones(size(P,1),1) * Rn);

if isnan(g.powbase)
  myprintf(g.verbose,'\nComputing the mean baseline spectrum\n');
  mbase = mean(P(:,baseln),2)';
else
  myprintf(g.verbose,'Using the input baseline spectrum\n');
  mbase = g.powbase;
end
if ~isnan( g.baseline ) & ~isnan( mbase )
    P = 10 * (log10(P) - repmat(log10(mbase(1:size(P,1)))',[1 g.timesout])); % convert to (10log10) dB
else
    P = 10 * log10(P);
end;

Rsign = sign(imag(R));
R = abs(R); % convert coherence vector to magnitude

if ~isnan(g.alpha) % if bootstrap analysis included . . .
    if Rbn>0
	   i = round(g.naccu*g.alpha);
       if isnan(g.pboot)
            Pboot = Pboot / Rbn; % normalize
            if ~isnan( g.baseline ) 
	            Pboot = 10 * (log10(Pboot) - repmat(log10(mbase)',[1 g.naccu]));
            else
                Pboot = 10 * log10(Pboot);
            end;  
            Pboot = sort(Pboot');
	        Pboot = [mean(Pboot(1:i,:)) ; mean(Pboot(g.naccu-i+1:g.naccu,:))];
       else
            Pboot = g.pboot;
       end
  
      if isnan(g.rboot)
	       Rboot = abs(Rboot) / Rbn;
	       Rboot = sort(Rboot');
	       Rboot = mean(Rboot(g.naccu-i+1:g.naccu,:));
      else
           Rboot = g.rboot;
      end
    else
      myprintf(g.verbose,'No valid bootstrap trials...!\n');
    end
end

switch lower(g.plotitc)
   case 'on',  
       switch lower(g.plotersp), 
          case 'on', ordinate1 = 0.67; ordinate2 = 0.1; height = 0.33; g.plot = 1;
          case 'off', ordinate2 = 0.1; height = 0.9; g.plot = 1;
       end;     
   case 'off', ordinate1 = 0.1; height = 0.9; 
       switch lower(g.plotersp), 
          case 'on', ordinate1 = 0.1; height = 0.9;  g.plot = 1;
          case 'off', g.plot = 0;
       end;     
end;    

if g.plot
    myprintf(g.verbose,'\nNow plotting...\n');
    set(gcf,'DefaultAxesFontSize',AXES_FONT)
    colormap(jet(256));
    pos = get(gca,'position');
    q = [pos(1) pos(2) 0 0];
    s = [pos(3) pos(4) pos(3) pos(4)];
end;

switch lower(g.plotersp)
 case 'on' 
    %
    %%%%%%% image the ERSP %%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    
	h(1) = subplot('Position',[.1 ordinate1 .9 height].*s+q);
	
	PP = P;
	if ~isnan(g.alpha) % zero out nonsignif. power differences
		PP(find((PP > repmat(Pboot(1,:)',[1 g.timesout])) ...
                    & (PP < repmat(Pboot(2,:)',[1 g.timesout])))) = 0;
	end

    if ERSP_CAXIS_LIMIT == 0
	   ersp_caxis = [-1 1]*1.1*max(max(abs(P(dispf,:))));
    else
       ersp_caxis = ERSP_CAXIS_LIMIT*[-1 1];
    end

    if ~isnan( g.baseline ) 
        imagesc(times,freqs(dispf),PP(dispf,:),ersp_caxis); 
    else
        imagesc(times,freqs(dispf),PP(dispf,:));
    end;
	if ~isempty(g.erspmax)
		caxis([-g.erspmax g.erspmax]);
	end;
    
	hold on
	plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth); % plot time 0
    if ~isnan(g.marktimes) % plot marked time
     for mt = g.marktimes(:)'
	   plot([mt mt],[0 freqs(max(dispf))],'--k','LineWidth',g.linewidth);
     end
    end
	hold off
	set(h(1),'YTickLabel',[],'YTick',[])
	set(h(1),'XTickLabel',[],'XTick',[])
	if ~isempty(g.vert)
		for index = 1:length(g.vert)
			line([g.vert(index), g.vert(index)], [min(freqs(dispf)) max(freqs(dispf))], 'linewidth', 1, 'color', 'm');
		end;
	end;

	h(2) = gca;
	h(3) = cbar('vert'); % ERSP colorbar axes
	set(h(2),'Position',[.1 ordinate1 .8 height].*s+q)
	set(h(3),'Position',[.95 ordinate1 .05 height].*s+q)
	title('ERSP (dB)')

	E = [min(P(dispf,:));max(P(dispf,:))];
	h(4) = subplot('Position',[.1 ordinate1-0.1 .8 .1].*s+q); % plot marginal ERSP means
                                                    % below the ERSP image
	plot(times,E,[0 0],...
	     [min(E(1,:))-max(max(abs(E)))/3 max(E(2,:))+max(max(abs(E)))/3], ...
             '--m','LineWidth',g.linewidth)
	axis([min(times) max(times) ...
             min(E(1,:))-max(max(abs(E)))/3 max(E(2,:))+max(max(abs(E)))/3])
	
	tick = get(h(4),'YTick');
	set(h(4),'YTick',[tick(1) ; tick(end)])
	set(h(4),'YAxisLocation','right')
    set(h(4),'TickLength',[0.020 0.025]);
	xlabel('Time (ms)')
	ylabel('dB')

	E = 10 * log10(mbase(dispf));
	h(5) = subplot('Position',[0 ordinate1 .1 height].*s+q); % plot mean spectrum
                                                    % to left of ERSP image
    plot(freqs(dispf),E,'LineWidth',g.linewidth)
	if ~isnan(g.alpha)
        hold on;
		plot(freqs(dispf),Pboot(:,dispf)+[E;E],'g', 'LineWidth',g.linewidth);
		plot(freqs(dispf),Pboot(:,dispf)+[E;E],'k:','LineWidth',g.linewidth)
	end

	axis([freqs(1) freqs(max(dispf)) min(E)-max(abs(E))/3 max(E)+max(abs(E))/3])
	tick = get(h(5),'YTick');
    if (length(tick)>1)
	   set(h(5),'YTick',[tick(1) ; tick(end-1)])
    end
    set(h(5),'TickLength',[0.020 0.025]);
	set(h(5),'View',[90 90])
	xlabel('Frequency (Hz)')
	ylabel('dB')
end;

switch lower(g.plotitc)
  case 'on'
    %
    %%%%%%%%%%%% Image the ITC %%%%%%%%%%%%%%%%%%
    %
	h(6) = subplot('Position',[.1 ordinate2 .9 height].*s+q); % ITC image

	RR = R;
	if ~isnan(g.alpha)
		RR(find(RR < repmat(Rboot(1,:)',[1 g.timesout]))) = 0;
	end

    if ITC_CAXIS_LIMIT == 0
	   coh_caxis = min(max(max(R(dispf,:))),1)*[-1 1]; % 1 WAS 0.4 !
    else
       coh_caxis = ITC_CAXIS_LIMIT*[-1 1];
    end

	if exist('Rsign') & strcmp(g.plotphase, 'on')
		imagesc(times,freqs(dispf),Rsign(dispf,:).*RR(dispf,:),coh_caxis); % <---
	else
		imagesc(times,freqs(dispf),RR(dispf,:),coh_caxis); % <---
	end
	if ~isempty(g.itcmax)
		caxis([-g.itcmax g.itcmax]);
	end;
	tmpcaxis = caxis;

	hold on
	plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth);
    if ~isnan(g.marktimes)
     for mt = g.marktimes(:)'
	   plot([mt mt],[0 freqs(max(dispf))],'--k','LineWidth',g.linewidth);
     end
    end
	hold off
	set(h(6),'YTickLabel',[],'YTick',[])
	set(h(6),'XTickLabel',[],'XTick',[])
	if ~isempty(g.vert)
		for index = 1:length(g.vert)
			line([g.vert(index), g.vert(index)], [min(freqs(dispf)) max(freqs(dispf))], 'linewidth', 1, 'color', 'm');
		end;
	end;

	h(7) = gca;
	h(8) = cbar('vert');
	%h(9) = get(h(8),'Children');
	set(h(7),'Position',[.1 ordinate2 .8 height].*s+q)
	set(h(8),'Position',[.95 ordinate2 .05 height].*s+q)
	set(h(8),'YLim',[0 tmpcaxis(2)]); 
	title('ITC')

    %
    %%%%% plot the ERP below the ITC image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
	% E = mean(R(dispf,:));

    ERPmax = max(ERP);
    ERPmin = min(ERP);
    ERPmax = ERPmax + 0.1*(ERPmax-ERPmin);
    ERPmin = ERPmin - 0.1*(ERPmax-ERPmin);
	h(10) = subplot('Position',[.1 ordinate2-0.1 .8 .1].*s+q); % ERP

    plot(ERPtimes,ERP(ERPindices),...
         [0 0],[ERPmin ERPmax],'--m','LineWidth',g.linewidth);
    hold on; plot([times(1) times(length(times))],[0 0], 'k');
    axis([min(ERPtimes) max(ERPtimes) ERPmin ERPmax]);

	tick = get(h(10),'YTick');
	set(h(10),'YTick',[tick(1) ; tick(end)])
    set(h(10),'TickLength',[0.02 0.025]);
	set(h(10),'YAxisLocation','right')
	xlabel('Time (ms)')
	ylabel('uV')

	E = mean(R(dispf,:)');
	h(11) = subplot('Position',[0 ordinate2 .1 height].*s+q); % plot the marginal mean
                                                    % ITC left of the ITC image
	if ~isnan(g.alpha)
		plot(freqs(dispf),E,'LineWidth',g.linewidth); hold on;
        plot(freqs(dispf),Rboot(dispf),'g', 'LineWidth',g.linewidth);
        plot(freqs(dispf),Rboot(dispf),'k:','LineWidth',g.linewidth);
        axis([freqs(1) freqs(max(dispf)) 0 max([E Rboot(dispf)])+max(E)/3])
	else
		plot(freqs(dispf),E,'LineWidth',g.linewidth)
		axis([freqs(1) freqs(max(dispf)) min(E)-max(E)/3 max(E)+max(E)/3])
	end

	tick = get(h(11),'YTick');
	set(h(11),'YTick',[tick(1) ; tick(length(tick))])
	set(h(11),'View',[90 90])
    set(h(11),'TickLength',[0.020 0.025]);
	xlabel('Frequency (Hz)')
	ylabel('ERP')
    %
    %%%%%%%%%%%%%%% plot a topoplot() %%%%%%%%%%%%%%%%%%%%%%%
    %
	if (~isempty(g.topovec))
		h(12) = subplot('Position',[-.1 .43 .2 .14].*s+q);
		if length(g.topovec) == 1
		  topoplot(g.topovec,g.elocs,'electrodes','off', ...
			   'style', 'blank', 'emarkersize1chan', 10);
		else
		  topoplot(g.topovec,g.elocs,'electrodes','off');
		end;
		    axis('square')
	end
end; %switch

if g.plot
	try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end;
    if (length(g.title) > 0)
	    axes('Position',pos,'Visible','Off');               
	    h(13) = text(-.05,1.01,g.title);
	    set(h(13),'VerticalAlignment','bottom')     
	    set(h(13),'HorizontalAlignment','left') 
        set(h(13),'FontSize',TITLE_FONT);
    end

    axcopy(gcf);
end;

% syemtric hanning function
function w = hanning(n)
if ~rem(n,2)
   w = .5*(1 - cos(2*pi*(1:n/2)'/(n+1)));
   w = [w; w(end:-1:1)];
else
   w = .5*(1 - cos(2*pi*(1:(n+1)/2)'/(n+1)));
   w = [w; w(end-1:-1:1)];
end

function myprintf(verbose, varargin)
    if strcmpi(verbose, 'on')
        fprintf(varargin{:});
    end;