crossf.m

Cross frecuency relations

   Tf.win = dftfilt(winsize,maxfreq/srate,cycles,padratio,cyclesfact);
   Tf.nb_points = size(Tf.win,2);
end;
Tf.tmpalltimes = zeros(Tf.nb_points, timesout);
trials = length(X)/frame;
if saveall
   Tf.tmpall     = repmat(nan,[trials timesout Tf.nb_points]);
else
   Tf.tmpall = [];
end
Tf.tmpallbool = zeros(trials,timesout);
Tf.ITCdone = 0;
if Tf.subitc
   Tf.ITC  = zeros(Tf.nb_points, timesout);
   switch Tf.type,
	   case { 'coher' 'phasecoher2' }
		Tf.ITCcumul  = zeros(Tf.nb_points, timesout);
   end;
end;

% function for itc
% ----------------
function [Tf, itcvals] = tfitc(Tf, trials, times);
Tf = tfcomp(Tf, trials, times);
switch Tf.type
   case 'coher',
      Tf.ITC(:,times)      = Tf.ITC(:,times) + Tf.tmpalltimes; % complex coher.
      Tf.ITCcumul(:,times) = Tf.ITCcumul(:,times)+abs(Tf.tmpalltimes).^2;
   case 'phasecoher2',
      Tf.ITC(:,times)      = Tf.ITC(:,times) + Tf.tmpalltimes; % complex coher.
      Tf.ITCcumul(:,times) = Tf.ITCcumul(:,times)+abs(Tf.tmpalltimes);
   case 'phasecoher',
      Tf.ITC(:,times)      = Tf.ITC(:,times) + Tf.tmpalltimes ./ abs(Tf.tmpalltimes); 
                                                            % complex coher.
end % ~any(isnan())
return;

function [Tf, itcvals] = tfitcpost(Tf, trials);
switch Tf.type
   case 'coher',       Tf.ITC = Tf.ITC ./ sqrt(trials * Tf.ITCcumul);
   case 'phasecoher2', Tf.ITC = Tf.ITC ./ Tf.ITCcumul;
   case 'phasecoher',  Tf.ITC = Tf.ITC / trials; % complex coher.
end % ~any(isnan())

if Tf.saveall
  Tf.ITC = transpose(Tf.ITC); % do not use ' otherwise conjugate

	%imagesc(abs(Tf.ITC)); colorbar; figure;
	%squeeze(Tf.tmpall(1,1,1:Tf.nb_points))
	%squeeze(Tf.ITC   (1,1,1:Tf.nb_points))
	%Tf.ITC = shiftdim(Tf.ITC, -1);

	Tf.ITC = repmat(shiftdim(Tf.ITC, -1), [trials 1 1]);
	Tf.tmpall = (Tf.tmpall - abs(Tf.tmpall) .* Tf.ITC) ./ abs(Tf.tmpall);

  %	for index = 1:trials
  %		imagesc(squeeze(abs(Tf.tmpall(index,:,:)))); drawnow; figure;
  %		Tf.tmpall(index,:,:) = (Tf.tmpall(index,:,:) - Tf.tmpall(index,:,:) .* Tf.ITC)./Tf.tmpall(index,:,:);
  %		imagesc(squeeze(abs(Tf.tmpall(index,:,:)))); drawnow;
  %		subplot(10,10, index); imagesc(squeeze(abs(Tf.tmpall(index,:,:)))); caxis([0 1]); drawnow;
  %	end;
  %	squeeze(Tf.tmpall(1,1,1:Tf.nb_points))
  %	figure; axcopy;

end;
Tf.ITCdone = 1;
return;

