crossf.m

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

   for index = 1:g.shuffle
      XX = shuffle(XX,3);
      X = [X XX(:,:)];
      Y = [Y YY];
   end;
end;

% 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);
   Y = reshape(Y, g.frame, length(Y)/g.frame);
   Y = Y - mean(Y,2)*ones(1, length(Y(:))/g.frame);
end;        

% time limits
wintime = 500*g.winsize/g.srate;
times = [g.tlimits(1)+wintime:(g.tlimits(2)-g.tlimits(1)-2*wintime)/(g.timesout-1):g.tlimits(2)-wintime];

%%%%%%%%%%
% baseline
%%%%%%%%%%
if ~isnan(g.baseline)
   baseln = find(times < g.baseline); % subtract means of pre-0 (centered) windows
   if isempty(baseln)
      baseln = 1:length(times); % use all times as baseline
      disp('Bootstrap baseline empty, using the whole epoch');
   end;
   baselength = length(baseln);
else
   baseln = 1:length(times); % use all times as baseline
   baselength = length(times); % used for bootstrap
end;

%%%%%%%%%%%%%%%%%%%%
% Initialize objects
%%%%%%%%%%%%%%%%%%%%
tmpsaveall = (~isnan(g.alpha) & isnan(g.rboot) & strcmp(g.memory, 'high')) ...
                 | (strcmp(g.subitc, 'on') & strcmp(g.memory, 'high'));
trials = length(X)/g.frame;
if ~strcmp(lower(g.compute), 'c')
   Tfx = tfinit(X, g.timesout, g.winsize, g.cycles, g.frame, g.padratio, g.detret, g.srate, g.maxfreq, g.subitc, g.type, g.cyclesfact, tmpsaveall);
   Tfy = tfinit(Y, g.timesout, g.winsize, g.cycles, g.frame, g.padratio, g.detret, g.srate, g.maxfreq, g.subitc, g.type, g.cyclesfact, tmpsaveall);
   Coher     = coherinit(Tfx.nb_points, trials, g.timesout, g.type);
   Coherboot = coherinit(Tfx.nb_points, trials, g.naccu   , g.type);
   Boot      = bootinit( Coherboot, Tfx.nb_points, g.timesout, g.naccu, baselength, g.baseboot, g.boottype, g.alpha, g.rboot);
   freqs = Tfx.freqs;
   dispf = find(freqs <= g.maxfreq);
   freqs = freqs(dispf);
else
   freqs = g.srate*g.cycles/g.winsize*[2:2/g.padratio:g.winsize]/2;
end;
dispf     = find(Tfx.freqs <= g.maxfreq);

%%%%%%%%%%%%%%
% reserve space
%%%%%%%%%%%%%%
%R  = zeros(tfx.nb_points,g.timesout);       % mean coherence
%RR = repmat(nan,tfx.nb_points,g.timesout); % initialize with nans
%Rboot = zeros(tfx.nb_points,g.naccu);  % summed bootstrap coher
%switch g.type
%ase 'coher',
%   cumulXY = zeros(tfx.nb_points,g.timesout);
%   cumulXYboot = zeros(tfx.nb_points,g.naccu);
%end;        
%if g.bootsub > 0
%   Rboottrial = zeros(tfx.nb_points, g.timesout, g.bootsub); % summed bootstrap coher
%   cumulXYboottrial = zeros(tfx.nb_points, g.timesout, g.bootsub);
%nd;
%if ~isnan(g.alpha) & isnan(g.rboot)
%   tf.tmpalltimes = repmat(nan,tfx.nb_points,g.timesout);
%nd
   
% display text to user
% --------------------
fprintf('\nComputing the Event-Related \n');
switch g.type
    case 'phasecoher',  fprintf('Phase Coherence (ITC) images based on %d trials\n',length(X)/g.frame);
    case 'phasecoher2', fprintf('Phase Coherence 2 (ITC) images based on %d trials\n',length(X)/g.frame);
    case 'coher',       fprintf('Linear Coherence (ITC) images based on %d trials\n',length(X)/g.frame);
