crossf.m

Cross frecuency relations

%
% Revision 1.18  2002/04/24 21:45:18  scott
% topovec bug -sm
%
% Revision 1.17  2002/04/24 21:43:00  scott
% editing topovec code -sm
%
% Revision 1.16  2002/04/24 21:02:28  scott
% added topoplots of two heads -sm
%
% Revision 1.15  2002/04/24 02:43:18  arno
% debugging amplitude coherence
%
% Revision 1.14  2002/04/20 00:53:14  arno
% restorings some outputs options
%
% Revision 1.13  2002/04/19 23:20:23  arno
% changing trial bootstrap, not optimal, waiting for further inputs
%
% Revision 1.12  2002/04/19 19:46:28  arno
% crossf with new trial coherence bootstrap (minus mean)
%
% Revision 1.11  2002/04/12 18:10:55  scott
% added note
%
% Revision 1.10  2002/04/12 01:30:43  arno
% compatibility for returning frequencies with timef
%
% Revision 1.9  2002/04/12 01:13:40  arno
% debuging no ploting option
%
% Revision 1.8  2002/04/12 01:08:13  arno
% change plotamps to plotamp in help message
%
% Revision 1.7  2002/04/12 00:41:37  arno
% programming baseboot
%
% Revision 1.6  2002/04/11 02:39:34  arno
% updated header message
%
% Revision 1.5  2002/04/10 01:29:45  arno
% adding vert optional input
%
% Revision 1.4  2002/04/09 19:36:38  arno
% corrected bootstrap optional input
%
% Revision 1.3  2002/04/09 18:59:06  arno
% corrected typo in header that made the function to crash
%
% Revision 1.2  2002/04/07 02:24:36  scott
% worked on hlpe message, changed some defaults -sm
%
% Revision 1.1  2002/04/05 17:36:45  jorn
% Initial revision
%

% 11-20-98 defined g.linewidth constant -sm
% 04-01-99 made number of frequencies consistent -se
% 06-29-99 fixed constant-Q freq indexing -se
% 08-13-99 added cohangle plotting -sm
% 08-20-99 made bootstrap more efficient -sm
% 08-24-99 allow nan values introduced by possible eventlock() preproc. -sm
% 03-05-2007 eventlock.m deprecated to eegalign.m. -tf
% 03-16-00 added lead/lag interpretation to help msg - sm & eric visser
% 03-16-00 added axcopy() feature -sm & tpj
% 04-20-00 fixed Rangle sign for wavelets, added verts array -sm
% 01-22-01 corrected help msg when nargin<2 -sm & arno delorme
% 01-25-02 reformated help & license, added links -ad 
% 03-09-02 function restructuration -ad
%  add 'key', val arguments (+ external baseboot, baseline, color axis, angleunit...)
%  add detrending (across time and trials) + 'coher' option for amplitude coherence
%  significance only if alpha is given, ploting options in 'plotamp' and 'plotphase'
% 03-16-02 timeout automatically adjusted if too high -ad 
% 04-03-02 added new options for bootstrap -ad 

% Note: 3 "objects" (Tf, Coher and Boot) are handled by specific functions under Matlab
%    (Tf) function Tf = tfinit(...) - create object Time Frequency (Tf) associated with some data
%    (Tf) function [Tf, itcvals] = tfitc(...) - compute itc for the selected data
%    (Tf) function [Tf, itcvals] = tfitcpost(Tf, trials) - itc normlisation 
%    (Tf) function [Tf, tmpX] = tfcomp(Tf, trials, times) - compute time freq. decomposition
%    (Coher) function Coher = coherinit(...) - initialize coherence object
%    (Coher) function Coher = cohercomp(Coher, tmpX, tmpY, trial, time) - compute coherence
%    (Coher) function Coher = cohercomppost(Coher, trials) - coherence normalization
%    (Boot) function Boot = bootinit(...) - intialize bootstrap object
%    (Boot) function Boot = bootcomp(...) - compute bootstrap
%    (Boot) function [Boot, Rbootout] = bootcomppost(...) - bootstrap normalization
% and by real objects under C++ (C++ code, incomplete)