% function for time freq decomposition
% ------------------------------------
function [Tf, tmpX] = tfcomp(Tf, trials, times);
% tf is an structure containing all the information about the decomposition
for trial = trials
   for index = times
      if ~Tf.tmpallbool(trial, index) % already computed
         tmpX = Tf.X([1:Tf.winsize]+floor((index-1)*Tf.stp)+(trial-1)*Tf.frame);
         
         if ~any(isnan(tmpX)) % perform the decomposition
            tmpX = tmpX - mean(tmpX);
            switch Tf.detret, case 'on', 
               tmpX = detrend(tmpX); 
            end;
            
            if Tf.cycles == 0 % use FFTs
               tmpX = Tf.win .* tmpX(:);
               tmpX = fft(tmpX,Tf.padratio*Tf.winsize);
               tmpX = tmpX(2:Tf.padratio*Tf.winsize/2+1);
            else 
               tmpX = transpose(Tf.win) * tmpX(:);
            end
         else
            tmpX = NaN;
         end;
         if Tf.ITCdone
            tmpX = (tmpX - abs(tmpX) .* Tf.ITC(:,index)) ./ abs(tmpX);
         end;
         Tf.tmpalltimes(:,index) = tmpX;
         if Tf.saveall
            Tf.tmpall(trial, index,:) = tmpX;
            Tf.tmpallbool(trial, index) = 1;
         end
	  else
		  Tf.tmpalltimes(:,index) = Tf.tmpall(trial, index,:);
     end;
   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    COHERENCE    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for coherence initialisation
% -------------------------------------
function Coher = coherinit(nb_points, trials, timesout, type);
Coher.R  = zeros(nb_points,timesout);       % mean coherence
% Coher.RR = repmat(nan,nb_points,timesout); % initialize with nans
Coher.type = type;
Coher.Rn=zeros(trials,timesout);
switch type
 case 'coher',
  Coher.cumulX = zeros(nb_points,timesout);
  Coher.cumulY = zeros(nb_points,timesout);
 case 'phasecoher2',
  Coher.cumul  = zeros(nb_points,timesout);
end;

% function for coherence calculation
% -------------------------------------
% function Coher = cohercomparray(Coher, tmpX, tmpY, trial);
% switch Coher.type
%   case 'coher',
%      Coher.R = Coher.R + tmpX.*conj(tmpY); % complex coher.
%      Coher.cumulXY = Coher.cumulXY + abs(tmpX).*abs(tmpY);
%   case 'phasecoher',
%      Coher.R = Coher.R + tmpX.*conj(tmpY) ./ (abs(tmpX).*abs(tmpY)); % complex coher.
%      Coher.Rn(trial,:) = 1;
% end % ~any(isnan())

function [Coher,tmptrialcoh] = cohercomp(Coher, tmpX, tmpY, trial, time);
tmptrialcoh = tmpX.*conj(tmpY);
switch Coher.type
   case 'coher',
      Coher.R(:,time) = Coher.R(:,time) + tmptrialcoh; % complex coher.
      Coher.cumulX(:,time) = Coher.cumulX(:,time) + abs(tmpX).^2;
      Coher.cumulY(:,time) = Coher.cumulY(:,time) + abs(tmpY).^2;
 case 'phasecoher2',
      Coher.R(:,time) = Coher.R(:,time) + tmptrialcoh; % complex coher.
      Coher.cumul(:,time) = Coher.cumul(:,time) + abs(tmptrialcoh);
   case 'phasecoher',
      Coher.R(:,time) = Coher.R(:,time) + tmptrialcoh ./ abs(tmptrialcoh); % complex coher.
	  %figure; imagesc(abs(tmpX.*conj(tmpY) ./ (abs(tmpX).*abs(tmpY))));
      Coher.Rn(trial,time) = Coher.Rn(trial,time)+1;
end % ~any(isnan())

% function for post coherence calculation
% ---------------------------------------
function Coher = cohercomppost(Coher, trials);
switch Coher.type
 case 'coher',
   Coher.R = Coher.R ./ sqrt(Coher.cumulX) ./ sqrt(Coher.cumulY);
 case 'phasecoher2',
   Coher.R = Coher.R ./ Coher.cumul;
 case 'phasecoher',
   Coher.Rn = sum(Coher.Rn, 1);
   Coher.R  = Coher.R ./ (ones(size(Coher.R,1),1)*Coher.Rn); % coherence magnitude
end;

