crossfold.m

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

fprintf('\nComputing the Event-Related Cross-Coherence image\n');
fprintf(' based on %d trials of %d frames sampled at %g Hz.\n',...
                       trials,epoch,Fs);
fprintf('Trial timebase is %d ms before to %d ms after the stimulus\n',...
                       timelim(1),timelim(2));
fprintf('The frequency range displayed is %g-%g Hz.\n',min(dispf),maxfreq);
if varwin==0
  fprintf('The data window size is %d samples (%g ms).\n',winsize,2*wintime);
  fprintf('The FFT length is %d samples\n',winsize*oversmp);
else
  fprintf('The window size is %d cycles.\n',varwin);
  fprintf('The maximum window size is %d samples (%g ms).\n',winsize,2*wintime);
end
fprintf('The window is applied %d times\n',nwin);
fprintf(' with an average step size of %g samples (%g ms).\n',...
                       stp,1000*stp/Fs);
fprintf('Results are oversampled %d times.\n',oversmp);
if nargout>4
  fprintf('Bootstrap confidence limits will be computed based on alpha = %g\n',...
              alpha);
else
  fprintf('Bootstrap confidence limits will NOT be computed.\n'); 
end
if nargout>5
  fprintf(['Coherence angles will be imaged in ',ANGLEUNITS,' and saved.\n']);
end

fprintf('\nProcessing trial (of %d):',trials);
firstboot = 1;
Rn=zeros(1,nwin);
X = X(:)'; % make X and Y column vectors
Y = Y(:)';
for t=1:trials,
	if (rem(t,10) == 0)
		fprintf(' %d',t);
	end
    if rem(t,120) == 0
        fprintf('\n');
    end

	for j=1:nwin, % for each time window
		tmpX = X([1:winsize]+floor((j-1)*stp)+(t-1)*epoch);
		tmpY = Y([1:winsize]+floor((j-1)*stp)+(t-1)*epoch);

        if ~any(isnan(tmpX))
		  tmpX = tmpX - mean(tmpX);
		  tmpY = tmpY - mean(tmpY);

		  if varwin == 0 % use FFTs
			tmpX = win .* tmpX(:);
			tmpY = win .* tmpY(:);
			tmpX = fft(tmpX,oversmp*winsize);
			tmpY = fft(tmpY,oversmp*winsize);
			tmpX = tmpX(2:oversmp*winsize/2+1);
			tmpY = tmpY(2:oversmp*winsize/2+1);
		  else 
			tmpX = win' * tmpX(:);
			tmpY = win' * tmpY(:);
		  end

          if nargout > 4
           if firstboot==1
             tmpsX = repmat(nan,length(tmpX),nwin);
             tmpsY = repmat(nan,length(tmpY),nwin);
             firstboot = 0;
           end
           tmpsX(:,j) = tmpX;
           tmpsY(:,j) = tmpY;
          end
		
		  RR(:,j) = tmpX.*conj(tmpY) ./ (abs(tmpX).*abs(tmpY)); % complex coher.
	      R(:,j) = R(:,j) + RR(:,j);
          Rn(j) = Rn(j)+1;
        end % ~any(isnan())
	end % time window
	
	if (nargout > 4) % get NACCU bootstrap estimates for each trial
        j=1;
		while j<=NACCU
           s = ceil(rand([1 2])*nwin); % random ints [1,nwin]
           tmpX = tmpsX(:,s(1));
           tmpY = tmpsY(:,s(2));
           if ~any(isnan(tmpX)) & ~any(isnan(tmpY))
		      RR = tmpX.*conj(tmpY) ./ (abs(tmpX).*abs(tmpY)); % complex coher.
		      Rboot(:,j) = Rboot(:,j) + RR;
              j = j+1;
           end
        end
	end
	
end % t = trial

fprintf('\nNow plotting...\n');

Rangle = angle(R);
if varwin ~= 0
   Rangle = -Rangle; % make lead/lag the same for FFT and wavelet analysis
