erpimage.m

含有多种ICA算法的eeglab工具箱

% 08/04/99 added new times spec., 'srate' arg, made 'phase' and 'allamps'
%          work together, plot re-aligned time zeros  -sm
% 06/26/99 debugged cbar; added vert lines at aligntime to plot1erp() axes -sm
% 09/29/99 fixed srate computation from times -sm & se
% 01/18/00 output outsort without clipping -sm
% 02/29/00 added 'vert' arg, fixed xticklabels, added ampsig -sm
% 03/03/00 added 'vert' arg lines to erp/amp/coher axes -sm
% 03/17/00 added axcopy -sm & tpj
% 03/29/00 added 'topo' option -sm 
% 05/05/00 fixed y-axis label bug when time limits given -sm
% 06/01/00 added topphase arg to 'phase' option for phasemovie.m -sm
% 07/12/00 adjusted prctle()
% 07/12/00 added 'spec' option -sm
% 08/22/00 added coherfreq to limits output -sm
% 09/13/00 added hard limit (1) to maxcoh -sm
% 09/14/00 made topoplot() and psd() plots relative to gca, not gcf -sm
% 10/10/00 added NoTimeflag -sm
% 11/03/00 changed method of rejecting small amplitude trials for phase sort -sm
% 11/03/00 added number_of_trials_out option for decfactor -sm
% 11/16/00 added ampoffset to center sig lines around baseline mean amp (0) -sm
% 01/06/01 edited help message; adjusted ampsig plot limits;initialized outputs -sm
%          rm'd 'allcoher' from help message - not fully implemented -sm
% 01/09/01 documented limits arg 'bamp' (baseline amplitude) -sm
% 02/13/01 debugged use of stored baseamp, ampsig parameters -sm
% 03/28/01 made erpimage(data) possible. Debugged ampsig change limits -sm
% 08/31/01 fixed allamps bug -sm
% 09/03/01 added 'auxvar' plotting -sm
% 01-25-02 reformated help & license, added links -ad
% 02-16-02 added matrix option to arg 'vert' -sm
% 04-05-02 corrected zero alignment problem (display only) -ad
%
% Known Bugs:
% 'limits', [lotime hitime] may not work with 'erp'
% 'limits', [... loerp hierp] (still??) may leave "ghost" grey numbers 
%       on the coher axis when printed (-djpeg or -depsc)
% 'allcohers' - not fully implemented, and has been omitted from the help msg

function [data,outsort,outtrials,limits,axhndls,erp,amps,cohers,cohsig,ampsig,allamps,phaseangles,phsamp,sortidx,erpsig] = erpimage(data,sortvar,times,titl,avewidth,decfactor,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9,arg10,arg11,arg12,arg13,arg14,arg15,arg16,arg17,arg18,arg19,arg20,arg21,arg22,arg23,arg24,arg25,arg26)

%
%%%%%%%%%%%%%%%%%%% Define defaults %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Initialize optional output variables:
erp = []; amps = []; cohers = []; cohsig = []; ampsig = []; 
allamps = []; phaseangles = []; phsamp = []; sortidx = [];
auxvar = []; erpsig = []; winloc = [];

YES = 1;  % logical variables
NO  = 0;

DEFAULT_BASELINE_END = 0; % ms
TIMEX = 1;          % 1 -> plot time on x-axis; 
                    % 0 -> plot trials on x-axis

BACKCOLOR = [0.8 0.8 0.8]; % grey background
try, icadefs; catch, end;
                    % read BACKCOLOR for plot from defs file (edit this)
                    % read DEFAULT_SRATE for coher,phase,allamps, etc.

% Fix plotting text and line style parameters
SORTWIDTH = 2.5;    % Linewidth of plotted sortvar
ZEROWIDTH = 3.0;    % Linewidth of vertical 0 line
VERTWIDTH = 2.5;    % Linewidth of optional vertical lines
HORZWIDTH = 2.1;    % Linewidth of optional vertical lines
SIGNIFWIDTH = 1.9;  % Linewidth of red significance lines for amp, coher
DOTSTYLE   = 'k--'; % line style to use for vetical dotted/dashed lines
LINESTYLE = '-';    % solid line
LABELFONT = 14;     % font sizes for axis labels, tick labels
TICKFONT  = 11;

PLOT_HEIGHT = 0.2;  % fraction of y dim taken up by each time series axes
YGAP = 0.03;        % fraction gap between time axes
YEXPAND = 1.3;      % expansion factor for y-axis about erp, amp data limits

DEFAULT_AVEWIDTH  = 1; % smooth trials with this window size by default
DEFAULT_DECFACTOR = 1; % decimate by this factor by default
DEFAULT_CYCLES    = 3; % use this many cycles in amp,coher computation window
DEFAULT_CBAR      = NO;% do not plot color bar by default
DEFAULT_PHARGS = [0 25 8 13]; % Default arguments for phase sorting
DEFAULT_ALPHA     = 0.01;
alpha     = 0;      % default alpha level for coherence significance

MIN_ERPALPHA = 0.001; % significance bounds for ERP 
MAX_ERPALPHA = 0.1; 