function [R,mbase,times,freqs,Rbootout,Rangle, trialcoher, Tfx, Tfy] = crossf(X, Y, frame, tlimits, Fs, varwin, varargin)

%varwin,winsize,nwin,oversmp,maxfreq,alpha,verts,caxmax)

% ------------------------
% Commandline arg defaults:
% ------------------------
DEFAULT_ANGLEUNIT = 'deg'; % angle plotting units - 'rad', 'ms', or 'deg'
DEFAULT_EPOCH	= 750;			% Frames per epoch
DEFAULT_TIMELIM = [-1000 2000];	% Time range of epochs (ms)
DEFAULT_FS		= 250;			% Sampling frequency (Hz)
DEFAULT_NWIN	= 200;			% Number of windows = horizontal resolution
DEFAULT_VARWIN	= 0;			% Fixed window length or base on cycles.

% =0: fix window length to nwin
% >0: set window length equal varwin cycles
%     bounded above by winsize, also determines
%     the min. freq. to be computed.

DEFAULT_OVERSMP	= 2;			% Number of times to oversample = vertical resolution
DEFAULT_MAXFREQ = 50;			% Maximum frequency to display (Hz)
DEFAULT_TITLE	= 'Event-Related Coherence';			% Figure title
DEFAULT_ALPHA   = NaN;			% Default two-sided significance probability threshold

if (nargin < 2)
   help crossf
   return
end

if ~iscell(X)
	if (min(size(X))~=1 | length(X)<2)
		fprintf('crossf(): X must be a row or column vector.\n');
		return
	elseif (min(size(Y))~=1 | length(Y)<2)
		fprintf('crossf(): Y must be a row or column vector.\n');
		return
	elseif (length(X) ~= length(Y))
		fprintf('crossf(): X and Y must have same length.\n');
		return
	end
end;

if (nargin < 3)
   frame = DEFAULT_EPOCH;
elseif (~isnumeric(frame) | length(frame)~=1 | frame~=round(frame))
   fprintf('crossf(): Value of frames must be an integer.\n');
   return
elseif (frame <= 0)
   fprintf('crossf(): Value of frames must be positive.\n');
   return
elseif ~iscell(X) & (rem(length(X),frame) ~= 0)
   fprintf('crossf(): Length of data vectors must be divisible by frames.\n');
   return
end

if (nargin < 4)
   tlimits = DEFAULT_TIMELIM;
elseif (~isnumeric(tlimits) | sum(size(tlimits))~=3)
   error('crossf(): Value of tlimits must be a vector containing two numbers.');
elseif (tlimits(1) >= tlimits(2))
   error('crossf(): tlimits interval must be [min,max].');
end

if (nargin < 5)
   Fs = DEFAULT_FS;
elseif (~isnumeric(Fs) | length(Fs)~=1)
   error('crossf(): Value of srate must be a number.');
elseif (Fs <= 0)
   error('crossf(): Value of srate must be positive.');
end

if (nargin < 6)
   varwin = DEFAULT_VARWIN;
elseif (~isnumeric(varwin) | length(varwin)>2)
   error('crossf(): Value of cycles must be a number or a (1,2) vector.');
elseif (varwin < 0)
   error('crossf(): Value of cycles must be either zero or positive.');
end

% consider structure for these arguments
% --------------------------------------
vararginori = varargin;
for index=1:length(varargin)
	if iscell(varargin{index}), varargin{index} = { varargin{index} }; end;
end;
if ~isempty(varargin)
   try, g = struct(varargin{:}); 
   catch, error('Argument error in the {''param'', value} sequence'); end; 
else 
	g = [];
end;