end;
fprintf('Trial timebase is %d ms before to %d ms after the stimulus\n', g.tlimits(1),g.tlimits(2));
if ~isnan(g.baseline)
   if length(baseln) == length(times)
      fprintf('Using full time range as baseline\n');
   else
      fprintf('Using times in under %d ms for baseline\n', g.baseline);
   end;
else 
   fprintf('No baseline time range specified.\n');	
end;
fprintf('The frequency range displayed is %g-%g Hz.\n',min(dispf),g.maxfreq);
if g.cycles==0
   fprintf('The data window size is %d samples (%g ms).\n',g.winsize,2*wintime);
   fprintf('The FFT length is %d samples\n',g.winsize*g.padratio);
else
   fprintf('The window size is %d cycles.\n',g.cycles);
   fprintf('The maximum window size is %d samples (%g ms).\n',g.winsize,2*wintime);
end
fprintf('The window is applied %d times\n',g.timesout);
fprintf(' with an average step size of %g samples (%g ms).\n', Tfx.stp,1000*Tfx.stp/g.srate);
fprintf('Results are oversampled %d times.\n',g.padratio);
if ~isnan(g.alpha)
   fprintf('Bootstrap confidence limits will be computed based on alpha = %g\n', g.alpha);
else
   fprintf('Bootstrap confidence limits will NOT be computed.\n'); 
end
switch g.plotphase
case 'on', fprintf(['Coherence angles will be imaged in ',g.angleunit,'\n']);
end;
fprintf('\nProcessing trial (of %d):\n',trials);

%firstboot = 1;
%Rn=zeros(trials,g.timesout);
%X = X(:)'; % make X and Y column vectors
%Y = Y(:)';
%tfy = tfx;

if strcmp(lower(g.compute), 'c')
   % C PART
   filename = [ 'tmpcrossf' num2str(round(rand(1)*1000)) ];
   f = fopen([ filename '.in'], 'w');
   fwrite(f, tmpsaveall, 'int32');
   fwrite(f, g.detret, 'int32');
   fwrite(f, g.srate, 'int32');
   fwrite(f, g.maxfreq, 'int32');
   fwrite(f, g.padratio, 'int32');
   fwrite(f, g.cycles, 'int32');
   fwrite(f, g.winsize, 'int32');
   fwrite(f, g.timesout, 'int32');
   fwrite(f, g.subitc, 'int32');
   fwrite(f, g.type, 'int32');
   fwrite(f, trials, 'int32');
   fwrite(f, g.naccu, 'int32');
   fwrite(f, length(X), 'int32');
   fwrite(f, X, 'double');
   fwrite(f, Y, 'double');
   fclose(f);
   
   command = [ '!cppcrosff ' filename '.in ' filename '.out' ];
   eval(command);
   
   f = fopen([ filename '.out'], 'r');
   size1 = fread(f, 'int32', 1);
   size2 = fread(f, 'int32', 1);
   Rreal = fread(f, 'double', [size1 size2]);
   Rimg  = fread(f, 'double', [size1 size2]);
   Coher.R = Rreal + j*Rimg;
   Boot.Coherboot.R = [];
   Boot.Rsignif = [];