Noshow    = NO;     % show sortvar by default
Nosort    = NO;     % sort on sortvar by default
Caxflag   = NO;     % use default caxis by default
Caxis     = [];
caxfraction = [];
Coherflag = NO;     % don't compute or show amp,coher by default
Cohsigflag= NO;     % default: do not compute coherence significance
Allampsflag=NO;     % don't image the amplitudes by default
Allcohersflag=NO;   % don't image the coherence amplitudes by default
Topoflag  = NO;     % don't plot a topoplot in upper left
Specflag  = NO;     % don't plot a spectrum in upper right
Erpflag   = NO;     % don't show erp average by default
Erpstdflag= NO; 
Erpalphaflag= NO; 
Alignflag = NO;     % don't align data to sortvar by default
Colorbar  = NO;     % if YES, plot a colorbar to right of erp image
Limitflag = NO;     % plot whole times range by default
Phaseflag = NO;     % don't sort by phase
Ampflag   = NO;     % don't sort by amplitude
Valflag   = NO;     % don't sort by value
Srateflag = NO;     % srate not given
Vertflag  = NO;
Horzflag  = NO;
Noshowflag  = NO;
Renormflag = NO;
Showwin = NO;
% yerplabel = '\muV';
yerplabel = 'ERP';
yerplabelflag = NO;
verttimes = [];
horzepochs = [];
NoTimeflag= NO;     % by default DO print "Time (ms)" below bottom axis
Signifflag= NO;     % compute significance instead of receiving it
Auxvarflag= NO;
Cycleflag = NO;
signifs   = NaN;
coherfreq = nan;    % amp/coher-calculating frequency
freq = 0;           % phase-sorting frequency
srate = DEFAULT_SRATE; % from icadefs.m
aligntime = nan;
timelimits= nan;
topomap   = [];     % topo map vector
lospecHz  = [];     % spec lo frequency
topphase = 180;     % default top phase for 'phase' option
renorm    = 'no';
noshow    = 'no';
Rmerp     = 'no';

minerp = NaN; % default limits
maxerp = NaN;
minamp = NaN;
maxamp = NaN;
mincoh = NaN;
maxcoh = NaN;
baseamp =NaN;
allamps = []; % default return 

ax1    = NaN; % default axes handles
axcb   = NaN;
ax2    = NaN;
ax3    = NaN;
ax4    = NaN;

%
%%%%%%%%%%%%%%%%%%% Test, fill in commandline args %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if nargin < 1
  help erpimage
  return
end

if nargin < 3 | isempty(times)
  if size(data,1)==1 | size(data,2)==1
   fprintf('erpimage(): either input a times vector or make data size = (frames,trials).\n')
   return
  end
  times = 1:size(data,1);
  NoTimesPassed= 1;
end

if nargin < 2 | isempty(sortvar)
  sortvar = 1:size(data,2);
  Noshow = 1; % don't plot the dummy sortvar
end

framestot = size(data,1)*size(data,2);
ntrials = length(sortvar);
if ntrials < 2
  help erpimage
  fprintf('\nerpimage(): too few trials.\n');
  return
end

frames = floor(framestot/ntrials);
if frames*ntrials ~= framestot
  help erpimage
  fprintf(...
    '\nerpimage(); length of sortvar doesn''t divide number of data elements??\n')
  return
end

if nargin < 6
  decfactor = 0;
end
if nargin < 5
  avewidth = 0;
end
if nargin<4
  titl = ''; % default no title
end
if nargin<3
  times = NO;
end
if length(times) == 1 | times == NO,  % make default times
   times = 0:frames-1;
   srate = 1000*(length(times)-1)/(times(length(times))-times(1));
   fprintf('Using sampling rate %g Hz.\n',srate);
elseif length(times) == 3
   mintime = times(1);
   frames = times(2);
   srate = times(3);
   times = mintime:1000/srate:mintime+(frames-1)*1000/srate;
   fprintf('Using sampling rate %g Hz.\n',srate);
else
   % NB: might use default srate read from icadefs here...
   srate = 1000*(length(times)-1)/(times(end)-times(1));
end
if length(times) ~= frames
   fprintf(...
'erpimage(): length(data)(%d) ~= length(sortvar)(%d) * length(times)(%d).\n\n',...
                  framestot,              length(sortvar),   length(times));
   return
end
if avewidth == 0,
  avewidth = DEFAULT_AVEWIDTH;
end
if decfactor == 0,
  decfactor = DEFAULT_DECFACTOR;
end
if avewidth < 1
  help erpimage
  fprintf('\nerpimage(): Variable avewidth cannot be < 1.\n')
  return
end
if avewidth > ntrials
  fprintf('Setting variable avewidth to max %d.\n',ntrials)
  avewidth = ntrials;  
end
if decfactor > ntrials
  fprintf('Setting variable decfactor to max %d.\n',ntrials)
  decfactor = ntrials;  