% function for 2 conditions coherence calculation
% -----------------------------------------------
function [coherimage, coherimage1, coherimage2] = coher2conddiff( allsavedcoher, alltrials, cond1trials, type, tfx, tfy);
	t1s = alltrials(1:cond1trials);
	t2s = alltrials(cond1trials+1:end);
	switch type
	 case 'coher',
	  coherimage1 = sum(allsavedcoher(:,:,t1s),3) ./ sqrt(sum(tfx(:,:,t1s),3)) ./ sqrt(sum(tfy(:,:,t1s),3));
	  coherimage2 = sum(allsavedcoher(:,:,t2s),3) ./ sqrt(sum(tfx(:,:,t2s),3)) ./ sqrt(sum(tfy(:,:,t1s),3));
	 case 'phasecoher2',
	  coherimage1 = sum(allsavedcoher(:,:,t1s),3) ./ sum(abs(allsavedcoher(:,:,t1s)),3);
	  coherimage2 = sum(allsavedcoher(:,:,t2s),3) ./ sum(abs(allsavedcoher(:,:,t2s)),3);
	 case 'phasecoher',
	  coherimage1 = sum(allsavedcoher(:,:,t1s),3) / cond1trials;
	  coherimage2 = sum(allsavedcoher(:,:,t2s),3) / (size(allsavedcoher,3)-cond1trials);
	end;
	coherimage = coherimage2 - coherimage1;
	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BOOTSTRAP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for bootstrap initialisation
% -------------------------------------
function Boot = bootinit(Coherboot, nb_points, timesout, naccu, baselength, baseboot, boottype, alpha, rboot);
Boot.Rboot       = zeros(nb_points,naccu);  % summed bootstrap coher
Boot.boottype    = boottype;
Boot.baselength  = baselength;
Boot.baseboot    = baseboot;
Boot.Coherboot   = Coherboot;
Boot.naccu       = naccu;
Boot.alpha       = alpha;
Boot.rboot       = rboot;

% function for bootstrap computation
% ----------------------------------
function Boot = bootcomp(Boot, Rn, tmpalltimesx, tmpalltimesy);
if ~isnan(Boot.alpha) & isnan(Boot.rboot)
   if strcmp(Boot.boottype, 'times') % get g.naccu bootstrap estimates for each trial
      goodbasewins = find(Rn==1);
      if Boot.baseboot % use baseline windows only
         goodbasewins = find(goodbasewins<=Boot.baselength); 
      end
      ngdbasewins = length(goodbasewins);
      j=1;
      tmpsX = zeros(size(tmpalltimesx,1), Boot.naccu);
      tmpsY = zeros(size(tmpalltimesx,1), Boot.naccu);
      if ngdbasewins > 1
         while j<=Boot.naccu
            s = ceil(rand([1 2])*ngdbasewins); % random ints [1,g.timesout]
            s = goodbasewins(s);
            if ~any(isnan(tmpalltimesx(:,s(1)))) & ~any(isnan(tmpalltimesy(:,s(2))))
               tmpsX(:,j) = tmpalltimesx(:,s(1));
               tmpsY(:,j) = tmpalltimesy(:,s(2));
               j = j+1;
            end
         end
         Boot.Coherboot = cohercomp(Boot.Coherboot, tmpsX, tmpsY, 1, 1:Boot.naccu);
      end;
   end
end;