end
R = abs(R) ./ (ones(size(R,1),1)*Rn);               % coherence magnitude
Rraw = R;						% raw coherence values
mbase = mean(R(:,baseln)');     % mean baseline coherence magnitude
% R = R - repmat(mbase',[1 nwin]);% remove baseline mean

if (nargout > 4) % bootstrap
	i = round(NACCU*alpha);
	Rboot = abs(Rboot) / trials; % normalize bootstrap magnitude to [0,1]
	Rboot = sort(Rboot');  
	Rsignif = mean(Rboot(NACCU-i+1:NACCU,:)); % significance levels for Rraw
%	Rboot = [mean(Rboot(1:i,:))-mbase ; mean(Rboot(NACCU-i+1:NACCU,:))-mbase];
 	Rboot = [mean(Rboot(1:i,:)) ; mean(Rboot(NACCU-i+1:NACCU,:))];
end % NOTE: above, mean ?????

set(gcf,'DefaultAxesFontSize',AXES_FONT)
colormap(jet(256));

pos = get(gca,'position'); % plot relative to current axes
q = [pos(1) pos(2) 0 0];
s = [pos(3) pos(4) pos(3) pos(4)];
axis('off')

if nargout>5 % image coherence lag as well as coherence magnitude
  ybase = 0.5;
  subheight = 0.4;
  MARGIN = MARGIN*0.75;
else
  ybase = 0.0;
  subheight = 0.9;
end

%
% Image the coherence [% perturbations]
%
RR = R;
if (nargout > 4) % zero out (and 'green out') nonsignif. R values
	RR(find((RR > repmat(Rboot(1,:)',[1 nwin])) ...
              & (RR < repmat(Rboot(2,:)',[1 nwin])))) = 0;
end
if (nargout > 5) % zero out nonsignif. Rraw values
	Rraw(find(repmat(Rsignif',[1,size(Rraw,2)])>=Rraw))=0;
end

if COH_CAXIS_LIMIT == 0
    coh_caxis = max(max(R(dispf,:)))*[-1 1];
else
    coh_caxis = COH_CAXIS_LIMIT*[-1 1];
end

h(6) = axes('Units','Normalized',...
               'Position',[MARGIN ybase+MARGIN 0.9-MARGIN subheight].*s+q);

map=hsv(300); % install circular color map - green=0, yellow, orng, red, violet = max
              %                                         cyan, blue, violet = min
map = flipud([map(251:end,:);map(1:250,:)]);
map(151,:) = map(151,:)*0.9; % tone down the (0=) green!
colormap(map);

imagesc(times,freqs(dispf),RR(dispf,:),coh_caxis); % plot the coherence image

set(h(6),'Units','Normalized',...
               'Position',[MARGIN ybase+MARGIN 0.9-MARGIN subheight].*s+q);
hold on
plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',LINEWIDTH)
for i=1:length(verts)
  plot([verts(i) verts(i)],[0 freqs(max(dispf))],'--m','LineWidth',LINEWIDTH);
end;
hold off
set(h(6),'YTickLabel',[],'YTick',[])
set(h(6),'XTickLabel',[],'XTick',[])
title('Event-Related Coherence')
h(8) = axes('Position',[.92 ybase+MARGIN .05 subheight].*s+q);
cbar(h(8),151:300,[0 coh_caxis(2)]); % use only positive colors (gyorv) 
%                                      for coherences
%
% Plot delta-mean min and max coherence at each time point on bottom of image
%
h(10) = axes('Units','Normalized','Position',[MARGIN ybase 0.9-MARGIN MARGIN].*s+q); 
Emax = max(R(dispf,:)); % mean coherence at each time point
Emin = min(R(dispf,:)); % mean coherence at each time point
plot(times,Emax,'b');
hold on
plot(times,Emin,'b');
plot([times(1) times(length(times))],[0 0],'LineWidth',0.7);
plot([0 0],[-500 500],'--m','LineWidth',LINEWIDTH);
for i=1:length(verts)
  plot([verts(i) verts(i)],[-500 500],'--m','LineWidth',LINEWIDTH);
end;
axis([min(times) max(times) 0 max(Emax)*1.2])
tick = get(h(10),'YTick');
set(h(10),'YTick',[tick(1) ; tick(length(tick))])
set(h(10),'YAxisLocation','right')
midpos = get(h(10),'Position');
if nargout<6
 xlabel('Time (ms)')
end
ylabel('coh.')

%
% Plot mean baseline coherence at each freq on left side of image
%

h(11) = axes('Units','Normalized','Position',[0 ybase+MARGIN MARGIN subheight].*s+q);
E = mbase(dispf); % baseline mean coherence at each frequency
if (nargout > 4) % plot bootstrap significance limits (base mean +/-)
	plot(freqs(dispf),E,'m','LineWidth',LINEWIDTH); % plot mbase
    hold on
	% plot(freqs(dispf),Rboot(:,dispf)+[E;E],'g','LineWidth',LINEWIDTH);
	plot(freqs(dispf),Rboot([1 2],dispf),'g','LineWidth',LINEWIDTH);
	plot(freqs(dispf),Rsignif(dispf),'k:','LineWidth',LINEWIDTH);
	axis([freqs(1) freqs(max(dispf)) 0 max([E Rsignif])*1.2]);
else             % plot marginal mean coherence only
	plot(freqs(dispf),E,'LineWidth',LINEWIDTH);
	% axis([freqs(1) freqs(max(dispf)) min(E)-max(E)/3 max(E)+max(E)/3]);
	if ~isnan(max(E))
	   axis([freqs(1) freqs(max(dispf)) 0 max(E)*1.2]);
	end;   
end

tick = get(h(11),'YTick');
set(h(11),'YTick',[tick(1) ; tick(length(tick))])
set(h(11),'View',[90 90])
xlabel('Freq. (Hz)')
ylabel('coh.')

if (length(ftitle) > 0) % plot title
	axes('Position',pos,'Visible','Off');               
	h(12) = text(-.05,1.01,ftitle);
	set(h(12),'VerticalAlignment','bottom')
	set(h(12),'HorizontalAlignment','left')
	set(h(12),'FontSize',TITLE_FONT)
end
%
% Plot coherence time lags in bottom panel
%
if nargout>5
   h(13) = axes('Units','Normalized','Position',[MARGIN MARGIN 0.9-MARGIN subheight].*s+q);
   if strcmp(ANGLEUNITS,'ms')  % convert to ms
     Rangle = (Rangle/(2*pi)).*repmat(1000./freqs(dispf)',1,length(times)); 
     maxangle = max(max(abs(Rangle)));
   else
     Rangle = Rangle*180/pi; % convert to degrees
     maxangle = 180; % use full-cycle plotting 
   end
   Rangle(find(Rraw==0)) = 0; % set angle at non-signif coher points to 0

   imagesc(times,freqs(dispf),Rangle(dispf,:),[-maxangle maxangle]); % plot the 
   hold on                                             % coherence phase angles
   plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',LINEWIDTH); % zero-time line
   for i=1:length(verts)
     plot([verts(i) verts(i)],[0 freqs(max(dispf))],'--m','LineWidth',LINEWIDTH);
   end;

   pos13 = get(h(13),'Position');
   set(h(13),'Position',[pos13(1) pos13(2) midpos(3) pos13(4)]);
   ylabel('Freq. (Hz)')
   xlabel('Time (ms)')
   h(14)=axes('Position',[.92 MARGIN .05 subheight].*s+q);
   cbar(h(14),0,[-maxangle maxangle]); % two-sided colorbar

   if (length(ftitle) > 0) % plot title
	axes('Position',pos,'Visible','Off');               
	h(13) = text(-.05,1.01,ftitle);
	set(h(13),'VerticalAlignment','bottom')
	set(h(13),'HorizontalAlignment','left')
	set(h(13),'FontSize',TITLE_FONT)
   end
 end
axcopy(gcf);