end
%
%%%%%%%%%%%%%%%%% Collect optional args %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if nargin > 6
  flagargs = [];

  for a=7:nargin % for each remaining Arg

	  Arg = eval(['arg' int2str(a-6)]);
	  if Caxflag == YES
		  if size(Arg,1) ~= 1 | size(Arg,2) > 2
			  help erpimage
			  fprintf('\nerpimage(): caxis arg must be a scalar or (1,2) vector.\n');
			  return
		  end
		  if size(Arg,2) == 2
			  Caxis = Arg;
		  else
			  caxfraction = Arg;
		  end
		  Caxflag = NO;
	  elseif Coherflag == YES
		  if length(Arg) > 3 | length(Arg) < 1
			  help erpimage
			  fprintf('\nerpimage(): coher arg must be size <= 3.\n');
			  return
		  end
		  coherfreq = Arg(1);
		  if size(Arg,2) == 1
			  coherfreq = Arg(1);
		  else
			  coherfreq = Arg(1:2);
		  end;
		  if size(Arg,2) == 3
			  Cohsigflag = YES;
			  alpha  = Arg(3);
			  if alpha < 0 | alpha > 0.1
				  fprintf('erpimage(): alpha value %g out of bounds.\n',alpha); 
				  return
			  end
		  end
		  Coherflag = NO;
		  Erpflag = YES;  % plot amp, coher below erp time series
	  elseif Topoflag == YES;
		  if length(Arg) ~= 2
			  help erpimage
			  fprintf('\nerpimage(): topo arg must be a list of length 2.\n');
			  return
		  end
		  topomap = Arg{1};
		  eloc_file = Arg{2};
		  Topoflag = NO;
	  elseif Specflag == YES;
		  if length(Arg) ~= 2
			  help erpimage
			  fprintf('\nerpimage(): spec arg must be a list of length 2.\n');
			  return
		  end
		  lospecHz = Arg(1);
		  hispecHz = Arg(2);
		  Specflag = NO;
	  elseif Renormflag == YES
		  renorm = Arg;
		  Renormflag = NO;
	  elseif Noshowflag == YES
		  noshow = Arg;
		  Noshowflag = NO;
	  elseif Alignflag == YES
		  aligntime = Arg;
		  Alignflag = NO;
	  elseif Limitflag == YES
		  %  [lotime hitime loerp hierp loamp hiamp locoher hicoher]
		  if size(Arg,1) ~= 1 | size(Arg,2) < 2 ...
				  | size(Arg,2) > 9 ...
				  help erpimage
			  fprintf('\nerpimage(): limits arg must be a vector sized (1,2<->9).\n');
			  return
		  end
		  if  ~isnan(Arg(1)) & (Arg(2) <= Arg(1))
			  help erpimage
			  fprintf('\nerpimage(): time limits out of order or out of range.\n');
			  return
		  end
		  if Arg(1) < min(times)
			  Arg(1) = min(times);
			  fprintf('Adjusting mintime limit to first data value %g\n',min(times));
		  end
		  if Arg(2) > max(times)
			  Arg(2) = max(times);
			  fprintf('Adjusting maxtime limit to last data value %g\n',max(times));
		  end
		  timelimits = Arg(1:2);
		  if length(Arg)> 2
			  minerp = Arg(3);
		  end
		  if length(Arg)> 3
			  maxerp = Arg(4);
		  end
		  if ~isnan(maxerp) & maxerp <= minerp
			  help erpimage
			  fprintf('\nerpimage(): erp limits args out of order.\n');
			  return
		  end
		  if length(Arg)> 4
			  minamp = Arg(5);
		  end
		  if length(Arg)> 5
			  maxamp = Arg(6);
		  end
		  if maxamp <= minamp
			  help erpimage
			  fprintf('\nerpimage(): amp limits args out of order.\n');
			  return
		  end
		  if length(Arg)> 6
			  mincoh = Arg(7);
		  end
		  if length(Arg)> 7
			  maxcoh = Arg(8);
		  end
		  if maxcoh <= mincoh
			  help erpimage
			  fprintf('\nerpimage(): coh limits args out of order.\n');
			  return
		  end
		  if length(Arg)>8
			  baseamp = Arg(9);    % for 'allamps'
		  end
		  Limitflag = NO;
		  
	  elseif Srateflag == YES
          srate = Arg(1);
          Srateflag = NO;
	  elseif Cycleflag == YES
          DEFAULT_CYCLES = Arg;
          Cycleflag = NO;
	  elseif Auxvarflag == YES;
          if isa(Arg,'cell')==YES & length(Arg)==2
			  auxvar = Arg{1};
			  auxcolors = Arg{2};
          elseif isa(Arg,'cell')==YES
			  fprintf('erpimage(): auxvars argument must be a matrix or length-2 cell array.\n');
			  return
          else
			  auxvar = Arg; % no auxcolors specified
          end
          Auxvarflag = NO;
	  elseif Vertflag == YES
            verttimes = Arg;
            Vertflag = NO;
	  elseif Horzflag == YES
            horzepochs = Arg;
            Horzflag = NO;
	  elseif yerplabelflag == YES
            yerplabel = Arg;
            yerplabelflag = NO;
	  elseif Signifflag == YES
		  signifs = Arg; % [low_amp hi_amp coher]