% handle other trial bootstrap types
% ----------------------------------
function [Boot, Rbootout] = bootcomppost(Boot, allRn, alltmpsX, alltmpsY);
trials    = size(alltmpsX, 1);
times     = size(alltmpsX, 2);
nb_points = size(alltmpsX, 3);
if ~isnan(Boot.alpha) & isnan(Boot.rboot)
   if strcmp(Boot.boottype, 'trials') % get g.naccu bootstrap estimates for each trial
      fprintf('\nProcessing trial bootstrap (of %d):',times(end));
      tmpsX = zeros(size(alltmpsX,3), Boot.naccu);
      tmpsY = zeros(size(alltmpsY,3), Boot.naccu );
      Boot.fullcoherboot = zeros(nb_points, Boot.naccu, times);
      
      for index=1:times
         if rem(index,10) == 0,  fprintf(' %d',index); end
         if rem(index,120) == 0, fprintf('\n'); end
         for allt=1:trials
            j=1;
            while j<=Boot.naccu
               t = ceil(rand([1 2])*trials); % random ints [1,g.timesout]
               if (allRn(t(1),index) == 1) & (allRn(t(2),index) == 1)
                  tmpsX(:,j) = squeeze(alltmpsX(t(1),index,:));
                  tmpsY(:,j) = squeeze(alltmpsY(t(2),index,:));
                  j = j+1;
               end
            end
            Boot.Coherboot = cohercomp(Boot.Coherboot, tmpsX, tmpsY, 1, 1:Boot.naccu);
         end;
         Boot.Coherboot = cohercomppost(Boot.Coherboot);  % CHECK IF NECSSARY FOR ALL BOOT TYPE
         Boot.fullcoherboot(:,:,index) = Boot.Coherboot.R; 
         Boot.Coherboot = coherinit(nb_points, trials, Boot.naccu, Boot.Coherboot.type);
      end;
      Boot.Coherboot.R = Boot.fullcoherboot;
      Boot = rmfield(Boot, 'fullcoherboot');
   elseif strcmp(Boot.boottype, 'timestrials') % handle timestrials bootstrap
      fprintf('\nProcessing time and trial bootstrap (of %d):',trials);
      tmpsX = zeros(size(alltmpsX,3), Boot.naccu);
      tmpsY = zeros(size(alltmpsY,3), Boot.naccu );
      for allt=1:trials
         if rem(allt,10) == 0,  fprintf(' %d',allt); end
         if rem(allt,120) == 0, fprintf('\n'); end
         j=1;
         while j<=Boot.naccu
            t = ceil(rand([1 2])*trials); % random ints [1,g.timesout]
            goodbasewins = find((allRn(t(1),:) & allRn(t(2),:)) ==1);
            if Boot.baseboot % use baseline windows only
               goodbasewins = find(goodbasewins<=baselength); 
            end
            ngdbasewins = length(goodbasewins);
            
            if ngdbasewins>1
               s = ceil(rand([1 2])*ngdbasewins); % random ints [1,g.timesout]
               s=goodbasewins(s);
               
               if all(allRn(t(1),s(1)) == 1) & all(allRn(t(2),s(2)) == 1)
                  tmpsX(:,j) = squeeze(alltmpsX(t(1),s(1),:));
                  tmpsY(:,j) = squeeze(alltmpsY(t(2),s(2),:));
                  j = j+1;
               end
            end
         end
         Boot.Coherboot = cohercomp(Boot.Coherboot, tmpsX, tmpsY, 1, 1:Boot.naccu);
      end
      Boot.Coherboot = cohercomppost(Boot.Coherboot);
   elseif strcmp(Boot.boottype, 'times') % boottype is 'times'
      Boot.Coherboot = cohercomppost(Boot.Coherboot);
   end;
end;

% test if precomputed
if ~isnan(Boot.alpha) & isnan(Boot.rboot) % if bootstrap analysis included . . .
   % 'boottype'='times' or 'timestrials', size(R)=nb_points*naccu
   % 'boottype'='trials',                 size(R)=nb_points*naccu*times
   Boot.Coherboot.R = abs (Boot.Coherboot.R);
   Boot.Coherboot.R = sort(Boot.Coherboot.R,2);

   % compute bootstrap significance level
   i = round(Boot.naccu*Boot.alpha);
   Boot.Rsignif = mean(Boot.Coherboot.R(:,Boot.naccu-i+1:Boot.naccu),2); % significance levels for Rraw
   Boot.Coherboot.R = squeeze(mean(Boot.Coherboot.R(:,Boot.naccu-i+1:Boot.naccu),2));
   if size(Boot.Coherboot.R, 2) == 1
	   Rbootout(:,2) = Boot.Coherboot.R;
   else
	   Rbootout(:,:,2) = Boot.Coherboot.R;
   end;
   % BEFORE
   %Rboot = [mean(Rboot(1:i,:)) ; mean(Rboot(g.naccu-i+1:g.naccu,:))];
elseif ~isnan(Boot.rboot)
	Boot.Coherboot.R = Boot.rboot;
	Boot.Rsignif     = Boot.rboot
	Rbootout         = Boot.rboot;
else 
	Boot.Coherboot.R = [];
	Boot.Rsignif     = [];
	Rbootout         = [];
end % NOTE: above, mean ?????

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