% 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);      % Note 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);% --------------------------------------% Convert output Rangle to degrees or ms - Disabled to keep original default: radians output% --------------------------------------% Rangle = angle(Coher.R); % returns radians% if strcmp(g.angleunit,'ms')  % convert to ms%    Rangle = (Rangle/(2*pi)).*repmat(1000./freqs(dispf)',1,length(times)); % elseif strcmp(g.angleunit,'deg')  % convert to deg%    Rangle = Rangle*180/pi; % convert to degrees% else % angleunit is 'rad'%    % Rangle = Rangle;% end% Rangle(find(Rraw==0)) = 0; % mask for significance - set angle at non-signif coher points to 0R = abs(Coher.R);Rsignif = Boot.Rsignif;Tfx = permute(Tfx.tmpall, [3 2 1]); % from [trials timesout nb_points] %                                     to   [nb_points timesout trials]Tfy = permute(Tfy.tmpall, [3 2 1]);return; % end crossf() *************************************************%% crossf() 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 cross-spectral angles% -----------------------------Rangle = angle(R); % returns radians%% Optionally convert Rangle to degrees or ms% ------------------------------------------if strcmp(g.angleunit,'ms')  % convert to ms   Rangle = (Rangle/(2*pi)).*repmat(1000./freqs(dispf)',1,length(times));    maxangle = max(max(abs(Rangle)));elseif strcmp(g.angleunit,'deg')  % convert to degrees   Rangle = Rangle*180/pi; % convert to degrees   maxangle = 180; % use full-cycle plotting else   maxangle = pi;  % radiansendR = abs(R);% if ~isnan(g.baseline)% 	R = R - repmat(mbase',[1 g.timesout]); % remove baseline mean% end;Rraw = R; % raw coherence (e.g., coherency) magnitude values outputif 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);   Rangle(find(Rraw==0)) = 0; % when plotting, mask for significance                               % = 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)')      h(14)=axes('Position',[.95 ordinate2 .05 height].*s+q);   cbar(h(14),0,[-maxangle maxangle]); % two-sided colorbarendif g.plot	try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end;    if (length(g.title) > 0) % plot title        if h(6) ~= 0, axes(h(6)); else axes(h(13)); end;        %h = subplot('Position',[0 0  1 1].*s+q, 'Visible','Off');                       %h(13) = text(-.05,1.01,g.title);        h(13) = title(g.title);        %set(h(13),'VerticalAlignment','bottom')        %set(h(13),'HorizontalAlignment','left')        set(h(13),'FontSize',g.TITLE_FONT);    end   %   %%%%%%%%%%%%%%% plot topoplot() %%%%%%%%%%%%%%%%%%%%%%%   %   if (~isempty(g.topovec))      h(15) = subplot('Position',[-.1 .43 .2 .14].*s+q);      if size(g.topovec,2) <= 2         topoplot(g.topovec(1),g.elocs,'electrodes','off', ...            'style', 'blank', 'emarkersize1chan', 10, 'chaninfo', g.chaninfo);      else         topoplot(g.topovec(1,:),g.elocs,'electrodes','off', 'chaninfo', g.chaninfo);      end;      axis('square')            h(16) = subplot('Position',[.9 .43 .2 .14].*s+q);      if size(g.topovec,2) <= 2         topoplot(g.topovec(2),g.elocs,'electrodes','off', ...            'style', 'blank', 'emarkersize1chan', 10, 'chaninfo', g.chaninfo);      else         topoplot(g.topovec(2,:),g.elocs,'electrodes','off', 'chaninfo', g.chaninfo);      end;      axis('square')   end      axcopy(gcf);end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  TIME FREQUENCY   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function for time freq initialisation% -------------------------------------function Tf = tfinit(X, timesout, winsize, ...   cycles, frame, padratio, detret, srate, maxfreq, subitc, type, cyclesfact, saveall);Tf.X         = X(:)'; % make X column vectorsTf.winsize   = winsize;Tf.cycles    = cycles;Tf.frame     = frame;Tf.padratio  = padratio;Tf.detret    = detret;Tf.stp       = (frame-winsize)/(timesout-1);Tf.subitc    = subitc; % for ITCTf.type      = type; % for ITCTf.saveall   = saveall;if (Tf.cycles == 0) %%%%%%%%%%%%%% constant window-length FFTs %%%%%%%%%%%%%%%%   % Tf.freqs = srate/winsize*[1:2/padratio:winsize]/2; % incorect for padratio > 2   Tf.freqs = linspace(0, srate/2, length([1:2/padratio:winsize])+1);   Tf.freqs = Tf.freqs(2:end);   Tf.win   = hanning(winsize);   Tf.nb_points = padratio*winsize/2;   else % %%%%%%%%%%%%%%%%%% Constant-Q (wavelet) DFTs %%%%%%%%%%%%%%%%%%%%%%%%%%%%   Tf.freqs = srate*cycles/winsize*[2:2/padratio:winsize]/2;