crossf.m

Cross frecuency relations

end;    
if (~isnumeric(g.shuffle))
   error('Shuffle argument type must be numeric');
end;
switch g.memory
   case { 'low', 'high' },;
   otherwise error('memory must be either ''low'' or ''high''');
end;
if strcmp(g.memory, 'low') & ~strcmp(g.boottype, 'times')
   error(['Bootstrap type ''' g.boottype ''' cannot be used in low memory mode']);
end;

switch g.compute
   case { 'matlab', 'c' },;
   otherwise error('compute must be either ''matlab'' or ''c''');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compare 2 conditions 
%%%%%%%%%%%%%%%%%%%%%%%%%%%
if iscell(X)
	if length(X) ~= 2 | length(Y) ~= 2
		error('crossf: to compare conditions, X and Y input must be 2-elements cell arrays');
	end;
	if ~strcmp(g.boottype, 'times')
		disp('crossf warning: The significance bootstrap type is irrelevant when comparing conditions');
	end;
	for index = 1:2:length(vararginori)
		if index<=length(vararginori) % needed: if elemenets are deleted
			%if strcmp(vararginori{index}, 'alpha'), vararginori(index:index+1) = [];
			if strcmp(vararginori{index}, 'title'), vararginori(index:index+1) = []; 
			end;
		end;
	end;
	if iscell(g.title) 
		if length(g.title) <= 2,
			g.title{3} = 'Condition 2 - condition 1';
		end;
	else
		g.title = { 'Condition 1', 'Condition 2', 'Condition 2 - condition 1' };
	end;
	
	fprintf('Running crossf on condition 1 *********************\n');
	fprintf('Note: If an out-of-memory error occurs, try reducing the\n');
	fprintf('      number of time points or number of frequencies\n');
	if ~strcmp(g.type, 'coher')
	   fprintf('Note: Type ''coher'' takes 3 times as much memory as other options!)\n');
        end
	figure; 
	subplot(1,3,1); title(g.title{1});
	if ~strcmp(g.type, 'coher')
		[R1,mbase,times,freqs,Rbootout1,Rangle1, savecoher1] = crossf(X{1}, Y{1}, ...
								frame, tlimits, Fs, varwin, 'savecoher', 1, 'title', ' ',vararginori{:});
	else
		[R1,mbase,times,freqs,Rbootout1,Rangle1, savecoher1, Tfx1, Tfy1] = crossf(X{1}, Y{1}, ...
								frame, tlimits, Fs, varwin, 'savecoher', 1,'title', ' ',vararginori{:});
	end;
	R1 = R1.*exp(j*Rangle1); % output Rangle is in radians
	
	% Asking user for memory limitations
	% if ~strcmp(g.noinput, 'yes')
	%	  tmp = whos('Tfx1');
	%	  fprintf('This function will require an additional %d bytes, do you wish\n', ...
    %     tmp.bytes*6+size(savecoher1,1)*size(savecoher1,2)*g.naccu*8);
	%	  res = input('to continue (y/n) (use the ''noinput'' option to disable this message):', 's');
	%  	if res == 'n', return; end;
	% end;

	fprintf('\nRunning crossf on condition 2 *********************\n');
	subplot(1,3,2); title(g.title{2});
	if ~strcmp(g.type, 'coher')
		  [R2,mbase,times,freqs,Rbootout2,Rangle2, savecoher2] = crossf(X{2}, Y{2}, ...
								frame, tlimits, Fs, varwin,'savecoher', 1, 'title', ' ',vararginori{:});
	else
		  [R2,mbase,times,freqs,Rbootout2,Rangle2, savecoher2, Tfx2, Tfy2] = crossf(X{2}, Y{2}, ...
								frame, tlimits, Fs, varwin,'savecoher', 1, 'title', ' ',vararginori{:});
	end;
	R2 = R2.*exp(j*Rangle2); % output Rangle is in radians

	subplot(1,3,3); title(g.title{3});
	if isnan(g.alpha)
		plotall(R2-R1, [], [], times, freqs, mbase,  find(freqs <= g.maxfreq), g);
	else 
		% accumulate coherence images (all arrays [nb_points * timesout * trials])
		% ---------------------------
		allsavedcoher = zeros(size(savecoher1,1), ...
                          size(savecoher1,2), ...
                          size(savecoher1,3)+size(savecoher2,3));
		allsavedcoher(:,:,1:size(savecoher1,3))     = savecoher1;
		allsavedcoher(:,:,size(savecoher1,3)+1:end) = savecoher2;
		clear savecoher1 savecoher2;
		
		if strcmp(g.type, 'coher')
			alltfx = zeros(size(Tfx1,1), size(Tfx2,2), size(Tfx1,3)+size(Tfx2,3));
			alltfx(:,:,1:size(Tfx1,3))     = Tfx1;
			alltfx(:,:,size(Tfx1,3)+1:end) = Tfx2;
			clear Tfx1 Tfx2;
			
			alltfy = zeros(size(Tfy1,1), size(Tfy2,2), size(Tfy1,3)+size(Tfy2,3));
			alltfy(:,:,1:size(Tfy1,3))   = Tfy1;
			alltfy(:,:,size(Tfy1,3)+1:end) = Tfy2;
			clear Tfy1 Tfy2;
		end;
		
		coherimages = zeros(size(allsavedcoher,1), size(allsavedcoher,2), g.naccu);
		cond1trials = length(X{1})/g.frame;
		cond2trials = length(X{2})/g.frame;
		alltrials = [1:cond1trials+cond2trials];
		fprintf('Accumulating bootstrap:');
		
		% preprocess data
		% ---------------
		switch g.type
		 case 'coher', % take the square of alltfx and alltfy
		  alltfx = alltfx.^2;
		  alltfy = alltfy.^2;
		 case 'phasecoher', % normalize
		  allsavedcoher = allsavedcoher ./ abs(allsavedcoher);
		 case 'phasecoher2', % don't do anything
		end;
		
		if strcmp(g.type, 'coher')
			[coherdiff coher1 coher2] = coher2conddiff( allsavedcoher, alltrials, ...
                                                        cond1trials, g.type, alltfx, alltfy);
		else
			[coherdiff coher1 coher2] = coher2conddiff( allsavedcoher, alltrials, ...
                                                        cond1trials, g.type);
		end;
		%figure; g.alpha = NaN; & to check that the new images are the same as the original
		%subplot(1,3,1); plotall(coher1, [], [], times, freqs, mbase, find(freqs <= g.maxfreq), g);
		%subplot(1,3,2); plotall(coher2, [], [], times, freqs, mbase, find(freqs <= g.maxfreq), g);
		%return;

		for index=1:g.naccu
			if rem(index,10) == 0,  fprintf(' %d',index); end
			if rem(index,120) == 0, fprintf('\n'); end
			
			if strcmp(g.type, 'coher')
				coherimages(:,:,index) = coher2conddiff( allsavedcoher, shuffle(alltrials), ...
                                                        cond1trials, g.type, alltfx, alltfy);
			else
				coherimages(:,:,index) = coher2conddiff( allsavedcoher, shuffle(alltrials), ...
                                                        cond1trials, g.type);
			end;
		end;
		fprintf('\n');

		% create articially a Bootstrap object to compute significance
		Boot = bootinit( [], size(allsavedcoher,1), g.timesout, g.naccu, 0, g.baseboot, ...
														'noboottype', g.alpha, g.rboot);
		Boot.Coherboot.R = coherimages;
		Boot = bootcomppost(Boot, [], [], []);
		g.title = '';
		plotall(coherdiff, Boot.Coherboot.R, Boot.Rsignif, times, freqs, mbase, ...
														find(freqs <= g.maxfreq), g);
	end;
	return; % ********************************** END PROCESSING TWO CONDITIONS
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% shuffle trials if necessary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if g.shuffle ~= 0
   fprintf('x and y data trials being shuffled %d times\n',g.shuffle);
   XX = reshape(X, 1, frame, length(X)/g.frame);
   YY = Y;
   X = [];
   Y = [];
   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
% case '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);
% end;
% if ~isnan(g.alpha) & isnan(g.rboot)
%    tf.tmpalltimes = repmat(nan,tfx.nb_points,g.timesout);
% end
   
% --------------------
% Display text to user
% --------------------
fprintf('\nComputing Event-Related ');
switch g.type
    case 'phasecoher',  fprintf('Phase Coherence (ITC) images for %d trials.\n',length(X)/g.frame);
    case 'phasecoher2', fprintf('Phase Coherence 2 (ITC) images for %d trials.\n',length(X)/g.frame);
    case 'coher',       fprintf('Linear Coherence (ITC) images for %d trials.\n',length(X)/g.frame);
end;
fprintf('The trial latency range is from %4.5g ms before to %4.5g ms after\n     the time-locking event.\n', 1000*g.tlimits(1),1000*g.tlimits(2));
fprintf('The frequency range displayed will be %g-%g Hz.\n',min(freqs),g.maxfreq);
if ~isnan(g.baseline)
   if length(baseln) == length(times)
      fprintf('Using the full trial latency range as baseline.\n');
   else
      fprintf('Using trial latencies from %4.5g ms to %4.5g ms as baseline.\n', 1000*g.tlimits,g.baseline);
   end;
else 
   fprintf('No baseline time range was specified.\n');	
end;
if g.cycles==0
   fprintf('The data window size will be %d samples (%g ms).\n',g.winsize,2*wintime);
   fprintf('The FFT length will be %d samples\n',g.winsize*g.padratio);
else
   fprintf('The window size will be %2.3g cycles.\n',g.cycles);
   fprintf('The maximum window size will be %d samples (%g ms).\n',g.winsize,2*wintime);
end
fprintf('The window will be applied %d times\n',g.timesout);
fprintf('     with an average step size of %2.2g samples (%2.4g ms).\n', Tfx.stp,1000*Tfx.stp/g.srate);
fprintf('Results will be 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', 
    if strcmp(g.angleunit,'deg')
       fprintf(['Coherence angles will be imaged in degrees.\n']);
    elseif strcmp(g.angleunit,'rad')
       fprintf(['Coherence angles will be imaged in radians.\n']);
    elseif strcmp(g.angleunit,'ms')
       fprintf(['Coherence angles will be imaged in ms.\n']);
    end
end;
fprintf('\nProcessing trial (of %d): ',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