try, g.shuffle;    catch, g.shuffle = 0; end;
try, g.title;      catch, g.title = DEFAULT_TITLE; end;
try, g.winsize;    catch, g.winsize = max(pow2(nextpow2(frame)-3),4); end;
try, g.pad;        catch, g.pad = max(pow2(nextpow2(g.winsize)),4); end;
try, g.timesout;   catch, g.timesout = DEFAULT_NWIN; end;
try, g.padratio;   catch, g.padratio = DEFAULT_OVERSMP; end;
try, g.maxfreq;    catch, g.maxfreq = DEFAULT_MAXFREQ; end;
try, g.topovec;    catch, g.topovec = []; end;
try, g.elocs;      catch, g.elocs = ''; end;
try, g.alpha;      catch, g.alpha = DEFAULT_ALPHA; end;  
try, g.marktimes;  catch, g.marktimes = []; end; % default no vertical lines
try, g.marktimes = g.vert;       catch, g.vert = []; end; % default no vertical lines
try, g.powbase;    catch, g.powbase = nan; end;
try, g.rboot;      catch, g.rboot = nan; end;
try, g.plotamp;    catch, g.plotamp = 'on'; end;
try, g.plotphase;  catch, g.plotphase  = 'on'; end;
try, g.plotbootsub;  catch, g.plotbootsub  = 'on'; end;
try, g.detrep;     catch, g.detrep = 'off'; end;
try, g.detret;     catch, g.detret = 'off'; end;
try, g.baseline;   catch, g.baseline = NaN; end;
try, g.baseboot;   catch, g.baseboot = 0; end;
try, g.linewidth;  catch, g.linewidth = 2; end;
try, g.naccu;      catch, g.naccu = 200; end;
try, g.angleunit;  catch, g.angleunit = DEFAULT_ANGLEUNIT; end;
try, g.cmax;       catch, g.cmax = 0; end; % 0=use data limits
try, g.type;       catch, g.type = 'phasecoher'; end; 
try, g.boottype;   catch, g.boottype = 'times'; end; 
try, g.subitc;     catch, g.subitc = 'off'; end;
try, g.memory;     catch, g.memory = 'high'; end;
try, g.compute;    catch, g.compute = 'matlab'; end;
try, g.maxamp;     catch, g.maxamp = []; end;
try, g.savecoher;  catch, g.savecoher = 0; end;
try, g.noinput;    catch, g.noinput = 'no'; end;
try, g.chaninfo;   catch, g.chaninfo = []; end;