else
   % MATLAB PART
   % compute ITC if necessary
   % ------------------------
   if strcmp(g.subitc, 'on')
      for t=1:trials
         if rem(t,10) == 0,  fprintf(' %d',t); end
         if rem(t,120) == 0, fprintf('\n'); end
         Tfx = tfitc( Tfx, t, 1:g.timesout); 
         Tfy = tfitc( Tfy, t, 1:g.timesout); 
      end; 
	  fprintf('\n');
      Tfx = tfitcpost( Tfx, trials); 
      Tfy = tfitcpost( Tfy, trials); 
   end;
   
   % main loop
   % ---------
   if g.savecoher,
	   trialcoher = zeros(Tfx.nb_points, g.timesout, trials);   
   else  
	   trialcoher = [];
   end;
   for t=1:trials
      if rem(t,10) == 0,  fprintf(' %d',t); end
      if rem(t,120) == 0, fprintf('\n'); end
      
      Tfx = tfcomp( Tfx, t, 1:g.timesout); 
      Tfy = tfcomp( Tfy, t, 1:g.timesout);
	  if g.savecoher
		  [Coher trialcoher(:,:,t)] = cohercomp( Coher, Tfx.tmpalltimes, Tfy.tmpalltimes, t, 1:g.timesout);      
      else
		  Coher = cohercomp( Coher, Tfx.tmpalltimes, Tfy.tmpalltimes, t, 1:g.timesout);      
	  end;
	  
      Boot = bootcomp( Boot, Coher.Rn(t,:), Tfx.tmpalltimes, Tfy.tmpalltimes);
   end % t = trial
   [Boot Rbootout] = bootcomppost(Boot, Coher.Rn, Tfx.tmpall, Tfy.tmpall);
   % not that the bootstrap thresholding is actually performed in the display subfunction plotall()
   Coher  = cohercomppost(Coher, trials);
end;

% if coherence, perform the division
% ----------------------------------
%switch g.type
%case 'coher',
%   R = R ./ cumulXY;
%   if ~isnan(g.alpha) & isnan(g.rboot)
%      Rboot = Rboot ./ cumulXYboot;  
%   end;
%   if g.bootsub > 0
%      Rboottrial = Rboottrial ./ cumulXYboottrial;
%   end;
%case 'phasecoher',
%   Rn = sum(Rn, 1);
%   R = R ./ (ones(size(R,1),1)*Rn);               % coherence magnitude
%   if ~isnan(g.alpha) & isnan(g.rboot)
%      Rboot = Rboot / trials;  
%   end;
%   if g.bootsub > 0
%      Rboottrial = Rboottrial / trials;
%   end;
%end;

