crossf.m

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

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_ANGLEUNITS; 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;

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' };
	  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: value of g.maxfreq greater that Nyquist rate\n\n');
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 { 'ms', 'deg' },;
otherwise error('Angleunit must be either ''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''');
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 part
%%%%%%%%%%%%%%%%%%%%%%%%%%%
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: 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');
	fprintf('      (the ''coher'' options takes 3 times more memory than other options)\n');
	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/180*pi);
	
	% 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 ''noinput'' option to disable 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/180*pi);

	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 x timesout x 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):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):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;
		  alltfx = alltfx.^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 FOR SEVERAL 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 = [];