allfields = fieldnames(g);
for index = 1:length(allfields)
	switch allfields{index}
	 case { 'shuffle' 'title' 'winsize' 'pad' 'timesout' 'padratio' 'maxfreq' 'topovec' 'elocs' 'alpha' ...
		  'marktimes' 'vert' 'powbase' 'rboot' 'plotamp' 'plotphase' 'plotbootsub' 'detrep' 'detret' ...
		  'baseline' 'baseboot' 'linewidth' 'naccu' 'angleunit' 'cmax' 'type' 'boottype' 'subitc' ...
		  'memory' 'compute' 'maxamp' 'savecoher' 'noinput' 'chaninfo' };
	  case {'plotersp' 'plotitc' }, disp(['crossf warning: timef option ''' allfields{index} ''' ignored']);
	 otherwise disp(['crossf error: unrecognized option ''' allfields{index} '''']); beep; return;
	end;
end;

g.tlimits = tlimits;
g.frame   = frame;
g.srate   = Fs;
g.cycles  = varwin(1);
if length(varwin)>1
	g.cyclesfact = varwin(2);
else 
	g.cyclesfact = 1;
end;
g.type       = lower(g.type);
g.boottype   = lower(g.boottype);
g.detrep     = lower(g.detrep);
g.detret     = lower(g.detret);
g.plotphase  = lower(g.plotphase);
g.plotbootsub = lower(g.plotbootsub);
g.subitc     = lower(g.subitc);
g.plotamp    = lower(g.plotamp);
g.shuffle    = lower(g.shuffle);
g.compute    = lower(g.compute);
g.AXES_FONT  = 10;
g.TITLE_FONT = 14;

% testing arguments consistency
% -----------------------------
if (~ischar(g.title))
   error('Title must be a string.');
end

if (~isnumeric(g.winsize) | length(g.winsize)~=1 | g.winsize~=round(g.winsize))
   error('Value of winsize must be an integer number.');
elseif (g.winsize <= 0)
   error('Value of winsize must be positive.');
elseif (g.cycles == 0 & pow2(nextpow2(g.winsize)) ~= g.winsize)
   error('Value of winsize must be an integer power of two [1,2,4,8,16,...]');
elseif (g.winsize > g.frame)
   error('Value of winsize must be less than frame length.');
end

if (~isnumeric(g.timesout) | length(g.timesout)~=1 | g.timesout~=round(g.timesout))
   error('Value of timesout must be an integer number.');
elseif (g.timesout <= 0)
   error('Value of timesout must be positive.');
end
if (g.timesout > g.frame-g.winsize)
   g.timesout = g.frame-g.winsize;
   disp(['Value of timesout must be <= frame-winsize, timeout adjusted to ' int2str(g.timesout) ]);
end

if (~isnumeric(g.padratio) | length(g.padratio)~=1 | g.padratio~=round(g.padratio))
   error('Value of padratio must be an integer.');
elseif (g.padratio <= 0)
   error('Value of padratio must be positive.');
elseif (pow2(nextpow2(g.padratio)) ~= g.padratio)
   error('Value of padratio must be an integer power of two [1,2,4,8,16,...]');
end

if (~isnumeric(g.maxfreq) | length(g.maxfreq)~=1)
   error('Value of g.maxfreq must be a number.');
elseif (g.maxfreq <= 0)
   error('Value of g.maxfreq must be positive.');
elseif (g.maxfreq > Fs/2)
   fprintf('Warning: input value of g.maxfreq larger that Nyquist frequency %3.4 Hz\n\n',Fs/2);
end

if isempty(g.topovec)
   g.topovec = [];
elseif min(size(g.topovec))==1
   g.topovec = g.topovec(:);
   if size(g.topovec,1)~=2
      error('topovec must be a row or column vector.');
   end
end;

if isempty(g.elocs)
   g.elocs = '';
elseif (~ischar(g.elocs)) & ~isstruct(g.elocs)
   error('Channel location file must be a valid text file.');
end

if (~isnumeric(g.alpha) | length(g.alpha)~=1)
   error('timef(): Value of g.alpha must be a number.\n');
elseif (round(g.naccu*g.alpha) < 2)
   fprintf('Value of g.alpha is out of the normal range [%g,0.5]\n',2/g.naccu);
   g.naccu = round(2/g.alpha);
   fprintf('  Increasing the number of bootstrap iterations to %d\n',g.naccu);
end
if g.alpha>0.5 | g.alpha<=0
   error('Value of g.alpha is out of the allowed range (0.00,0.5).');
end
if ~isnan(g.alpha)
   if g.baseboot > 0
      fprintf('Bootstrap analysis will use data in baseline (pre-0) subwindows only.\n')
   else
      fprintf('Bootstrap analysis will use data in all subwindows.\n')
   end
end
switch g.angleunit
   case { 'rad', 'ms', 'deg' },;
   otherwise error('Angleunit must be either ''rad'', ''deg'', or ''ms''');
end;    
switch g.type
   case { 'coher', 'phasecoher' 'phasecoher2' },;
   otherwise error('Type must be either ''coher'' or ''phasecoher''');
end;    
switch g.boottype
   case { 'times' 'timestrials' 'trials'},;
   otherwise error('Boot type must be either ''times'', ''trials'' or ''timestrials''');