% compute baseline
% ----------------
mbase = mean(abs(Coher.R(:,baseln)'));     % mean baseline coherence magnitude

% plot everything
% ---------------
plotall(Coher.R, Boot.Coherboot.R, Boot.Rsignif, times, freqs, mbase, dispf, g);
Rangle = angle(Coher.R);
R = abs(Coher.R);
Rsignif = Boot.Rsignif;
Tfx = permute(Tfx.tmpall, [3 2 1]); % from [trials timesout nb_points] -> [nb_points timesout trials]
Tfy = permute(Tfy.tmpall, [3 2 1]);

return; % ***********************************************************************

% ------------------
% plotting functions
% ------------------
function plotall(R, Rboot, Rsignif, times, freqs, mbase, dispf, g) 

switch lower(g.plotphase)
case 'on',  
   switch lower(g.plotamp), 
   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.plotamp), 
   case 'on', ordinate1 = 0.1; height = 0.9;  g.plot = 1;
   case 'off', g.plot = 0;
   end;     
end; 

% compute angles
% --------------
Rangle = angle(R);
R = abs(R);
% if ~isnan(g.baseline)
% 	R = R - repmat(mbase',[1 g.timesout]); % remove baseline mean
% end;
Rraw =R; % raw coherence values

if g.plot
   fprintf('\nNow plotting...\n');
   set(gcf,'DefaultAxesFontSize',g.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')
end;

switch lower(g.plotamp)
case 'on' 
   %
   % Image the coherence [% perturbations] 
   %
   RR = R;
   if ~isnan(g.alpha) % zero out (and 'green out') nonsignif. R values
      RR(find(RR < repmat(Rboot(:),[1 g.timesout]))) = 0;
      Rraw(find(repmat(Rsignif(:),[1,size(Rraw,2)])>=Rraw))=0;
   end
   
   if g.cmax == 0
      coh_caxis = max(max(R(dispf,:)))*[-1 1];
   else
      coh_caxis = g.cmax*[-1 1];
   end
   
   h(6) = axes('Units','Normalized', 'Position',[.1 ordinate1 .8 height].*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
   if ~isempty(g.maxamp)
	   caxis([-g.maxamp g.maxamp]);
   end;
   tmpscale = caxis;
   
   hold on
   plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth)
   for i=1:length(g.marktimes)
      plot([g.marktimes(i) g.marktimes(i)],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth);
   end;
   hold off
   set(h(6),'YTickLabel',[],'YTick',[])
   set(h(6),'XTickLabel',[],'XTick',[])
   %title('Event-Related Coherence')
   
   h(8) = axes('Position',[.95 ordinate1 .05 height].*s+q);
   cbar(h(8),151:300, [0 tmpscale(2)]); % use only positive colors (gyorv) 
   
   %
   % Plot delta-mean min and max coherence at each time point on bottom of image
   %
   h(10) = axes('Units','Normalized','Position',[.1 ordinate1-0.1 .8 .1].*s+q); % plot marginal means below
   Emax = max(R(dispf,:)); % mean coherence at each time point
   Emin = min(R(dispf,:)); % mean coherence at each time point
   plot(times,Emin, times, Emax, 'LineWidth',g.linewidth); hold on;
   plot([times(1) times(length(times))],[0 0],'LineWidth',0.7);
   plot([0 0],[-500 500],'--m','LineWidth',g.linewidth);
   for i=1:length(g.marktimes)
       plot([g.marktimes(i) g.marktimes(i)],[-500 500],'--m','LineWidth',g.linewidth);
   end;
   if ~isnan(g.alpha) & strcmp(g.boottype, 'trials') 
       % plot bootstrap significance limits (base mean +/-)
      plot(times,mean(Rboot(dispf,:)),'g','LineWidth',g.linewidth); hold on;
      plot(times,mean(Rsignif(dispf,:)),'k:','LineWidth',g.linewidth);
      axis([min(times) max(times) 0 max([Emax(:)' Rsignif(:)'])*1.2])
   else
      axis([min(times) max(times) 0 max(Emax)*1.2])
   end;
   tick = get(h(10),'YTick');
   set(h(10),'YTick',[tick(1) ; tick(length(tick))])
   set(h(10),'YAxisLocation','right')
   xlabel('Time (ms)')
   ylabel('coh.')
   
   %
   % Plot mean baseline coherence at each freq on left side of image
   %
   
   h(11) = axes('Units','Normalized','Position',[0 ordinate1 .1 height].*s+q); % plot mean spectrum
   E = abs(mbase(dispf)); % baseline mean coherence at each frequency
   plot(freqs(dispf),E,'LineWidth',g.linewidth); % plot mbase
   if ~isnan(g.alpha) % plot bootstrap significance limits (base mean +/-)
      hold on
      % plot(freqs(dispf),Rboot(:,dispf)+[E;E],'g','LineWidth',g.linewidth);
      plot(freqs(dispf),mean(Rboot  (dispf,:),2),'g','LineWidth',g.linewidth);
      plot(freqs(dispf),mean(Rsignif(dispf,:),2),'k:','LineWidth',g.linewidth);
      axis([freqs(1) freqs(max(dispf)) 0 max([E Rsignif(:)'])*1.2]);
   else             % plot marginal mean coherence only
      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.')
end;

switch lower(g.plotphase)
case 'on'
   %
   % Plot coherence phase lags in bottom panel
   %
   h(13) = axes('Units','Normalized','Position',[.1 ordinate2 .8 height].*s+q);
   if strcmp(g.angleunit,'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',g.linewidth); % zero-time line
   for i=1:length(g.marktimes)
      plot([g.marktimes(i) g.marktimes(i)],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth);
   end;
   
   ylabel('Freq. (Hz)')
   xlabel('Time (ms)')