functions/sigprocfunc/erpimage.m

% ERPIMAGE - Plot a colored image of a collection of single-trial data epochs, optionally
%              sorted on and/or aligned to an input sorting variable and smoothed across
%              trials with a Gaussian weighted moving-average. (To return event-aligned data
%              without plotting, use EEGALIGN).  Optionally sort trials on value, amplitude
%              or phase within a specified latency window. Optionally plot the ERP mean
%              and std. dev.and moving-window spectral amplitude and inter-trial coherence
%              at aselected or peak frequency. Optionally 'time warp' the single trial
%              time-domain (potential) or power data to align the plotted data to a series
%              of events with varying latencies that occur in each trial. Click on
%              individual figures parts to examine them separately and zoom (using AXCOPY).
% Usage:
%            >> figure; erpimage(data,[],times); % image trials as colored lines in input order
%
%            >> figure; [outdata,outvar,outtrials,limits,axhndls, ...
%                        erp,amps,cohers,cohsig,ampsig,outamps,...
%                        phsangls,phsamp,sortidx,erpsig] ...
%                            = erpimage(data,sortvar,times,'title',avewidth,decimate,...
%                                  'key', 'val', ...); % use options
% Required input:
%   data     = [vector or matrix] Single-channel input data to image.
%               Formats (1,frames*trials) or (frames,trials)
%
% Optional ordered inputs {with defaults}:
%
%   sortvar  = [vector | []] Variable to sort epochs on (length(sortvar) = nepochs)
%              Example: sortvar may by subject response time in each epoch (in ms)
%              {default|[]: plot in input order}
%   times    = [vector | []] vector of latencies (in ms) for each epoch time point.
%               Else [startms ntimes srate] = [start latency (ms), time points
%               (=frames) per epoch, sampling rate (Hz)]. Else [] -> 0:nframes-1
%               {default: []}
%  'title'   = ['string'] Plot title {default: none}
%   avewidth = [positive scalar (may be non-integer)]. If avg_type is set to 'boxcar'
%               (the default), this is the number of trials used to smooth
%               (vertically) with a moving-average. If avg_type is set to
%               'Gaussian,' this is the standard deviation (in units of
%               trials) of the Gaussian window used to smooth (vertically)
%               with a moving-average.  Gaussian window extends three
%               standard deviations below and three standard deviations above window
%               center (trials beyond window are not incorporated into average). {default: no
%               smoothing}
%   decimate = Factor to decimate|interpolate ntrials by (may be non-integer)
%               Else, if this is large (> sqrt(ntrials)), output this many epochs.
%               {default|0->1}
%
% Optional unordered 'keyword',argument pairs:
%
% Re-align data epochs:
%   'align'  = [latency] Time-lock data to sortvar. Plot sortvar at given latency
%               (in ms). Else Inf -> plot sortvar at median sortvar latency
%               {default: do not align}
%   'timewarp' = {[events], [warpms], {colors}} Time warp ERP, amplitude and phase
%               time-courses before smoothing. 'events' is a matrix whose columns
%               specify the latencies (in ms) at which a series of successive events occur
%               in each trial. 'warpms' is an optional vector of latencies (in ms) to which
%               the series of events should be time locked. (Note: Epoch start and end
%               should not be declared as events or warpms}. If 'warpms' is absent or [],
%               the median of each 'events' column will be used. {colors} contains a
%               list of Matlab linestyles to use for vertical lines marking the occurrence
%               of the time warped events. If '', no line will be drawn for this event
%               column. If fewer colors than event columns, cycles through the given color
%               labels.  Note: Not compatible with 'vert' (below).
%   'renorm' = ['yes'|'no'| formula] Normalize sorting variable to epoch
%               latency range and plot. 'yes'= autoscale. formula must be a linear
%               transformation in the format 'a*x+b'
%               Example of formula: '3*x+2'. {default: 'no'}
%               If sorting by string values like event type, suggested formulas for:
%                 letter string: '1000*x', number string: '30000*x-1500'
%   'noplot' = ['on'|'off'] Do not plot sortvar {default: Plot sortvar if in times range}
%   'NoShow' = ['on'|'off'] Do not plot erpimage, simply return outputs {default: 'off'}
%
% Sort data epochs:
% 'nosort'       = ['on'|'off'] Do not sort data epochs. {default: Sort data epochs by
%                  sortvar (see sortvar input above)}
% 'replace_ties' = ['yes'|'no'] Replace trials with the same value of
%                  sortvar with the mean of those trials.  Only works if sorting trials
%                  by sortvar. {default: 'no'}
% 'valsort'      = [startms endms direction] Sort data on (mean) value
%                  between startms and (optional) endms. Direction is 1 or -1.
%                  If -1, plot max-value epoch at bottom {default: sort on sortvar}
% 'phasesort'    = [ms_center prct freq maxfreq topphase] Sort epochs by phase in
%                  a 3-cycle window centered at latency ms_center (ms).
%                  Percentile (prct) in range [0,100] gives percent of trials
%                  to reject for (too) low amplitude. Else, if in range [-100,0],
%                  percent of trials to reject for (too) high amplitude;
%                  freq (Hz) is the phase-sorting frequency. With optional
%                  maxfreq, sort by phase at freq of max power in the data in
%                  the range [freq,maxfreq] (Note: 'phasesort' arg freq overrides
%                  the frequency specified in 'coher'). With optional topphase,
%                  sort by phase, putting topphase (degrees, in range [-180,180])
%                  at the top of the image. Note: 'phasesort' now uses circular
%                  smoothing. Use 'cycles' (below) for wavelet length.
%                  {default: [0 25 8 13 180]}
%  'ampsort'     = [center_ms prcnt freq maxfreq]  Sort epochs by amplitude.
%                  (See 'phasesort' above). If ms_center is 'Inf', then sorting
%                  is by mean power across the time window specified by 'sortwin'
%                  below. If third arg, freq, is < 0, sort by mean power in the range
%                  [ abs(freq)   maxfreq ].
%  'sortwin'     = [start_ms end_ms] If center_ms == Inf in 'ampsort' arg (above), sorts
%                  by mean amplitude across window centers shifted from start_ms
%                  to end_ms by 10 ms.
%  'showwin'     = ['on'|'off'] Show sorting window behind ERP trace. {default: 'off'}
%
% Plot time-varying spectral amplitude instead of potential:
% 'plotamps'     = ['on'|'off'] Image amplitudes at each trial and latency instead of 
%                  potential values. Note: Currently requires 'coher' (below) with alpha signif.
%                  Use 'cycles' (see below) > (its default) 3 for better frequency specificity,
%                  {default: plot potential, not amplitudes, with no minimum}. The average power
%                  (in log space) before time 0 is automatically removed. Note that the 
%                  'baseline' parameter has no effect on 'plotamps'. Instead use
%                  change "baselinedb" or "basedB" in the 'limits' parameter. By default
%                  the baseline is removed before time 0.
%
% Specify plot parameters:
%   'limits'         = [lotime hitime minerp maxerp lodB hidB locoher hicoher basedB]
%                      Plot axes limits. Can use NaN (or nan, but not Nan) for missing items
%                      and omit late items. Use last input, basedB, to set the
%                      baseline dB amplitude in 'plotamps' plots {default: from data}
%   'sortvar_limits' = [min max] minimum and maximum sorting variable
%                      values to image. This only affects visualization of
%                      ERPimage and ERPs (not smoothing).  Cannot be used
%                      if sorting by any factor besides sortvar (e.g.,
%                      phase).
%   'signif'         = [lo_dB, hi_dB, coher_signif_level] Use precomputed significance
%                      thresholds (as from outputs ampsig, cohsig) to save time. {default: none}
%   'caxis'          = [lo hi] Set color axis limits. Else [fraction] Set caxis limits at
%                      (+/-)fraction*max(abs(data)) {default: symmetrical in dB, based on data limits}
%
% Add epoch-mean ERP to plot:
%   'erp'      = ['on'|'off'|1|2|3|4] Plot ERP time average of the trials below the
%                image.  If 'on' or 1, a single ERP (the mean of all trials) is shown.  If 2,
%                two ERPs (super and sub median trials) are shown.  If 3, the trials are split into
%                tertiles and their three ERPs are shown.  If 4, the trials are split into quartiles
%                and four ERPs are shown. Note, if you want negative voltage plotted up, change YDIR
%                to -1 in icadefs.m.  If 'erpalpha' option is used, any values of 'erp' greater than
%                1 will be reset to 1. {default no ERP plotted}
%   'erpalpha' = [alpha] Visualizes two-sided significance threshold (i.e., a two-tailed test) for the
%                null hypothesis of a zero mean, symmetric distribution (range: [.001 0.1]). Thresholds
%                are determined via a permutation test. Requires 'erp' to be a value other than 'off'.
%                If 'erp' is set  to a value greater than 1, it is reset to 1 to increase plot readability.
%                {default: no alpha significance thresholds plotted}
%   'erpstd'   = ['on'|'off'] Plot ERP +/- stdev. Requires 'erp' {default: no std. dev. plotted}
%   'erp_grid' = If 'erp_grid' is added as an option voltage axis dashed grid lines will be
%                 added to the ERP plot to facilitate judging ERP amplitude
%   'rmerp'    = ['on'|'off'] Subtract the average ERP from each trial before processing {default: no}
%
% Add time/frequency information:
%   'coher'  = [freq] Plot ERP average plus mean amplitude & coherence at freq (Hz)
%               Else [minfrq maxfrq] = same, but select frequency with max power in
%               given range. (Note: the 'phasesort' freq (above) overwrites these
%               parameters). Else [minfrq maxfrq alpha] = plot coher. signif. level line
%               at probability alpha (range: [0,0.1]) {default: no coher, no alpha level}
%   'srate'  = [freq] Specify the data sampling rate in Hz for amp/coher (if not
%               implicit in third arg, times) {default: as defined in icadefs.m}
%   'cycles' = [float] Number of cycles in the wavelet time/frequency decomposition {default: 3}
%
% Add plot features:
%   'cbar'           = ['on'|'off'] Plot color bar to right of ERP-image {default no}
%   'cbar_title'     = [string] The title for the color bar (e.g., '\muV' for
%                       microvolts).
%   'topo'           = {map,chan_locs,eloc_info} Plot a 2-D scalp map at upper left of image.
%                       map may be a single integer, representing the plotted data channel,
%                       or a vector of scalp map channel values. chan_locs may be a channel locations
%                       file or a chanlocs structure (EEG.chanlocs). See '>> topoplot example'
%                       eloc_info (EEG.chaninfo), if empty ([]) or absent, implies the 'X' direction
%                       points towards the nose and all channels are plotted {default: no scalp map}
%   'spec'           = [loHz,hiHz] Plot the mean data spectrum at upper right of image.
%   'specaxis'       = ['log'|'lin] Use 'lin' for linear spectrum frequency scaling, 
%                       else 'log' for log scaling {default: 'log'}
%   'horz'           = [epochs_vector] Plot horizontal lines at specified epoch numbers.
%   'vert'           = [times_vector] Plot vertical dashed lines at specified latencies
%   'auxvar'         = [size(nvars,ntrials) matrix] Plot auxiliary variable(s) for each trial
%                       as separate traces. Else, 'auxvar',{[matrix],{colorstrings}}
%                       to specify N trace colors.  Ex: colorstrings = {'r','bo-','','k:'}
%                       (see also: 'vert' and 'timewarp' above). {default: none}
%   'sortvarpercent' = [float vector] Plot percentiles for the sorting variable
%                       for instance, [0.1 0.5 0.9] plots the 10th percentile, the median
%                       and the 90th percentile.
% Plot options:
% 'noxlabel'          = ['on'|'off'] Do not plot "Time (ms)" on the bottom x-axis
% 'yerplabel'         = ['string'] ERP ordinate axis label (default is ERP). Print uV with '\muV'
% 'avg_type'          = ['boxcar'|'Gaussian'] The type of moving average used to smooth
%                        the data. 'Boxcar' smoothes the data by simply taking the mean of
%                        a certain number of trials above and below each trial.
%                        'Gaussian' does the same but first weights the trials
%                        according to a Gaussian distribution (e.g., nearby trials
%                        receive greater weight).  The Gaussian is better than the
%                        boxcar in that it rather evenly filters out high frequency
%                        vertical components in the ERPimage. See 'avewidth' argument
%                        description for more information. {default: boxcar}
% 'img_trialax_label' = ['string'] The label of the axis corresponding to trials in the ERPimage
%                        (e.g., 'Reaction Time').  Note, if img_trialax_label is set to something
%                        besides 'Trials' or [], the tick marks on this axis will be set in units
%                        of the sorting variable.  This is a useful alternative to plotting the
%                        sorting variable when the sorting variable is not in milliseconds. This
%                        option is not effective if sorting by amplitude, phase, or EEG value. {default: 'Trials'}
% 'img_trialax_ticks' = Vector of sorting variable values at which tick marks (e.g., [300 350 400 450]
%                        for reaction time in msec) will appear on the trial axis of the erpimage. Tick mark
%                        values should be given in units img_trialax_label (e.g., 'Trials' or msec).
%                        This option is not effective if sorting by amplitude, phase, or EEG value.
%                        {default: automatic}
% 'baseline'          = [low_boundary high_boundary] a time window (in msec) whose mean amplitude in
%                        each trial will be removed from each trial (e.g., [-100 0]) after filtering.
%                        Useful in conjunction with 'filt' option to re-basline trials after they have been
%                        filtered. Not necessary if data have already been baselined and erpimage
%                        processing does not affect baseline amplitude {default: no further baselining
%                        of data}.
% 'baselinedb'        = [low_boundary high_boundary] a time window (in msec) whose mean amplitude in
%                        each trial will be removed from each trial (e.g., [-100 0]). Use basedB in limits
%                        to remove a fixed value. Default is time before 0. If you do not want to use a 
%                        baseline for amplitude plotting, enter a NaN value.
% 'filt'              = [low_boundary high_boundary] a two element vector indicating the frequency
%                        cut-offs for a 3rd order Butterworth filter that will be applied to each
%                        trial of data.  If low_boundary=0, the filter is a low pass filter.  If
%                        high_boundary=srate/2, then the filter is a high pass filter.  If both
%                        boundaries are between 0 and srate/2, then the filter is a bandpass filter.
%                        If both boundaries are between 0 and -srate/2, then the filter is a bandstop
%                        filter (with boundaries equal to the absolute values of low_boundary and
%                        high_boundary).  Note, using this option requires the 'srate' option to be
%                        specified and the signal processing toolbox function butter.m.  You should
%                        probably use the 'baseline' option as well since the mean prestimulus baseline
%                        may no longer be 0 after the filter is applied {default: no filtering}
%
% Optional outputs:
%    outdata  = (times,epochsout) data matrix (after smoothing)
%     outvar  = (1,epochsout) actual values trials are sorted on (after smoothing).
%               if 'sortvarpercent' is used, this variable contains a cell array with
%               { sorted_values { sorted_percent1 ... sorted_percentN } }
%   outtrials = (1,epochsout)  smoothed trial numbers
%     limits  = (1,10) array, 1-9 as in 'limits' above, then analysis frequency (Hz)
%    axhndls  = vector of 1-7 plot axes handles (img,cbar,erp,amp,coh,topo,spec)
%        erp  = plotted ERP average
%       amps  = mean amplitude time course
%      coher  = mean inter-trial phase coherence time course
%     cohsig  = coherence significance level
%     ampsig  = amplitude significance levels [lo high]
%    outamps  = matrix of imaged amplitudes (from option 'plotamps')
%   phsangls  = vector of sorted trial phases at the phase-sorting frequency
%     phsamp  = vector of sorted trial amplitudes at the phase-sorting frequency
%    sortidx  = indices of input data epochs in the sorting order
%     erpsig  = trial average significance levels [2,frames]
%
% Example:  >> figure;
%              erpimage(data,RTs,[-400 256 256],'Test',1,1,...
%                            'erp','cbar','vert',-350);
%    Plots an ERP-image of 1-s data epochs sampled at 256 Hz, sorted by RTs, with
%    title ('Test'), and sorted epochs not smoothed or decimated (1,1). Overplots
%    the (unsmoothed) RT latencies on the colored ERP-image. Also plots the
%    epoch-mean (ERP), a color bar, and a dashed vertical line at -350 ms.
%
% Authors: Scott Makeig, Tzyy-Ping Jung & Arnaud Delorme,
%          CNL/Salk Institute, La Jolla, 3-2-1998 -
%
% See also: phasecoher, rmbase, cbar, movav

% Copyright (C) Scott Makeig & Tzyy-Ping Jung, CNL / Salk Institute, La Jolla 3-2-98
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% Uses external toolbox functions: PHASECOHER, RMBASE, CBAR, MOVAV
% Uses included functions:         PLOT1TRACE, PHASEDET

% UNIMPLEMENTED - 'allcohers',[data2] -> image the coherences at each latency & epoch.
%                 Requires arg 'coher' with alpha significance.
%                 Shows projection on grand mean coherence vector at each latency
%                 and trial. {default: no}

%% LOG COMMENTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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,arg27,arg28,arg29,arg30,arg31,arg32,arg33,arg34,arg35,arg36,arg37,arg38,arg39,arg40,arg41,arg42,arg43,arg44,arg45,arg46,arg47,arg48,arg49,arg50,arg51,arg52,arg53,arg54,arg55,arg56,arg57,arg58,arg59,arg60,arg61,arg62,arg63,arg64,arg65,arg66)

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

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.
% read YDIR for plotting ERP
% 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 vertical dotted/dashed lines
LINESTYLE = '-';    % solid line
LABELFONT = 10;     % font sizes for axis labels, tick labels
TICKFONT  = 10;

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_SDEV  = 1/7; % smooth trials with this window size by default if Gaussian window
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
cycles = DEFAULT_CYCLES;
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;

NoShowVar = NO;     % show sortvar by default
Nosort    = NO;     % sort on sortvar by default
Caxflag   = NO;     % use default caxis by default

timestretchflag = NO; % Added -JH

mvavg_type='boxcar'; % use the original rectangular moving average -DG
erp_grid = NO; % add y-tick grids to ERP plot -DG
cbar_title = []; % title to add above ERPimage color bar (e.g., '\muV') -DG
img_ylab = 'Trials'; % make the ERPimage y-axis in units of the sorting variable -DG
img_ytick_lab = []; % the values at which tick marks will appear on the trial axis of the ERPimage (the y-axis by default).
%Note, this is in units of the sorting variable if img_ylab~='Trials', otherwise it is in units of trials -DG
baseline   = []; %time window of each trial whose mean amp will be used to baseline the trial -DG
baselinedb = []; %time window of each trial whose mean power will be used to baseline the trial
flt=[]; %frequency domain filter parameters -DG
sortvar_limits=[]; %plotting limits for sorting variable/trials; limits only affect visualization, not smoothing -DG
replace_ties = NO; %if YES, trials with the exact same value of a sorting variable will be replaced by their average -DG
erp_vltg_ticks=[]; %if non-empty, these are the voltage axis ticks for the ERP plot

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
SpecAxisflag = NO;  % don't change the default spectrum axis type from default
SpecAxis = 'log';   % default log frequency spectrum axis (if Specflag)
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
Sortwinflag = NO;   % sort by amplitude over a window
Valflag   = NO;     % don't sort by value
Srateflag = NO;     % srate not given
Vertflag  = NO;
Horzflag  = NO;
titleflag = NO;
NoShowflag  = NO;
Renormflag = NO;
Showwin = NO;
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;
plotmodeflag= NO;
plotmode = 'normal';
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';
percentiles = [];
percentileflag = NO;
erp_ptiles = 1;

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;


timeStretchRef = [];
timeStretchMarks = [];
tsurdata = []; % time-stretched data before smoothing (for timestretched-erp computation)
% If time-stretching is off, this variable remains empty
%
%% %%%%%%%%%%%%%%%%% Test, fill in commandline args %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if nargin < 1
    help erpimage
    return
end

data = squeeze(data);
if nargin < 3 || isempty(times)
    if size(data,1)==1 || size(data,2)==1
        fprintf('\nerpimage(): 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);
    NoShowVar = 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 = DEFAULT_AVEWIDTH;
end
if nargin<4
    titl = ''; % default no title
end
if nargin<3
    times = NO;
end
if (length(times) == 1) || (length(times) == 1 && times(1) == 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
    % Note: might use default srate read from icadefs here...
    srate = 1000*(length(times)-1)/(times(end)-times(1));
end
if length(times) ~= frames
    fprintf(...
        '\nerpimage(): length(data)(%d) ~= length(sortvar)(%d) * length(times)(%d).\n\n',...
        framestot,              length(sortvar),   length(times));
    return
end

if decfactor == 0
    decfactor = DEFAULT_DECFACTOR;
elseif decfactor > ntrials
    fprintf('Setting variable decfactor to max %d.\n',ntrials)
    decfactor = ntrials;
end
%
%% %%%%%%%%%%%%%%% Collect optional args %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if nargin > 6
    flagargs = [];
    
    a = 6;
    while a < nargin % for each remaining Arg
        a = a + 1;
        
        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 timestretchflag == YES % Added -JH
            timeStretchMarks = Arg{1};
            timeStretchMarks = round(1+(timeStretchMarks-times(1))*srate/1000); % convert from ms to frames -sm
            [smc, smr] = find(diff(timeStretchMarks') < 0);
            if ~isempty(smr)
                fprintf('\nerpimage(): Timewarp event latencies not in ascending order in trial %d.\n',smr)
                return
            end
            
            timeStretchMarks = [ ...
                repmat(1, [size(timeStretchMarks,1), 1]), ...% Epoch begins
                timeStretchMarks, ...
                repmat(length(times), [size(timeStretchMarks,1), 1])]; % Epoch ends
            if length(Arg) < 2 || isempty(Arg{2})
                timeStretchRef = median(timeStretchMarks);
            else
                timeStretchRef = Arg{2};
                timeStretchRef = round(1+(timeStretchRef-times(1))*srate/1000); % convert from ms to frames -sm
                timeStretchRef = [1 timeStretchRef length(times)];  % add epoch beginning, end
            end
            if length(Arg) < 3 || isempty(Arg{3})
                timeStretchColors = {};
            else
                timeStretchColors = Arg{3};
            end
            fprintf('The %d events specified in each trial will be time warped to latencies:',length(timeStretchRef)-2);
            fprintf(' %.0f', times(1)+1000*(timeStretchRef(2:end-1)-1)/srate); % converted from frames to ms -sm
            fprintf(' ms\n');
            timestretchflag = 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('\nerpimage(): 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 or 3.\n');
                return
            end
            topomap = Arg{1};
            eloc_file = Arg{2};
            if length(Arg) > 2, eloc_info = Arg{3};
            else                eloc_info = [];
            end
            Topoflag = NO;
        elseif Specflag == YES
            if length(Arg) ~= 2 
                error('\nerpimage(): ''spec'' flag argument must be a numeric array of length 2.\n');
                return
            end
            lospecHz = Arg(1);
            hispecHz = Arg(2);
            Specflag = NO;
        elseif SpecAxisflag == YES
            SpecAxis = Arg;
            if ~strcmpi(SpecAxis,'lin') && ~strcmpi(SpecAxis,'log')
              error('\nerpimage(): spectrum axis type must be ''lin'' or ''log''.\n');
              return
            end
            if strcmpi(SpecAxis,'lin')
               SpecAxis = 'linear';   % convert to MATLAB Xscale keyword
            end
            SpecAxisflag = NO;
        elseif Renormflag == YES
            renorm = Arg;
            Renormflag = NO;
        elseif NoShowflag == YES
            NoShow = Arg;
            if strcmpi(NoShow, 'off'), NoShow = 'no'; end
            NoShowflag = NO;
        elseif Alignflag == YES
            aligntime = Arg;
            Alignflag = NO;
        elseif percentileflag == YES
            percentiles = Arg;
            percentileflag = 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
            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('\nerpimage(): auxvars argument must be a matrix or length-2 cell array.\n');
                return
            else
                auxvar = Arg; % no auxcolors specified
            end
            [xr,xc] = size(auxvar);
            lns = length(sortvar);
            if xr ~= lns && xc ~= lns
                error('\nerpimage(): auxvar columns different from the number of epochs in data');
            elseif xr == lns && xc ~= lns
                auxvar = auxvar';   % exchange rows/cols
            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]
            if length(signifs) ~= 3
                fprintf('\nerpimage(): signif arg [%g] must have 3 values\n',Arg);
                return
            end
            Signifflag = NO;
        elseif Allcohersflag == YES
            data2=Arg;
            if size(data2) ~= size(data)
                fprintf('\nerpimage(): allcohers data matrix must be the same size as data.\n');
                return
            end
            Allcohersflag = NO;
        elseif Phaseflag == YES
            n = length(Arg);
            if n > 5
                error('\nerpimage(): Too many arguments for keyword ''phasesort''');
            end
            phargs = Arg;
            
            if phargs(3) < 0
                error('\nerpimage(): Invalid negative frequency argument for keyword ''phasesort''');
            end
            if n>=4
                if phargs(4) < 0
                    error('\nerpimage(): Invalid negative argument for keyword ''phasesort''');
                end
            end
            if min(phargs(1)) < times(1) || max(phargs(1)) > times(end)
                error('\nerpimage(): time for phase sorting filter out of bound.');
            end
            if phargs(2) >= 100 || phargs(2) < -100
                error('\nerpimage(): %-argument for keyword ''phasesort'' must be (-100;100)');
            end
            if length(phargs) >= 4 && phargs(3) > phargs(4)
                error('\nerpimage(): Phase sorting frequency range must be increasing.');
            end
            if length(phargs) == 5
                topphase = phargs(5);
            end
            Phaseflag = NO;
        elseif Sortwinflag == YES % 'ampsort' mean amplitude over a time window
            n = length(Arg);
            sortwinarg = Arg;
            if n > 2
                error('\nerpimage(): Too many arguments for keyword ''sortwin''');
            end
            if min(sortwinarg(1)) < times(1) || max(sortwinarg(1)) > times(end)
                error('\nerpimage(): start time for value sorting out of bounds.');
            end
            if n > 1
                if min(sortwinarg(2)) < times(1) || max(sortwinarg(2)) > times(end)
                    error('\nerpimage(): end time for value sorting out of bounds.');
                end
            end
            if n > 1 && sortwinarg(1) > sortwinarg(2)
                error('\nerpimage(): Value sorting time range must be increasing.');
            end
            Sortwinflag = NO;
        elseif Ampflag == YES % 'ampsort',[center_time,prcnt_reject,minfreq,maxfreq]
            n = length(Arg);
            if n > 4
                error('\nerpimage(): Too many arguments for keyword ''ampsort''');
            end
            ampargs = Arg;
            
            % if ampargs(3) < 0
            %    error('\nerpimage(): Invalid negative argument for keyword ''ampsort''');
            % end
            if n>=4
                if ampargs(4) < 0
                    error('\nerpimage(): Invalid negative argument for keyword ''ampsort''');
                end
            end
            
            if ~isinf(ampargs(1))
                if min(ampargs(1)) < times(1) || max(ampargs(1)) > times(end)
                    error('\nerpimage(): time for amplitude sorting filter out of bounds.');
                end
            end
            
            if ampargs(2) >= 100 || ampargs(2) < -100
                error('\nerpimage(): percentile argument for keyword ''ampsort'' must be (-100;100)');
            end
            
            if length(ampargs) == 4 && abs(ampargs(3)) > abs(ampargs(4))
                error('\nerpimage(): Amplitude sorting frequency range must be increasing.');
            end
            Ampflag = NO;
            
        elseif Valflag == YES % sort by potential value in a given window
            % Usage: 'valsort',[mintime,maxtime,direction]
            n = length(Arg);
            if n > 3
                error('\nerpimage(): Too many arguments for keyword ''valsort''');
            end
            valargs = Arg;
            
            if min(valargs(1)) < times(1) || max(valargs(1)) > times(end)
                error('\nerpimage(): start time for value sorting out of bounds.');
            end
            if n > 1
                if min(valargs(2)) < times(1) || max(valargs(2)) > times(end)
                    error('\nerpimage(): end time for value sorting out of bounds.');
                end
            end
            if n > 1 && valargs(1) > valargs(2)
                error('\nerpimage(): Value sorting time range must be increasing.');
            end
            if n==3 && (~isnumeric(valargs(3)) || valargs(3)==0)
                error('\nerpimage(): Value sorting direction must be +1 or -1.');
            end
            Valflag = NO;
        elseif plotmodeflag == YES
            plotmode = Arg; plotmodeflag = NO;
        elseif titleflag == YES
            titl = Arg; titleflag = NO;
        elseif Erpalphaflag == YES
            erpalpha = Arg(1);
            if erpalpha < MIN_ERPALPHA || erpalpha > MAX_ERPALPHA
                fprintf('\nerpimage(): erpalpha value is out of bounds [%g, %g]\n',...
                    MIN_ERPALPHA,MAX_ERPALPHA);
                return
            end
            Erpalphaflag = NO;
            % -----------------------------------------------------------------------
            % -----------------------------------------------------------------------
            % -----------------------------------------------------------------------
        elseif strcmpi(Arg,'avg_type')
            if a < nargin
                a=a+1;
                Arg = eval(['arg' int2str(a-6)]);
                if strcmpi(Arg, 'Gaussian'), mvavg_type='gaussian';
                elseif strcmpi(Arg, 'Boxcar'), mvavg_type='boxcar';
                else error('\nerpimage(): Invalid value for optional argument ''avg_type''.');
                end
            else
                error('\nerpimage(): Optional argument ''avg_type'' needs to be assigned a value.');
            end
        elseif strcmp(Arg,'nosort')
            Nosort = YES;
            if a < nargin
                Arg = eval(['arg' int2str(a+1-6)]);
                if strcmpi(Arg, 'on')
                    Nosort = YES;
                    a = a + 1;
                elseif strcmpi(Arg, 'off')
                    Nosort = NO;
                    a = a + 1;
                end
            end
        elseif strcmp(Arg,'showwin')
            Showwin = YES;
            if a < nargin
                Arg = eval(['arg' int2str(a+1-6)]);
                if strcmpi(Arg, 'on') 
                    Showwin = YES; a = a+1;
                elseif strcmpi(Arg, 'off') 
                    Showwin = NO;  a = a+1;
                end
            end
        elseif strcmp(Arg,'noplot') % elseif strcmp(Arg,'NoShow') % by Luca && Ramon
            NoShow = YES;
            if a < nargin
                Arg = eval(['arg' int2str(a+1-6)]);
                if strcmpi(Arg, 'on'),     NoShow = YES; a = a+1;
                elseif strcmpi(Arg, 'off'), NoShow = NO;  a = a+1;
                end
            end
        elseif strcmpi(Arg,'replace_ties')
            if a < nargin
                a = a+1;
                temp = eval(['arg' int2str(a-6)]);
                if strcmpi(temp,'on')
                    replace_ties = YES;
                elseif strcmpi(temp,'off') 
                    replace_ties = NO;
                else
                    error('\nerpimage(): Argument ''replace_ties'' needs to be followed by the string ''on'' or ''off''.');
                end
            else
                error('\nerpimage(): Argument ''replace_ties'' needs to be followed by the string ''on'' or ''off''.');
            end
        elseif strcmpi(Arg,'sortvar_limits')
            if a < nargin
                a = a+1;
                sortvar_limits = eval(['arg' int2str(a-6)]);
                if ischar(sortvar_limits) || length(sortvar_limits)~=2
                    error('\nerpimage(): Argument ''sortvar_limits'' needs to be followed by a two element vector.');
                end
            else
                error('\nerpimage(): Argument ''sortvar_limits'' needs to be followed by a two element vector.');
            end
        elseif strcmpi(Arg,'erp')
            Erpflag = YES;
            erp_ptiles=1;
            if a < nargin
                Arg = eval(['arg' int2str(a+1-6)]);
                if strcmpi(Arg, 'on'),      Erpflag = YES; erp_ptiles=1; a = a+1;
                elseif strcmpi(Arg, 'off'), Erpflag = NO;  a = a+1;
                elseif strcmpi(Arg,'1') || (Arg(1)==1), Erplag = YES; erp_ptiles=1; a=a+1;
                elseif strcmpi(Arg,'2') || (Arg(1)==2), Erplag = YES; erp_ptiles=2; a=a+1;
                elseif strcmpi(Arg,'3') || (Arg(1)==3), Erplag = YES; erp_ptiles=3; a=a+1;
                elseif strcmpi(Arg,'4') || (Arg(1)==4), Erplag = YES; erp_ptiles=4; a=a+1;
                end
            end
        elseif strcmpi(Arg,'rmerp')
            Rmerp = 'yes';
            if a < nargin
                Arg = eval(['arg' int2str(a+1-6)]);
                if strcmpi(Arg, 'on')   
                    Rmerp = 'yes'; a = a+1;
                elseif strcmpi(Arg, 'off') 
                    Rmerp = 'no';  a = a+1;
                end
            end
        elseif strcmp(Arg,'cbar') || strcmp(Arg,'colorbar')
            Colorbar = YES;
            if a < nargin
                Arg = eval(['arg' int2str(a+1-6)]);
                if strcmpi(Arg, 'on')   
                    Colorbar = YES; a = a+1;
                elseif strcmpi(Arg, 'off') 
                    Colorbar = NO;  a = a+1;
                end
            end
        elseif (strcmp(Arg,'allamps') || strcmp(Arg,'plotamps'))
            Allampsflag = YES;
            if a < nargin
                Arg = eval(['arg' int2str(a+1-6)]);
                if strcmpi(Arg, 'on')
                    Allampsflag = YES; a = a+1;
                elseif strcmpi(Arg, 'off') 
                    Allampsflag = NO;  a = a+1;
                end
            end
        elseif strcmpi(Arg,'erpstd')
            Erpstdflag = YES;
            if a < nargin
                Arg = eval(['arg' int2str(a+1-6)]);
                if strcmpi(Arg, 'on')
                    Erpstdflag = YES; a = a+1;
                elseif strcmpi(Arg, 'off') 
                    Erpstdflag = NO;  a = a+1;
                end
            end
        elseif strcmp(Arg,'noxlabel') || strcmp(Arg,'noxlabels') || strcmp(Arg,'nox')
            NoTimeflag = YES;
            if a < nargin
                Arg = eval(['arg' int2str(a+1-6)]);
                if strcmpi(Arg, 'on')
                    NoTimeflag = YES; a = a+1;
                elseif strcmpi(Arg, 'off') 
                    NoTimeflag = NO;  a = a+1;
                end
            end
        elseif strcmp(Arg,'plotmode')
            plotmodeflag = YES;
        elseif strcmp(Arg,'sortvarpercent')
            percentileflag = YES;
        elseif strcmp(Arg,'renorm')
            Renormflag = YES;
        elseif strcmp(Arg,'NoShow')
            NoShowflag = YES;
        elseif strcmp(Arg,'caxis')
            Caxflag = YES;
        elseif strcmp(Arg,'title')
            titleflag = YES;
        elseif strcmp(Arg,'coher')
            Coherflag = YES;
        elseif strcmp(Arg,'timestretch') || strcmp(Arg,'timewarp') % Added -JH
            timestretchflag = YES;
        elseif strcmp(Arg,'allcohers')
            Allcohersflag = YES;
        elseif strcmp(Arg,'topo') || strcmp(Arg,'topoplot')
            Topoflag = YES;
        elseif strcmp(Arg,'spec') || strcmp(Arg,'spectrum')
            Specflag = YES;
        elseif strcmp(Arg,'Specaxis') || strcmp(Arg,'specaxis') || strcmp(Arg,'SpecAxis')
            SpecAxisflag = YES;
        elseif strcmpi(Arg,'erpalpha')
            Erpalphaflag = YES;
        elseif strcmp(Arg,'align')
            Alignflag = YES;
        elseif strcmp(Arg,'limits')
            Limitflag = YES;
        elseif (strcmp(Arg,'phase') || strcmp(Arg,'phasesort'))
            Phaseflag = YES;
        elseif strcmp(Arg,'ampsort')
            Ampflag = YES;
        elseif strcmp(Arg,'sortwin')
            Sortwinflag = YES;
        elseif strcmp(Arg,'valsort')
            Valflag = YES;
        elseif strcmp(Arg,'auxvar')
            Auxvarflag = YES;
        elseif strcmp(Arg,'cycles')
            Cycleflag = YES;
        elseif strcmpi(Arg,'yerplabel')
            yerplabelflag = YES;
        elseif strcmpi(Arg,'srate')
            Srateflag = YES;
        elseif strcmpi(Arg,'erp_grid')
            erp_grid = YES;
        elseif strcmpi(Arg,'baseline')
            if a < nargin
                a = a+1;
                baseline = eval(['arg' int2str(a-6)]);
            else
                error('\nerpimage(): Argument ''baseline'' needs to be followed by a two element vector.');
            end
        elseif strcmpi(Arg,'baselinedb')
            if a < nargin
                a = a+1;
                baselinedb = eval(['arg' int2str(a-6)]);
                if length(baselinedb) > 2, error('''baselinedb'' need to be a 2 argument vector'); end
            else
                error('\nerpimage(): Argument ''baselinedb'' needs to be followed by a two element vector.');
            end
        elseif strcmpi(Arg,'filt')
            if a < nargin
                a = a+1;
                flt = eval(['arg' int2str(a-6)]);
            else
                error('\nerpimage(): Argument ''filt'' needs to be followed by a two element vector.');
            end
        elseif strcmpi(Arg,'erp_vltg_ticks')
            if a < nargin
                a = a+1;
                erp_vltg_ticks=eval(['arg' int2str(a-6)]);
            else
                error('\nerpimage(): Argument ''erp_vltg_ticks'' needs to be followed by a vector.');
            end
        elseif strcmpi(Arg,'img_trialax_label')
            if a < nargin
                a = a+1;
                img_ylab = eval(['arg' int2str(a-6)]);
            else
                error('\nerpimage(): Argument ''img_trialax_label'' needs to be followed by a string.');
            end
        elseif strcmpi(Arg,'img_trialax_ticks')
            if a < nargin
                a = a+1;
                img_ytick_lab = eval(['arg' int2str(a-6)]);
            else
                error('\nerpimage(): Argument ''img_trialax_ticks'' needs to be followed by a vector of values at which tick marks will appear.');
            end
        elseif strcmpi(Arg,'cbar_title')
            if a < nargin
                a = a+1;
                cbar_title = eval(['arg' int2str(a-6)]);
            else
                error('\nerpimage(): Argument ''cbar_title'' needs to be followed by a string.');
            end
        elseif strcmp(Arg,'vert') ||  strcmp(Arg,'verttimes')
            Vertflag = YES;
        elseif strcmp(Arg,'horz') ||  strcmp(Arg,'horiz') || strcmp(Arg,'horizontal')
            Horzflag = YES;
        elseif strcmp(Arg,'signif') || strcmp(Arg,'signifs') || strcmp(Arg,'sig') || strcmp(Arg,'sigs')
            Signifflag = YES;
        else
            help erpimage
            if ischar(Arg)
                fprintf('\nerpimage(): unknown arg %s\n',Arg);
            else
                fprintf('\nerpimage(): unknown arg %d, size(%d,%d)\n',a,size(Arg,1),size(Arg,2));
            end
            return
        end
    end % Arg
end

if exist('img_ylab','var') || exist('img_ytick_lab','var'),
    oops=0;
    if exist('phargs','var'),
        fprintf('********* Warning *********\n');
        fprintf('Options ''img_ylab'' and ''img_ytick_lab'' have no effect when sorting by phase.\n');
        oops=0;
    elseif exist('valargs','var'),
        fprintf('********* Warning *********\n');
        fprintf('Options ''img_ylab'' and ''img_ytick_lab'' have no effect when sorting by EEG voltage.\n');
        oops=0;
    elseif exist('ampargs','var'),
        fprintf('********* Warning *********\n');
        fprintf('Options ''img_ylab'' and ''img_ytick_lab'' have no effect when sorting by frequency amplitude.\n');
        oops=0;
    end
    if oops
        img_ylab=[];
        img_ytick_lab=[];
    end
end

if   Caxflag == YES ...
        |Coherflag == YES ...
        |Alignflag == YES ...
        |Limitflag == YES
    help erpimage
    fprintf('\nerpimage(): missing option arg.\n')
    return
end
if (Allampsflag || exist('data2')) && ( any(isnan(coherfreq)) || ~Cohsigflag )
    fprintf('\nerpimage(): allamps and allcohers flags require coher freq, srate, and cohsig.\n');
    return
end
if Allampsflag && exist('data2')
    fprintf('\nerpimage(): cannot image both allamps and allcohers.\n');
    return
end
if ~exist('srate') || srate <= 0
    fprintf('\nerpimage(): Data srate must be specified and > 0.\n');
    return
end
if ~isempty(auxvar)
    % whos auxvar
    if size(auxvar,1) ~= ntrials && size(auxvar,2) ~= ntrials
        fprintf('\nerpimage(): auxvar size should be (N,ntrials), e.g., (N,%d)\n',...
            ntrials);
        return
    end
    if size(auxvar,1) == ntrials && size(auxvar,2) ~= ntrials  % make (N,frames)
        auxvar = auxvar';
    end
    if size(auxvar,2) ~= ntrials
        fprintf('\nerpimage(): auxvar size should be (N,ntrials), e.g., (N,%d)\n',...
            ntrials);
        return
    end
    if exist('auxcolors')==YES % if specified
        if isa(auxcolors,'cell')==NO % if auxcolors is not a cell array
            fprintf(...
                '\nerpimage(): auxcolors argument to auxvar flag must be a cell array.\n');
            return
        end
    end
elseif exist('timeStretchRef') && ~isempty(timeStretchRef)
    if ~isnan(aligntime)
        fprintf(['\nerpimage(): options "align" and ' ...
            '"timewarp" are not compatible.\n']);
        return;
    end
    
    if ~isempty(timeStretchColors)
        if length(timeStretchColors) < length(timeStretchRef)
            nColors = length(timeStretchColors);
            for k=nColors+1:length(timeStretchRef)-2
                timeStretchColors = { timeStretchColors{:} timeStretchColors{1+rem(k-1,nColors)}};
            end
        end
        timeStretchColors = {'' timeStretchColors{:} ''};
    else
        timeStretchColors = { 'k--'};
        for k=2:length(timeStretchRef)
            timeStretchColors = { timeStretchColors{:} 'k--'};
        end
    end
    
    
    auxvarInd = 1-strcmp('',timeStretchColors); % indicate which lines to draw
    newauxvars = ((timeStretchRef(find(auxvarInd))-1)/srate+times(1)/1000) * 1000; % convert back to ms
    fprintf('Overwriting vert with auxvar\n');
    verttimes = [newauxvars'];
    verttimesColors = {timeStretchColors{find(auxvarInd)}};
    newauxvars = repmat(newauxvars, [1 ntrials]);
    
    if isempty(auxvar) % Initialize auxvar && auxcolors
        %  auxvar = newauxvars;
        auxcolors = {timeStretchColors{find(auxvarInd)}};
    else % Append auxvar & auxcolors
        if ~exist('auxcolors')
            % Fill with default color (k-- for now)
            auxcolors = {};
            for j=1:size(auxvar,1)
                auxcolors{end+1} = 'k--';
            end
        end
        for j=find(auxvarInd)
            auxcolors{end+1} = timeStretchColors{j};
        end
        auxvar = [auxvar; newauxvars];
    end
end

% No need to turn off ERP anymore; we now time-stretch potentials
% separately (see tsurdata below) so the (warped) ERP is actually meaningful

if exist('phargs')
    if phargs(3) > srate/2
        fprintf(...
            '\nerpimage(): Phase-sorting frequency (%g Hz) must be less than Nyquist rate (%g Hz).',...
            phargs(3),srate/2);
    end
    
    if frames < cycles*srate/phargs(3)
        fprintf('\nerpimage(): phase-sorting freq. (%g) too low: epoch length < %d cycles.\n',...
            phargs(3),cycles);
        return
    end
    if length(phargs)==4 && phargs(4) > srate/2
        phargs(4) = srate/2;
    end
    if length(phargs)==5 && (phargs(5)>180 || phargs(5) < -180)
        fprintf('\nerpimage(): coher topphase (%g) out of range.\n',topphase);
        return
    end
end
if exist('ampargs')
    if abs(ampargs(3)) > srate/2
        fprintf(...
            '\nerpimage(): amplitude-sorting frequency (%g Hz) must be less than Nyquist rate (%g Hz).',...
            abs(ampargs(3)),srate/2);
    end
    
    if frames < cycles*srate/abs(ampargs(3))
        fprintf('\nerpimage(): amplitude-sorting freq. (%g) too low: epoch length < %d cycles.\n',...
            abs(ampargs(3)),cycles);
        return
    end
    if length(ampargs)==4 && abs(ampargs(4)) > srate/2
        ampargs(4) = srate/2;
        fprintf('> Reducing max ''ampsort'' frequency to Nyquist rate (%g Hz)\n',srate/2)
    end
end
if ~any(isnan(coherfreq))
    if coherfreq(1) <= 0 || srate <= 0
        fprintf('\nerpimage(): coher frequency (%g) out of range.\n',coherfreq(1));
        return
    end
    if coherfreq(end) > srate/2 || srate <= 0
        fprintf('\nerpimage(): coher frequency (%g) out of range.\n',coherfreq(end));
        return
    end
    if frames < cycles*srate/coherfreq(1)
        fprintf('\nerpimage(): coher freq. (%g) too low:  epoch length < %d cycles.\n',...
            coherfreq(1),cycles);
        return
    end
end

if isnan(timelimits)
    timelimits = [min(times) max(times)];
end
if ~ischar(aligntime) && ~isnan(aligntime)
    if ~isinf(aligntime) ...
            & (aligntime < timelimits(1) | aligntime > timelimits(2))
        help erpimage
        fprintf('\nerpimage(): requested align time outside of time limits.\n');
        return
    end
end
%
%% %%%%%%%%%%%%%%  Replace nan's with 0s %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
nans = find(isnan(data));
if length(nans)
    fprintf('Replaced %d nan in data with 0s.\n');
    data(nans) = 0;
end
%
%% %%%%%%%%%%%% Reshape data to (frames,ntrials) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if size(data,2) ~= ntrials
    if size(data,1)>1
        % fprintf('frames %d, ntrials %d length(data) %d\n',frames,ntrials,length(data));
        data=reshape(data,1,frames*ntrials);
    end
    data=reshape(data,frames,ntrials);
end
fprintf('Plotting input data as %d epochs of %d frames sampled at %3.1f Hz.\n',...
    ntrials,frames,srate);
%
%% %%%%%%%%%%%% Reshape data2 to (frames,ntrials) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if exist('data2') == 1
    if size(data2,2) ~= ntrials
        if size(data2,1)>1
            data2=reshape(data2,1,frames*ntrials);
        end
        data2=reshape(data2,frames,ntrials);
    end
end
%
%% %%%%%%%%%%%%% if sortvar=NaN, remove lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ad
%
if any(isnan(sortvar))
    nanlocs = find(isnan(sortvar));
    fprintf('Removing %d trials with NaN sortvar values.\n', length(nanlocs));
    data(:,nanlocs) = [];
    sortvar(nanlocs) = [];
    if exist('data2') == 1
        data2(:,nanlocs) = [];
    end
    if ~isempty(auxvar)
        auxvar(:,nanlocs) = [];
    end
    if ~isempty(verttimes)
        if size(verttimes,1) == ntrials
            verttimes(nanlocs,:) = [];
        end
    end
    ntrials = size(data,2);
    if ntrials <= 1, close(gcf); error('\nerpimage(): Too few trials'); end
end

%% Create moving average window %%
if strcmpi(mvavg_type,'Gaussian'),
    %construct Gaussian window to weight trials
    if avewidth == 0,
        avewidth = DEFAULT_SDEV;
    elseif avewidth < 1,
        help erpimage
        fprintf('\nerpimage(): Variable avewidth cannot be < 1.\n')
        fprintf('\nerpimage(): avewidth needs to be a positive integer.\n')
        return
    end
    wt_wind=exp(-0.5*([-3*avewidth:3*avewidth]/avewidth).^2)';
    wt_wind=wt_wind/sum(wt_wind); %normalize to unit sum
    avewidth=length(wt_wind);
    
    if avewidth > ntrials
        avewidth=floor((ntrials-1)/6);
        if avewidth==0,
            avewidth=DEFAULT_SDEV; %should be a window with one time point (smallest possible)
        end
        wt_wind=exp(-0.5*([-3*avewidth:3*avewidth]/avewidth).^2)';
        wt_wind=wt_wind/sum(wt_wind);
        fprintf('avewidth is too big for this number of trials.\n');
        fprintf('Changing avewidth to maximum possible size: %d\n',avewidth);
        avewidth=length(wt_wind);
    end
else
    %construct rectangular "boxcar" window to equally weight trials within
    %window
    if avewidth == 0,
        avewidth = DEFAULT_AVEWIDTH;
    elseif avewidth < 1
        help erpimage
        fprintf('\nerpimage(): Variable avewidth cannot be < 1.\n')
        return
    elseif avewidth > ntrials
        fprintf('Setting variable avewidth to max %d.\n',ntrials)
        avewidth = ntrials;
    end
    wt_wind=ones(1,avewidth)/avewidth;
end


%% Filter data with Butterworth filter (if requested) %%%%%%%%%%%%
%
if ~isempty(flt)
    %error check
    if length(flt)~=2,
        error('\nerpimage(): ''filt'' parameter argument should be a two element vector.');
    elseif max(flt)>(srate/2),
        error('\nerpimage(): ''filt'' parameters need to be less than or equal to sampling rate/2 (i.e., %f).',srate/2);
    elseif (flt(2)==(srate/2)) && (flt(1)==0),
        error('\nerpimage(): If second element of ''filt'' parameter is srate/2, then the first element must be greater than 0.');
    elseif abs(flt(2))<=abs(flt(1)),
        error('\nerpimage(): Second element of ''filt'' parameters must be greater than first in absolute value.');
    elseif (flt(1)<0) || (flt(2)<0),
        if (flt(1)>=0) || (flt(2)>=0),
            error('\nerpimage(): BOTH parameters of ''filt'' need to be greater than or equal to zero OR need to be negative.');
        end
        if min(flt)<=(-srate/2),
            error('\nerpimage(): ''filt'' parameters need to be greater than sampling rate/2 (i.e., -%f) when creating a stop band.',srate/2);
        end
    end
    
    fprintf('\nFiltering data with 3rd order Butterworth filter: ');
    if (flt(1)==0),
        %lowpass filter the data
        [B A]=butter(3,flt(2)*2/srate,'low');
        fprintf('lowpass at %.0f Hz\n',flt(2));
    elseif (flt(2)==(srate/2)),
        %highpass filter the data
        [B A]=butter(3,flt(1)*2/srate,'high');
        fprintf('highpass at %.0f Hz\n',flt(1));
    elseif (flt(1)<0)
        %bandstop filter the data
        flt=-flt;
        [B A]=butter(3,flt*2/srate,'stop');
        fprintf('stopband from %.0f to %.0f Hz\n',flt(1),flt(2));
    else
        %bandpass filter the data
        [B A]=butter(3,flt*2/srate);
        fprintf('bandpass from %.0f to %.0f Hz\n',flt(1),flt(2));
    end
    s=size(data);
    for trial=1:s(2),
        data(:,trial)=filtfilt(B,A,double(data(:,trial)));
    end
    if isempty(baseline)
        fprintf('Note, you might want to re-baseline the data using the erpimage ''baseline'' option.\n\n');
    end
end

%% Mean Baseline Each Trial (if requested) %%
if ~isempty(baseline),
    %check argument values for errors
    if baseline(2)<baseline(1),
        error('\nerpimage(): First element of ''baseline'' argument needs to be less than or equal to second argument.');
    elseif baseline(2)<times(1),
        error('\nerpimage(): Second element of ''baseline'' argument needs to be greater than or equal to epoch start time %.1f.',times(1));
    elseif baseline(1)>times(end),
        error('\nerpimage(): First element of ''baseline'' argument needs to be less than or equal to epoch end time %.1f.',times(end));
    end
    
    %convert msec into time points
    if baseline(1)<times(1),
        strt_pt=1;
    else
        strt_pt=find_crspnd_pt(baseline(1),times,1:length(times));
        strt_pt=ceil(strt_pt);
    end
    if baseline(end)>times(end),
        end_pt=length(times);
    else
        end_pt=find_crspnd_pt(baseline(2),times,1:length(times));
        end_pt=floor(end_pt);
    end
    fprintf('\nRemoving pre-stimulus mean baseline from %.1f to %.1f msec.\n\n',times(strt_pt),times(end_pt));
    bsln_mn=mean(data(strt_pt:end_pt,:),1);
    data=data-repmat(bsln_mn,length(times),1);
end


%
%% %%%%%%%%%%%%%%%%% Renormalize sortvar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
switch lower(renorm)
    case 'yes',
        disp('\nerpimage warning: *** sorting variable renormalized ***');
        sortvar = (sortvar-min(sortvar)) / (max(sortvar) - min(sortvar)) * ...
            0.5 * (max(times) - min(times)) + min(times) + 0.4*(max(times) - min(times));
    case 'no',;
    otherwise,
        if ~isempty(renorm)
            locx = findstr('x', lower(renorm));
            if length(locx) ~= 1, error('\nerpimage: unrecognized renormalizing formula'); end
            eval( [ 'sortvar =' renorm(1:locx-1) 'sortvar' renorm(locx+1:end) ';'] );
        end
end
%
%% %%%%%%%%%%%%%%%%% Align data to sortvar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

if ischar(aligntime) || ~isnan(aligntime)
    if ~ischar(aligntime) && isinf(aligntime)
        aligntime= median(sortvar);
        fprintf('Aligning data to median sortvar.\n');
        % Alternative below: trimmed median - ignore top/bottom 5%
        %   ssv = sort(sortvar); % ssv = 'sorted sortvar'
        %   aligntime= median(ssv(ceil(ntrials/20)):floor(19*ntrials/20));
    end
    
    if ~ischar(aligntime)
        fprintf('Realigned sortvar plotted at %g ms.\n',aligntime);
        aligndata=zeros(frames,ntrials); % begin with matrix of zeros()
        shifts = zeros(1,ntrials);
        for t=1:ntrials, %%%%%%%%% foreach trial %%%%%%%%%
            shft = round((aligntime-sortvar(t))*srate/1000);
            shifts(t) = shft;
            if shft>0, % shift right
                if frames-shft > 0
                    aligndata(shft+1:frames,t)=data(1:frames-shft,t);
                else
                    fprintf('No aligned data for epoch %d - shift (%d frames) too large.\n',t,shft);
                end
            elseif shft < 0 % shift left
                if frames+shft > 0
                    aligndata(1:frames+shft,t)=data(1-shft:frames,t);
                else
                    fprintf('No aligned data for epoch %d - shift (%d frames) too large.\n',t,shft);
                end
            else % shft == 0
                aligndata(:,t) = data(:,t);
            end
        end % end trial
        if ~isempty(auxvar)
            auxvar = auxvar+shifts;
        end
        fprintf('Shifted epochs by %d to %d frames.\n',min(shifts),max(shifts));
        data = aligndata;                       % now data is aligned to sortvar
    else
        aligntime = str2num(aligntime);
        if isinf(aligntime),  aligntime= median(sortvar); end
    end
end

%
%% %%%%%%%%%%%%%%%%%%%% Remove the ERP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if strcmpi(Rmerp, 'yes')
    data = data - repmat(nan_mean(data')', [1 size(data,2)]);
end

%
%% %%%%%%%%%%%%% Sort the data trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if exist('phargs') == 1 % if phase-sort the data trials
    if length(phargs) >= 4 && phargs(3) ~= phargs(4) % find max frequency
        % in specified band
        if exist('pwelch') == 2 % requires Signal Processing Toolbox
            fprintf('Computing data spectrum using pwelch().\n');
            [pxx,freqs] = pwelch(data(:), frames, 0, max(1024, pow2(ceil(log2(frames)))), srate);
        else % EEGLAB native work-around
            fprintf('Computing data spectrum using spec().\n');
            [pxx,freqs] = spec(data(:),max(1024, pow2(ceil(log2(frames)))),srate,frames,0);
        end
        % gf = gcf; % figure;plot(freqs,pxx); %xx=axis; %axis([phargs(3) phargs(4) xx(3) xx(4)]); %figure(gf);
        pxx = 10*log10(pxx);
        n = find(freqs >= phargs(3) & freqs <= phargs(4));
        if ~length(n)
            freq = (phargs(3)+phargs(4))/2;
        end
        [dummy maxx] = max(pxx(n));
        freq = freqs(n(maxx));
    else
        freq = phargs(3); % else use specified frequency
    end
    fprintf('Sorting trials on phase at %.2g Hz.\n',freq);
    
    [amps, cohers, cohsig, ampsig, allamps, allphs] = ...
        phasecoher(data,length(times),srate,freq,cycles,0, ...
        [], [], timeStretchRef, timeStretchMarks);
    
    phwin = phargs(1);
    [dummy minx] = min(abs(times-phwin)); % closest time to requested
    winlen = floor(cycles*srate/freq);
    winloc = minx-linspace(floor(winlen/2), floor(-winlen/2), winlen+1);
    tmprange = find(winloc>0 & winloc<=frames);
    winloc = winloc(tmprange); % sorting window times
    winlocs = winloc;
    %
    %%%%%%%%%%%%%%%%%%%% Compute phsamp and phaseangles %%%%%%%%%%%%%%%%%%%%
    %
    % $$$     [phaseangles phsamp] = phasedet(data,frames,srate,winloc,freq);
    phaseangles = allphs(minx,:);
    phsamp = allamps(minx,:);
    % $$$     % hist(phaseangles,50);
    % $$$     % return
    % $$$
    % $$$     %
    % $$$     % Print facts on commandline %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % $$$     %
    % $$$     if length(tmprange) ~=  winlen+1
    % $$$       filtersize = cycles * length(tmprange) / (winlen+1);
    % $$$       timecenter = median(winloc)/srate*1000+times(1); % center of window in ms
    % $$$       phaseangles = phaseangles + 2*pi*(timecenter-phargs(1))*freq;
    % $$$       fprintf('Sorting data epochs by phase at frequency %2.1f Hz: \n', freq);
    % $$$       fprintf(...
    % $$$           '    Data time limits reached -> now uses a %1.1f cycles (%1.0f ms) window centered at %1.0f ms\n', ...
    % $$$           filtersize, 1000/freq*filtersize, timecenter);
    % $$$       fprintf(...
    % $$$           '    Filter length is %d; Phase has been linearly interpolated to latency at %1.0f ms.\n', ...
    % $$$           length(winloc), phargs(1));
    % $$$     else
    % $$$       fprintf(...
    % $$$           'Sorting data epochs by phase at %2.1f Hz in a %1.1f-cycle (%1.0f ms) window centered at %1.0f ms.\n',...
    % $$$           freq,cycles,1000/freq*cycles,times(minx));
    % $$$       fprintf('Phase is computed using a wavelet of %d frames.\n',length(winloc));
    % $$$     end
    %
    % Reject small (or large) phsamp trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    phargs(2) = phargs(2)/100; % convert rejection rate from % to fraction
    amprej = phargs(2);
    [tmp ampsortidx] = sort(phsamp); % sort amplitudes
    if amprej>=0
        ampsortidx = ampsortidx(ceil(amprej*length(ampsortidx))+1:end); % if amprej==0, select all trials
        fprintf('Retaining %d epochs (%g percent) with largest power at the analysis frequency,\n',...
            length(ampsortidx),100*(1-amprej));
    else % amprej < 0
        amprej = 1+amprej; % subtract from end
        ampsortidx = ampsortidx(1:floor(amprej*length(ampsortidx)));
        fprintf('Retaining %d epochs (%g percent) with smallest power at the analysis frequency,\n',...
            length(ampsortidx),amprej*100);
    end
    %
    % Remove low|high-amplitude trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    data = data(:,ampsortidx); % amp-sort the data, removing rejected-amp trials
    phsamp = phsamp(ampsortidx);           % amp-sort the amps
    phaseangles = phaseangles(ampsortidx); % amp-sort the phaseangles
    sortvar = sortvar(ampsortidx);         % amp-sort the trial indices
    ntrials = length(ampsortidx);          % number of trials retained
    if ~isempty(auxvar)
        auxvar = auxvar(:,ampsortidx);
    end
    if ~isempty(timeStretchMarks)
        timeStretchMarks =  timeStretchMarks(:,ampsortidx);
    end
    
    %
    % Sort remaining data by phase angle %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    phaseangles = -phaseangles;
    topphase = (topphase/360)*2*pi; % convert from degrees to radians
    ip = find(phaseangles>topphase);
    phaseangles(ip) = phaseangles(ip)-2*pi; % rotate so topphase at top of plot
    
    [phaseangles sortidx] = sort(phaseangles); % sort trials on (rotated) phase
    data    =  data(:,sortidx);                % sort data by phase
    phsamp  =  phsamp(sortidx);                % sort amps by phase
    sortvar = sortvar(sortidx);                % sort input sortvar by phase
    if ~isempty(auxvar)
        auxvar = auxvar(:,sortidx);
    end
    if ~isempty(timeStretchMarks)
        timeStretchMarks =  timeStretchMarks(:,sortidx);
    end
    phaseangles = -phaseangles; % Note: phsangles now descend from pi
    % TEST auxvar = 360 + (1000/256)*(256/5)*phaseangles/(2*pi); % plot phase+360 in ms for test
    
    fprintf('Size of data = [%d,%d]\n',size(data,1),size(data,2));
    sortidx = ampsortidx(sortidx); % return original trial indices in final sorted order
    %
    % %%%%%%%%%%%%%%% Sort data by amplitude %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
elseif exist('ampargs') == 1 % if amplitude-sort
    if length(ampargs) == 4 % find max frequency in specified band
        if exist('pwelch') == 2
            fprintf('Computing data spectrum using pwelch().\n');
            [pxx,freqs] = pwelch(data(:),frames,0,max(1024, pow2(ceil(log2(frames)))),srate);
        else
            fprintf('Computing data spectrum using spec().\n');
            [pxx,freqs] = spec(data(:),max(1024, pow2(ceil(log2(frames)))),srate,frames,0);
        end
        pxx = 10*log10(pxx);
        if ampargs(3) == ampargs(4)
            [freq n] = min(abs(freqs - ampargs(3)));
        else
            n = find(freqs >= abs(ampargs(3)) & freqs <= abs(ampargs(4)));
        end
        if ~length(n)
            freq = mean([abs(ampargs(3)),abs(ampargs(4))]);
        end
        if ampargs(3)>=0
            [dummy maxx] = max(pxx(n));
            freq = freqs(n(maxx)); % use the highest-power frequency
        else
            freq = freqs(n);  % use all frequencies in the specified range
        end
    else
        freq = abs(ampargs(3)); % else use specified frequency
    end
    if length(freq) == 1
        fprintf('Sorting data epochs by amplitude at frequency %2.1f Hz \n', freq);
    else
        fprintf('Sorting data epochs by amplitude at %d frequencies (%2.1f Hz to %.1f Hz) \n',...
            length(freq),freq(1),freq(end));
    end
    SPECWININCR = 10;	% make spectral sorting time windows increment by 10 ms
    if isinf(ampargs(1))
        ampwins = sortwinarg(1):SPECWININCR:sortwinarg(2);
    else
        ampwins = ampargs(1);
    end
    if ~isinf(ampargs(1)) % single time given
        if length(freq) == 1
            fprintf('   in a %1.1f-cycle (%1.0f ms) time window centered at %1.0f ms.\n',...
                cycles,1000/freq(1)*cycles,ampargs(1));
        else
            fprintf('   in %1.1f-cycle (%1.0f-%1.0f ms) time windows centered at %1.0f ms.\n',...
                cycles,1000/freq(1)*cycles,1000/freq(end)*cycles,ampargs(1));
        end
    else % range of times
        [dummy sortwin_st ] = min(abs(times-ampwins(1)));
        [dummy sortwin_end] = min(abs(times-ampwins(end)));
        if length(freq) == 1
            fprintf('   in %d %1.1f-cycle (%1.0f ms) time windows centered from %1.0f to  %1.0f ms.\n',...
                length(ampwins),cycles,1000/freq(1)*cycles,times(sortwin_st),times(sortwin_end));
        else
            fprintf('   in %d %1.1f-cycle (%1.0f-%1.0f ms) time windows centered from %1.0f to %1.0f ms.\n',...
                length(ampwins),cycles,1000/freq(1)*cycles,1000/freq(end)*cycles,times(sortwin_st),times(sortwin_end));
        end
    end
    
    phsamps = 0; %%%%%%%%%%%%%%%%%%%%%%%%%% sort by (mean) amplitude %%%%%%%%%%%%%%%%%%%%%%%%%%
    minxs = [];
    for f = 1:length(freq)  % use one or range of frequencies
        frq = freq(f);
        [amps, cohers, cohsig, ampsig, allamps, allphs] = ...
            phasecoher(data,length(times),srate,frq,cycles,0, ...
            [], [], timeStretchRef, timeStretchMarks);
        
        for ampwin = ampwins
            [dummy minx] = min(abs(times-ampwin)); % find nearest time point to requested
            minxs = [minxs minx];
            winlen = floor(cycles*srate/frq);
            % winloc = minx-[winlen:-1:0]; % ending time version
            winloc = minx-linspace(floor(winlen/2), floor(-winlen/2), winlen+1);
            tmprange = find(winloc>0 & winloc<=frames);
            winloc = winloc(tmprange); % sorting window frames
            if f==1
                winlocs = [winlocs;winloc];  % store tme windows
            end
            % $$$         [phaseangles phsamp] = phasedet(data,frames,srate,winloc,frq);
            % $$$         phsamps = phsamps+phsamps;  % accumulate amplitudes across 'sortwin'
            phaseangles = allphs(minx,:);
            phsamps = phsamps+allamps(minx,:);  % accumulate amplitudes across 'sortwin'
        end
    end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if length(tmprange) ~=  winlen+1 % ?????????
        filtersize = cycles * length(tmprange) / (winlen+1);
        timecenter = median(winloc)/srate*1000+times(1); % center of window in ms
        phaseangles = phaseangles + 2*pi*(timecenter-ampargs(1))*freq(end);
        fprintf(...
            '    Data time limits reached -> now uses a %1.1f cycles (%1.0f ms) window centered at %1.0f ms\n', ...
            filtersize, 1000/freq(1)*filtersize, timecenter);
        fprintf(...
            '    Wavelet length is %d; Phase has been linearly interpolated to latency et %1.0f ms.\n', ...
            length(winloc(1,:)), ampargs(1));
    end
    if length(freq) == 1
        fprintf('Amplitudes are computed using a wavelet of %d frames.\n',length(winloc(1,:)));
    end
    
    %
    % Reject small (or large) phsamp trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    ampargs(2) = ampargs(2)/100; % convert rejection rate from % to fraction
    [tmp n] = sort(phsamps); % sort amplitudes
    if ampargs(2)>=0
        n = n(ceil(ampargs(2)*length(n))+1:end); % if rej 0, select all trials
        fprintf('Retaining %d epochs (%g percent) with largest power at the analysis frequency,\n',...
            length(n),100*(1-ampargs(2)));
    else % ampargs(2) < 0
        ampargs(2) = 1+ampargs(2); % subtract from end
        n = n(1:floor(ampargs(2)*length(n)));
        fprintf(...
            'Retaining %d epochs (%g percent) with smallest power at the analysis frequency,\n',...
            length(n),ampargs(2)*100);
    end
    %
    % Remove low|high-amplitude trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    
    data = data(:,n); % amp-sort the data, removing rejected-amp trials
    phsamps = phsamps(n);           % amp-sort the amps
    phaseangles = phaseangles(n); % amp-sort the phaseangles
    sortvar = sortvar(n);         % amp-sort the trial indices
    ntrials = length(n);          % number of trials retained
    if ~isempty(auxvar)
        auxvar = auxvar(:,n);
    end
    %
    % Sort remaining data by amplitude %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    [phsamps sortidx] = sort(phsamps); % sort trials on amplitude
    data    =  data(:,sortidx);                % sort data by amp
    phaseangles  =  phaseangles(sortidx);      % sort angles by amp
    sortvar = sortvar(sortidx);                % sort input sortvar by amp
    if ~isempty(auxvar)
        auxvar = auxvar(:,sortidx);
    end
    fprintf('Size of data = [%d,%d]\n',size(data,1),size(data,2));
    
    sortidx = n(sortidx); % return original trial indices in final sorted order
    %
    %%%%%%%%%%%%%%%%%%%%%% Don't Sort trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
elseif Nosort == YES
    fprintf('Not sorting data on input sortvar.\n');
    sortidx = 1:ntrials;
    %
    %%%%%%%%%%%%%%%%%%%%%% Sort trials on (mean) value %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
elseif exist('valargs')
    [sttime stframe] = min(abs(times-valargs(1)));
    sttime = times(stframe);
    if length(valargs)>1
        [endtime endframe] = min(abs(times-valargs(2)));
        endtime = times(endframe);
    else
        endframe = stframe;
        endtime = times(endframe);
    end
    if length(valargs)==1 || sttime == endtime
        fprintf('Sorting data on value at time %4.0f ms.\n',sttime);
    elseif length(valargs)>1
        fprintf('Sorting data on mean value between %4.0f and %4.0f ms.\n',...
            sttime,endtime);
    end
    if endframe>stframe
        sortval = mean(data(stframe:endframe,:));
    else
        sortval = data(stframe,:);
    end
    [sortval,sortidx] = sort(sortval);
    if length(valargs)>2
        if valargs(3) <0
            sortidx = sortidx(end:-1:1); % plot largest values on top
            % if direction < 0
        end
    end
    data = data(:,sortidx);
    sortvar = sortvar(sortidx);                % sort input sortvar by amp
    if ~isempty(auxvar)
        auxvar = auxvar(:,sortidx);
    end
    if ~isempty(phaseangles)
        phaseangles  =  phaseangles(sortidx);      % sort angles by amp
    end
    winloc = [stframe,endframe];
    winlocs = winloc;
    %
    %%%%%%%%%%%%%%%%%%%%%% Sort trials on sortvar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
else
    fprintf('Sorting data on input sortvar.\n');
    [sortvar,sortidx] = sort(sortvar);
    data = data(:,sortidx);
    if ~isempty(auxvar)
        auxvar = auxvar(:,sortidx);
    end
    uni_svar=unique_bc(sortvar);
    n_ties=0;
    tie_dist=zeros(1,length(uni_svar));
    loop_ct=0;
    for tie_loop=uni_svar,
        ids=find(sortvar==tie_loop);
        n_ids=length(ids);
        if n_ids>1,
            if replace_ties==YES,
                mn=mean(data(:,ids),2);
                data(:,ids)=repmat(mn,1,n_ids);
                if ~isempty(auxvar)
                    mn=mean(auxvar(:,ids),2);
                    auxvar(:,ids) = repmat(mn,1,n_ids);
                end
            end
            n_ties=n_ties+n_ids;
            loop_ct=loop_ct+1;
            tie_dist(loop_ct)=n_ids;
        end
    end
    fprintf('%.2f%c of the trials (i.e., %d out of %d) have the same sortvar value as at least one other trial.\n', ...
        100*n_ties/length(sortvar),37,n_ties,length(sortvar));
    fprintf('Distribution of number ties per unique value of sortvar:\n');
    if exist('prctile')
        try
            fprintf('Min: %d, 25th ptile: %d, Median: %d, 75th ptile: %d, Max: %d\n',min(tie_dist),round(prctile(tie_dist,25)), ...
                round(median(tie_dist)),round(prctile(tie_dist,75)),max(tie_dist));
        catch
        end
    end
    if replace_ties==YES,
        fprintf('Trials with tied sorting values will be replaced by their mean.\n');
    end
    fprintf('\n');
end

%if max(sortvar)<0
%   fprintf('Changing the sign of sortvar: making it positive.\n');
%   sortvar = -sortvar;
%end
%
%% %%%%%%%%%%%%%%%%% Adjust decfactor if phargs or ampargs %%%%%%%%%%%%%%%%%%%%%
%
if decfactor < 0
    decfactor = -decfactor;
    invdec = 1;
else
    invdec = 0;
end
if decfactor > sqrt(ntrials) % if large, output this many trials
    n = 1:ntrials;
    if exist('phargs') && length(phargs)>1
        if phargs(2)>0
            n = n(ceil(phargs(2)*ntrials)+1:end); % trials after rejection
        elseif phargs(2)<0
            n = n(1:floor(phargs(2)*length(n)));  % trials after rejection
        end
    elseif exist('ampargs') && length(ampargs)>1
        if ampargs(2)>0
            n = n(ceil(ampargs(2)*ntrials)+1:end); % trials after rejection
        elseif ampargs(2)<0
            n = n(1:floor(ampargs(2)*length(n)));  % trials after rejection
        end
    end
    if invdec
        decfactor = (length(n)-avewidth)/decfactor;
    else
        decfactor = length(n)/decfactor;
    end
end
if ~isreal(decfactor), decfactor = imag(decfactor); end

%
%% %%%%%%%%%%%%%%%% Smooth data using moving average %%%%%%%%%%%%%%%%%%%%%%%%%%%
%
urdata = data; % save data to compute amp, coher on unsmoothed data

if ~Allampsflag && ~exist('data2') % if imaging potential,
    if length(timeStretchRef) > 0 && length(timeStretchMarks) > 0
        %
        % Perform time-stretching here -JH %%%%%%%%%%%%%%%%
        %
        for t=1:size(data,2)
            M = timewarp(timeStretchMarks(t,:)', timeStretchRef');
            data(:,t) = M*data(:,t);
        end
    end
    tsurdata = data;
    % Time-stretching ends here %%%%%%%%%%%%%%
    
    if avewidth > 1 || decfactor > 1
        if Nosort == YES
            fprintf('Smoothing the data using a window width of %g epochs ',avewidth);
        else
            fprintf('Smoothing the sorted epochs with a %g-epoch moving window.',...
                avewidth);
        end
        fprintf('\n');
        fprintf('  and a decimation factor of %g\n',decfactor);
        
        if ~exist('phargs') % if not phase-sorted trials
            [data,outtrials] = movav(data,1:ntrials,avewidth,decfactor,[],[],wt_wind);
            % Note: movav() here sorts using square window
            [outsort,outtrials] = movav(sortvar,1:ntrials,avewidth,decfactor,[],[],wt_wind);
            
        else % if phase-sorted trials, use circular / wrap-around smoothing
            backhalf  = floor(avewidth/2);
            fronthalf = floor((avewidth-1)/2);
            if avewidth > 2
                [data,outtrials] = movav([data(:,[(end-backhalf+1):end]),...
                    data,...
                    data(:,[1:fronthalf])],...
                    [1:(ntrials+backhalf+fronthalf)],avewidth,decfactor,[],[],wt_wind);
                [outsort,outtrials] = movav([sortvar((end-backhalf+1):end),...
                    sortvar,...
                    sortvar(1:fronthalf)],...
                    1:(ntrials+backhalf+fronthalf),avewidth,decfactor,[],[],wt_wind);
                % Shift elements of outtrials so the first element is 1
                outtrials = outtrials - outtrials(1) + 1;
            else % avewidth==2
                [data,outtrials] = movav([data(:,end),data],...
                    [1:(ntrials+1)],avewidth,decfactor,[],[],wt_wind);
                % Note: movav() here sorts using square window
                [outsort,outtrials] = movav([sortvar(end) sortvar],...
                    1:(ntrials+1),avewidth,decfactor,[],[],wt_wind);
                % Shift elements of outtrials so the first element is 1
                outtrials = outtrials - outtrials(1) + 1;
            end
        end
        for index=1:length(percentiles)
            outpercent{index} = compute_percentile( sortvar, percentiles(index), outtrials, avewidth);
        end
        if ~isempty(auxvar)
            if ~exist('phargs') % if not phase-sorted trials
                [auxvar,tmp] = movav(auxvar,1:ntrials,avewidth,decfactor,[],[],wt_wind);
            else % if phase-sorted trials
                if avewidth>2
                    [auxvar,tmp] = movav([auxvar(:,[(end-backhalf+1):end]),...
                        auxvar,...
                        auxvar(:,[1:fronthalf])],...
                        [1:(ntrials+backhalf+fronthalf)],avewidth,decfactor,[],[],wt_wind);
                    % Shift elements of tmp so the first element is 1
                    tmp = tmp - tmp(1) + 1;
                else % avewidth==2
                    [auxvar,tmp] = movav([auxvar(:,end),auxvar],[1:(ntrials+1)],avewidth,decfactor,[],[],wt_wind);
                    % Shift elements of tmp so the first element is 1
                    tmp = tmp - tmp(1) + 1;
                end
            end
        end
        %  if ~isempty(sortvar_limits),
        %      fprintf('Output data will be %d frames by %d smoothed trials.\n',...
        %          frames,length(outtrials));
        %      fprintf('Outtrials: %3.2f to %4.2f\n',min(outtrials),max(outtrials));
        %  end
    else % don't smooth
        outtrials = 1:ntrials;
        outsort = sortvar;
    end
    
    %
    %%%%%%%%%%%%%%%%%%%%%%%%% Find color axis limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    if ~isempty(Caxis)
        mindat = Caxis(1);
        maxdat = Caxis(2);
        fprintf('Using the specified caxis range of [%g,%g].\n', mindat, maxdat);
    else
        mindat = min(min(data));
        maxdat = max(max(data));
        maxdat =  max(abs([mindat maxdat])); % make symmetrical about 0
        mindat = -maxdat;
        if ~isempty(caxfraction)
            adjmax = (1-caxfraction)/2*(maxdat-mindat);
            mindat = mindat+adjmax;
            maxdat = maxdat-adjmax;
            fprintf(...
                'The caxis range will be %g times the sym. abs. data range -> [%g,%g].\n',...
                caxfraction,mindat,maxdat);
        else
            fprintf(...
                'The caxis range will be the sym. abs. data range -> [%g,%g].\n',...
                mindat,maxdat);
        end
    end
end % if ~Allampsflag & ~exist('data2')

%
%% %%%%%%%%%%%%%%%%%%%%%%%% Set time limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if isnan(timelimits(1))
    timelimits = [min(times) max(times)];
end
fprintf('Data will be plotted between %g and %g ms.\n',timelimits(1),timelimits(2));

%
%% %%%%%%%%%%% Image the aligned/sorted/smoothed data %%%%%%%%%%%%%%%%%%%%%%%%%%
%
if strcmpi(NoShow, 'no')
    if ~any(isnan(coherfreq))       % if plot three time axes
        image_loy = 3*PLOT_HEIGHT;
    elseif Erpflag == YES   % elseif if plot only one time axes
        image_loy = 1*PLOT_HEIGHT;
    else                    % else plot erp-image only
        image_loy = 0*PLOT_HEIGHT;
    end
    gcapos=get(gca,'Position');
    delete(gca)
    if isempty(topomap)
        image_top = 1;
    else
        image_top = 0.9;
    end
    ax1=axes('Position',...
        [gcapos(1) gcapos(2)+image_loy*gcapos(4) ...
        gcapos(3) (image_top-image_loy)*gcapos(4)]);
end
ind = isnan(data);    % find nan's in data
[i j]=find(ind==1);
if ~isempty(i)
    data(i,j) = 0;      % plot shifted nan data as 0 (=green)
end

%
%% %%%%%%%%%%% Determine coherence frequency %%%%%%%%%%%%%%%%%%%%%%%%%%
%
if length(coherfreq) == 2 && coherfreq(1) ~= coherfreq(2) && freq <= 0
    % find max frequency in specified band - should use Matlab pwelch()?
    if exist('pwelch') == 2 % from Signal Processing Toolbox
        [pxx,tmpfreq] = pwelch(urdata(:),frames,0,max(1024,pow2(ceil(log2(frames)))),srate);
    else % substitute from EEGLAB
        [pxx,tmpfreq] = spec(urdata(:),max(1024,pow2(ceil(log2(frames)))),srate,frames,0);
    end
    pxx = 10*log10(pxx);
    n = find(tmpfreq >= coherfreq(1) & tmpfreq <= coherfreq(2));
    % [tmpfreq(n) pxx(n)]
    % coherfreqs = coherfreq; % save for debugging spectrum plotting
    
    if ~length(n)
        coherfreq = coherfreq(1);
    end
    [dummy maxx] = max(pxx(n));
    coherfreq = tmpfreq(n(maxx));
else
    coherfreq = coherfreq(1);
end


if ~Allampsflag && ~exist('data2') %%%%%%%% Plot ERP image %%%%%%%%%%
    
    %Stretch the data array
    % $$$     keyboard;
    
    if strcmpi(NoShow, 'no')
        if TIMEX
            h_eim=imagesc(times,outtrials,data',[mindat,maxdat]);% plot time on x-axis
            set(gca,'Ydir','normal', 'tag', 'erpimage');
            axis([timelimits(1) timelimits(2) ...
                min(outtrials) max(outtrials)]);
        else
            h_eim=imagesc(outtrials,times,data,[mindat,maxdat]); % plot trials on x-axis
            set(gca, 'tag', 'erpimage');
            axis([min(outtrials) max(outtrials)...
                timelimits(1) timelimits(2)]);
        end
        try colormap(DEFAULT_COLORMAP); catch, end
        hold on
        drawnow
    end
    
elseif Allampsflag %%%%%%%%%%%%%%%% Plot allamps instead of data %%%%%%%%%%%%%%
    
    if freq > 0
        coherfreq = mean(freq); % use phase-sort frequency
    end
    
    if ~isnan(signifs) % plot received significance levels
        fprintf(['Computing and plotting received ERSP and ITC signif. ' ...
            'levels...\n']);
        [amps,cohers,cohsig,ampsig,allamps] = ...
            phasecoher(urdata,length(times),srate,coherfreq,cycles,0, ...
            [], [], timeStretchRef, timeStretchMarks);
        % Note: need to receive cohsig and ampsig to get allamps <---
        ampsig = signifs([1 2]); % assume these already in dB
        cohsig = signifs(3);
        
    elseif alpha>0 % compute significance levels
        fprintf('Computing and plotting %g ERSP and ITC signif. level...\n',alpha);
        [amps,cohers,cohsig,ampsig,allamps] = ...
            phasecoher(urdata,length(times),srate,coherfreq, ...
            cycles, alpha, [], [], ...
            timeStretchRef, timeStretchMarks');
        % Note: need to receive cohsig and ampsig to get allamps
        fprintf('Coherence significance level: %g\n',cohsig);
        
    else % no plotting of significance
        [amps,cohers,cohsig,ampsig,allamps] = ...
            phasecoher(urdata,length(times),srate,coherfreq, ...
            cycles,0,[], [], timeStretchRef, timeStretchMarks);
        % Note: need to receive cohsig and ampsig to get allamps
    end
    
    % fprintf('#1 Size of allamps = [%d %d]\n',size(allamps,1),size(allamps,2));
    
    base = find(times<=DEFAULT_BASELINE_END);
    if length(base)<2
        base = 1:floor(length(times)/4); % default first quarter-epoch
    end
    fprintf('Using %g to %g ms as amplitude baseline.\n',...
        times(1),times(base(end)));
    
    % fprintf('#2 Size of allamps = [%d %d]\n',size(allamps,1),size(allamps,2));
    
    % fprintf('Subtracting the mean baseline log amplitude \n');
    %fprintf('Subtracting the mean baseline log amplitude %g\n',baseall);
    % allamps = allamps./baseall;
    % fprintf('#3 Size of allamps = [%d %d]\n',size(allamps,1),size(allamps,2));
    
    if avewidth > 1 || decfactor > 1
        if Nosort == YES
            fprintf(...
                'Smoothing the amplitude epochs using a window width of %g epochs ',...
                avewidth);
        else % sort trials
            fprintf(...
                'Smoothing the sorted amplitude epochs with a %g-epoch moving window.',...
                avewidth);
        end
        fprintf('\n');
        fprintf('  and a decimation factor of %g\n',decfactor);
        
        %fprintf('4 Size of allamps = [%d %d]\n',size(allamps,1),size(allamps,2));
        
        if exist('phargs') % if phase-sorted trials, use circular/wrap-around smoothing
            backhalf  = floor(avewidth/2);
            fronthalf = floor((avewidth-1)/2);
            if avewidth > 2
                [allamps,outtrials] = movav([allamps(:,[(end-backhalf+1):end]),...
                    allamps,...
                    allamps(:,[1:fronthalf])],...
                    [1:(ntrials+backhalf+fronthalf)],avewidth,decfactor,[],[],wt_wind);
                % Note: sort using square window
                [outsort,outtrials] = movav([sortvar((end-backhalf+1):end),...
                    sortvar,...
                    sortvar(1:fronthalf)],...
                    1:(ntrials+backhalf+fronthalf),avewidth,decfactor,[],[],wt_wind);
                % Shift elements of outtrials so the first element is 1
                outtrials = outtrials - outtrials(1) + 1;
                if ~isempty(auxvar)
                    [auxvar,tmp] = movav([auxvar(:,[(end-backhalf+1):end]),...
                        auxvar,...
                        auxvar(:,[1:fronthalf])],...
                        [1:(ntrials+backhalf+fronthalf)],avewidth,decfactor,[],[],wt_wind);
                    % Shift elements of outtrials so the first element is 1
                    outtrials = outtrials - outtrials(1) + 1;
                end
            else % avewidth==2
                [allamps,outtrials] = movav([allamps(:,end),allamps],...
                    [1:(ntrials+1)],avewidth,decfactor,[],[],wt_wind);
                % Note: sort using square window
                [outsort,outtrials] = movav([sortvar(end) sortvar],...
                    1:(ntrials+1),avewidth,decfactor,[],[],wt_wind);
                % Shift elements of outtrials so the first element is 1
                outtrials = outtrials - outtrials(1) + 1;
                [auxvar,tmp] = movav([auxvar(:,end),auxvar],[1:(ntrials+1)],avewidth,decfactor,[],[],wt_wind);
                % Shift elements of tmp so the first element is 1
                tmp = tmp - tmp(1) + 1;
            end
        else % if trials not phase sorted, no wrap-around
            [allamps,outtrials] = movav(allamps,1:ntrials,avewidth,decfactor,[],[],wt_wind);
            %fprintf('5 Size of allamps = [%d %d]\n',size(allamps,1),size(allamps,2));
            [outsort,outtrials] = movav(sortvar,1:ntrials,avewidth,decfactor,[],[],wt_wind);
            if ~isempty(auxvar)
                [auxvar,tmp] = movav(auxvar,1:ntrials,avewidth,decfactor,[],[],wt_wind);
            end
        end
        for index=1:length(percentiles)
            outpercent{index} = compute_percentile( sortvar, percentiles(index), outtrials, avewidth);
        end
        fprintf('Output allamps data will be %d frames by %d smoothed trials.\n',...
            frames,length(outtrials));
        
    else % if no smoothing
        outtrials = 1:ntrials;
        outsort = sortvar;
    end
    
    allamps = 20*log10(allamps); % convert allamps to dB
    amps = 20*log10(amps);       % convert latency mean amps to dB
    ampsig = 20*log10(ampsig);   % convert amplitude signif thresholds to dB
    
    if alpha>0
        fprintf('Amplitude significance levels: [%g %g] dB\n',ampsig(1),ampsig(2));
    end
    
    if isnan(baseamp) % if not specified in 'limits'
        [amps,baseamp] = rmbase(amps,length(times),base); % subtract the dB baseline (baseamp)
        % amps are the means at each latency
        allamps = allamps - baseamp; % subtract dB baseline from allamps
        % amplitude
        ampsig = ampsig - baseamp; % subtract dB baseline from ampsig
        
    else % if baseamp specified in 'limits' (as last argument 'basedB', see help)
        amps = amps-baseamp; % use specified (log) baseamp
        allamps = allamps - baseamp; % subtract dB baseline
        if isnan(signifs);
            ampsig = ampsig-baseamp; % subtract dB baseline
        end
    end
    
    
    %
    %%%%%%%%%%%%%%%%%%%%%%%%% Find color axis limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    if ~isempty(Caxis)
        mindat = Caxis(1);
        maxdat = Caxis(2);
        fprintf('Using the specified caxis range of [%g,%g].\n',...
            mindat,maxdat);
    else
        
        % Changed -JH
        % 0.8 is Scott's suggestion to make the erp image show small
        % variations better
        maxdat = 0.8 * max(max(abs(allamps)));
        mindat = -maxdat;
        % $$$         mindat = min(min(allamps));
        % $$$         maxdat = max(max(allamps));
        % $$$         maxdat =  max(abs([mindat maxdat])); % make symmetrical about 0
        % $$$         mindat = -maxdat;
        
        if ~isempty(caxfraction)
            adjmax = (1-caxfraction)/2*(maxdat-mindat);
            mindat = mindat+adjmax;
            maxdat = maxdat-adjmax;
            fprintf(...
                'The caxis range will be %g times the sym. abs. data range -> [%g,%g].\n',...
                caxfraction,mindat,maxdat);
        else
            fprintf(...
                'The caxis range will be the sym. abs. data range -> [%g,%g].\n',...
                mindat,maxdat);
        end
    end
    %
    %%%%%%%%%%%%%%%%%%%%% Image amplitudes at coherfreq %%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    
    if strcmpi(NoShow, 'no')
        fprintf('Plotting amplitudes at freq %g Hz instead of potentials.\n',coherfreq);
        if TIMEX
            imagesc(times,outtrials,allamps',[mindat,maxdat]);% plot time on x-axis
            set(gca,'Ydir','normal');
            axis([timelimits(1) timelimits(2) ...
                min(outtrials) max(outtrials)]);
        else
            imagesc(outtrials,times,allamps,[mindat,maxdat]); % plot trials on x-axis
            axis([min(outtrials) max(outtrials)...
                timelimits(1) timelimits(2)]);
        end
        try colormap(DEFAULT_COLORMAP); catch, end
        drawnow
        hold on
    end
    data = allamps;
    
elseif exist('data2') %%%%%% Plot allcohers instead of data %%%%%%%%%%%%%%%%%%%
    %%%%%%%%% UNDOCUMENTED AND DEPRECATED OPTION %%%%%%%%%%%%
    if freq > 0
        coherfreq = mean(freq); % use phase-sort frequency
    end
    if alpha>0
        fprintf('Computing and plotting %g coherence significance level...\n',alpha);
        
        % [amps,cohers,cohsig,ampsig,allcohers] = ...
        %   crosscoher(urdata,data2,length(times),srate,coherfreq,cycles,alpha);
        
        fprintf('Inter-Trial Coherence significance level: %g\n',cohsig);
        fprintf('Amplitude significance levels: [%g %g]\n',ampsig(1),ampsig(2));
    else
        % [amps,cohers,cohsig,ampsig,allcohers] = ...
        %    crosscoher(urdata,data2,length(times),srate,coherfreq,cycles,0);
    end
    if ~exist('allcohers')
        fprintf('\nerpimage(): allcohers not returned....\n')
        return
    end
    allamps = allcohers; % output variable
    % fprintf('Size allcohers = (%d, %d)\n',size(allcohers,1),size(allcohers,2));
    % fprintf('#1 Size of allcohers = [%d %d]\n',size(allcohers,1),size(allcohers,2));
    base = find(times<=0);
    if length(base)<2
        base = 1:floor(length(times)/4); % default first quarter-epoch
    end
    
    amps = 20*(log10(amps) - log10(mean(amps))); % convert to dB
    %amps = 20*log10(amps); % convert to dB
    ampsig = 20*(log10(ampsig) - log10(mean(amps))); % convert to dB
    %ampsig = 20*log10(ampsig); % convert to dB
    
    if isnan(baseamp)
        [amps,baseamp] = rmbase(amps,length(times),base); % remove baseline
    else
        amps = amps - baseamp;
    end
    % fprintf('#2 Size of allcohers = [%d %d]\n',size(allcohers,1),size(allcohers,2));
    
    if avewidth > 1 || decfactor > 1
        if Nosort == YES
            fprintf(...
                'Smoothing the amplitude epochs using a window width of %g epochs '...
                ,avewidth);
        else
            fprintf(...
                'Smoothing the sorted amplitude epochs with a %g-epoch moving window.'...
                ,avewidth);
        end
        fprintf('\n');
        fprintf('  and a decimation factor of %g\n',decfactor);
        % fprintf('4 Size of allcohers = [%d %d]\n',size(allcohers,1),size(allcohers,2));
        
        [allcohers,outtrials] = movav(allcohers,1:ntrials,avewidth,decfactor,[],[],wt_wind);
        % fprintf('5 Size of allcohers = [%d %d]\n',size(allcohers,1),size(allcohers,2));
        [outsort,outtrials] = movav(sortvar,1:ntrials,avewidth,decfactor,[],[],wt_wind);
        %        fprintf('Output data will be %d frames by %d smoothed trials.\n',...
        %            frames,length(outtrials));
    else
        outtrials = 1:ntrials;
        outsort = sortvar;
    end
    %
    %%%%%%%%%%%%%%%%%%%%%%%%% Find color axis limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    if ~isempty(Caxis)
        mindat = Caxis(1);
        maxdat = Caxis(2);
        fprintf('Using the specified caxis range of [%g,%g].\n',...
            mindat,maxdat);
    else
        mindat = -1;
        maxdat = 1
        fprintf(...
            'The caxis range will be the sym. abs. data range [%g,%g].\n',...
            mindat,maxdat);
    end
    %
    %%%%%%%%%%%%%%%%%%%%% Image coherences at coherfreq %%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    if strcmpi(NoShow, 'no')
        fprintf('Plotting coherences at freq %g Hz instead of potentials.\n',coherfreq);
        if TIMEX
            imagesc(times,outtrials,allcohers',[mindat,maxdat]);% plot time on x-axis
            set(gca,'Ydir','normal');
            axis([timelimits(1) timelimits(2) ...
                min(outtrials) max(outtrials)]);
        else
            imagesc(outtrials,times,allcohers,[mindat,maxdat]); % plot trials on x-axis
            axis([min(outtrials) max(outtrials)...
                timelimits(1) timelimits(2)]);
        end
        try colormap(DEFAULT_COLORMAP); catch, end
        drawnow
        hold on
    end
    
    
end

%Change limits on ERPimage y-axis if requested
if ~isempty(sortvar_limits)
    if exist('phargs','var'),
        fprintf('********* Warning *********\n');
        fprintf('Specifying sorting variable limits has no effect when sorting by phase.\n');
    elseif exist('valargs','var'),
        fprintf('********* Warning *********\n');
        fprintf('Specifying sorting variable limits has no effect when sorting by mean EEG voltage.\n');
    elseif exist('ampargs','var'),
        fprintf('********* Warning *********\n');
        fprintf('Specifying sorting variable limits has no effect when sorting by frequency amp.\n');
    else
        v=axis;
        img_mn=find_crspnd_pt(sortvar_limits(1),outsort,outtrials);
        if isempty(img_mn),
            img_mn=1;
            sortvar_limits(1)=outsort(1);
        end
        img_mx=find_crspnd_pt(sortvar_limits(2),outsort,outtrials);
        if isempty(img_mx),
            img_mx=length(outsort);
            sortvar_limits(2)=outsort(img_mx);
        end
        axis([v(1:2) img_mn img_mx]);
        id1=find(sortvar>=sortvar_limits(1));
        id2=find(sortvar<=sortvar_limits(2));
        id=intersect_bc(id1,id2);
        fprintf('%d epochs fall within sortvar limits.\n',length(id));
        urdata=urdata(:,id);
        if ~isempty(tsurdata),
            tsurdata=tsurdata(:,id);
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%% End plot image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmpi(NoShow, 'no')
    v=axis;
    fprintf('Output data will be %d frames by %d smoothed trials.\n',...
        frames,v(4)-v(3)+1);
    fprintf('Outtrials: %3.2f to %4.2f\n',v(3),v(4));
end

%
%%%%%%%%%%%%%%%%%%%%% Compute y-axis tick values and labels (if requested) %%%%%%%%%%%%%%%%%%%%%%%%%%
%

if strcmpi(NoShow, 'no')
    if ~isempty(img_ylab) && ~strcmpi(img_ylab,'Trials')
        %make ERPimage y-tick labels in units of sorting variable
        if isempty(sortvar_limits),
            mn=min(outsort);
            mx=max(outsort);
        else
            mn=sortvar_limits(1);
            mx=sortvar_limits(2);
        end
        ord=orderofmag(mx-mn);
        rng_rnd=round([mn mx]/ord)*ord;
        if isempty(img_ytick_lab)
            img_ytick_lab=[rng_rnd(1):ord:rng_rnd(2)];
            in_range=find((img_ytick_lab>=mn) & (img_ytick_lab<=mx));
            img_ytick_lab=img_ytick_lab(in_range);
        else
            img_ytick_lab=unique_bc(img_ytick_lab); %make sure it is sorted
            in_range=find((img_ytick_lab>=mn) & (img_ytick_lab<=mx));
            if length(img_ytick_lab)~=length(in_range),
                fprintf('\n***Warning***\n');
                fprintf('''img_trialax_ticks'' exceed smoothed sorting variable values. Max/min values are %f/%f.\n\n',mn,mx);
                img_ytick_lab=img_ytick_lab(in_range);
            end
        end
        n_tick=length(img_ytick_lab);
        img_ytick=zeros(1,n_tick);
        for tickloop=1:n_tick,
            img_ytick(tickloop)=find_crspnd_pt(img_ytick_lab(tickloop),outsort,outtrials);
        end
    elseif ~isempty(img_ylab), %make ERPimage y-tick labels in units of Trials
        if isempty(img_ytick_lab)
            v=axis; %note: sorting variable limits have already been used to determine range of ERPimage y-axis
            mn=v(3);
            mx=v(4);
            ord=orderofmag(mx-mn);
            rng_rnd=round([mn mx]/ord)*ord;
            img_ytick=[rng_rnd(1):ord:rng_rnd(2)];
            in_range=find((img_ytick>=mn) & (img_ytick<=mx));
            img_ytick=img_ytick(in_range);
        else
            img_ytick=img_ytick_lab;
        end
    end
end

%% %%%%%%%%%%%%%%%%%%%%%%%%% plot vert lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isempty(verttimes)
    if size(verttimes,1) ~= 1 && size(verttimes,2) == 1 && size(verttimes,1) ~= ntrials
        verttimes = verttimes';
    end
    if size(verttimes,1) ~= 1 && size(verttimes,1) ~= ntrials
        fprintf('\nerpimage(): vert arg matrix must have 1 or %d rows\n',ntrials);
        return
    end
    if strcmpi(NoShow, 'no')
        if size(verttimes,1) == 1
            fprintf('Plotting %d lines at times: ',size(verttimes,2));
        else
            fprintf('Plotting %d traces starting at times: ',size(verttimes,2));
        end
        for vt = verttimes % for each column
            fprintf('%g ',vt(1));
            if isnan(aligntime) % if NOT re-aligned data
                if TIMEX          % overplot vt on image
                    if length(vt)==1
                        mydotstyle = DOTSTYLE;
                        if exist('auxcolors') && ...
                                length(verttimes) == length(auxcolors)
                            mydotstyle = auxcolors{find(verttimes == vt)};
                        end
                        plot([vt vt],[0 max(outtrials)],mydotstyle,'Linewidth',VERTWIDTH);
                    elseif length(vt)==ntrials
                        [outvt,ix] = movav(vt,1:ntrials,avewidth,decfactor,[],[],wt_wind);
                        plot(outvt,outtrials,DOTSTYLE,'Linewidth',VERTWIDTH);
                    end
                else
                    if length(vt)==1
                        plot([0 max(outtrials)],[vt vt],DOTSTYLE,'Linewidth',VERTWIDTH);
                    elseif length(vt)==ntrials
                        [outvt,ix] = movav(vt,1:ntrials,avewidth,decfactor,[],[],wt_wind);
                        plot(outtrials,outvt,DOTSTYLE,'Linewidth',VERTWIDTH);
                    end
                end
            else % re-aligned data
                if TIMEX          % overplot vt on image
                    if length(vt)==ntrials
                        [outvt,ix] = movav(vt,1:ntrials,avewidth,decfactor,[],[],wt_wind);
                        plot(aligntime+outvt-outsort,outtrials,DOTSTYLE,'LineWidth',VERTWIDTH);
                    elseif length(vt)==1
                        plot(aligntime+vt-outsort,outtrials,DOTSTYLE,'LineWidth',VERTWIDTH);
                    end
                else
                    if length(vt)==ntrials
                        [outvt,ix] = movav(vt,1:ntrials,avewidth,decfactor,[],[],wt_wind);
                        plot(outtrials,aligntime+outvt-outsort,DOTSTYLE,'LineWidth',VERTWIDTH);
                    elseif length(vt)==1
                        plot(outtrials,aligntime+vt-outsort,DOTSTYLE,'LineWidth',VERTWIDTH);
                    end
                end
            end
        end
        %end
        fprintf('\n');
    end
end

%
%% %%%%%%%%%%%%%%%%%%%%%%%%% plot horizontal ('horz') lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if ~isempty(horzepochs)
    if size(horzepochs,1) > 1 && size(horzepochs,1) > 1
        fprintf('\nerpimage(): horz arg must be a vector\n');
        return
    end
    if strcmpi(NoShow, 'no')
        if ~isempty(img_ylab) && ~strcmpi(img_ylab,'Trials'),
            %trial axis in units of sorting variable
            mx=max(outsort);
            mn=min(outsort);
            fprintf('Plotting %d lines at epochs corresponding to sorting variable values: ',length(horzepochs));
            for he = horzepochs % for each horizontal line
                fprintf('%g ',he);
                %find trial number corresponding to this value of sorting
                %variable:
                if (he>mn) && (he<mx)
                    he_ep=find_crspnd_pt(he,outsort,outtrials);
                    if TIMEX          % overplot he_ep on image
                        plot([timelimits(1) timelimits(2)],[he_ep he_ep],LINESTYLE,'Linewidth',HORZWIDTH);
                    else
                        plot([he_ep he_ep], [timelimits(1) timelimits(2)],LINESTYLE,'Linewidth',HORZWIDTH);
                    end
                end
            end
            
        else %trial axis in units of trials
            fprintf('Plotting %d lines at epochs: ',length(horzepochs));
            for he = horzepochs % for each horizontal line
                fprintf('%g ',he);
                if TIMEX          % overplot he on image
                    plot([timelimits(1) timelimits(2)],[he he],LINESTYLE,'Linewidth',HORZWIDTH);
                else
                    plot([he he], [timelimits(1) timelimits(2)],LINESTYLE,'Linewidth',HORZWIDTH);
                end
            end
        end
        fprintf('\n');
    end
end
if strcmpi(NoShow, 'no')
    set(gca,'FontSize',TICKFONT)
    hold on;
end
%
%% %%%%%%%%% plot vertical line at 0 or align time %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if strcmpi(NoShow, 'no')
    if ~isnan(aligntime) % if trials time-aligned
        if times(1) <= aligntime && times(frames) >= aligntime
            plot([aligntime aligntime],[min(outtrials) max(outtrials)],...
                'k','Linewidth',ZEROWIDTH); % plot vertical line at time 0
            % plot vertical line at aligntime
        end
    else % trials not time-aligned
        if times(1) <= 0 && times(frames) >= 0
            plot([0 0],[min(outtrials) max(outtrials)],...
                'k','Linewidth',ZEROWIDTH); % plot smoothed sortwvar
        end
    end
end

if strcmpi(NoShow, 'no') && ( min(outsort) < timelimits(1) ...
        |max(outsort) > timelimits(2))
    ur_outsort = outsort; % store the pre-adjusted values
    fprintf('Not all sortvar values within time vector limits: \n')
    fprintf('        outliers will be shown at nearest limit.\n');
    i = find(outsort< timelimits(1));
    outsort(i) = timelimits(1);
    i = find(outsort> timelimits(2));
    outsort(i) = timelimits(2);
end

if strcmpi(NoShow, 'no')
    if TIMEX
        if Nosort == YES
            
            l=ylabel(img_ylab);
            if ~isempty(img_ylab) && ~strcmpi(img_ylab,'Trials')
                set(gca,'ytick',img_ytick,'yticklabel',img_ytick_lab);
            end
        else
            if exist('phargs','var')
                l=ylabel('Phase-sorted Trials');
            elseif exist('ampargs','var')
                l=ylabel('Amplitude-sorted Trials');
            elseif exist('valargs','var')
                l=ylabel('Voltage-sorted Trials');
            else
                l=ylabel(img_ylab);
                if ~isempty(img_ylab)
                    set(gca,'ytick',img_ytick);
                end
                if ~strcmpi(img_ylab,'Trials')
                    set(gca,'yticklabel',img_ytick_lab);
                end
            end
        end
    else % if switch x<->y axes
        if Nosort == YES && NoTimeflag==NO
            
            l=xlabel(img_ylab);
            if ~isempty(img_ylab) && ~strcmpi(img_ylab,'Trials')
                set(gca,'xtick',img_ytick,'xticklabel',img_ytick_lab);
            end
        else
            if exist('phargs')
                l=ylabel('Phase-sorted Trials');
            elseif NoTimeflag == NO
                l=xlabel('Sorted Trials');
            else
                l=xlabel(img_ylab);
                if ~isempty(img_ylab) && ~strcmpi(img_ylab,'Trials')
                    set(gca,'xtick',img_ytick,'xticklabel',img_ytick_lab);
                end
            end
        end
    end
    set(l,'FontSize',LABELFONT);
    
    if ~strcmpi(plotmode, 'topo')
        t=title(titl);
        set(t,'FontSize',LABELFONT);
    else
        NAME_OFFSETX = 0.1;
        NAME_OFFSETY = 0.2;
        xx = xlim; xmin = xx(1); xdiff = xx(2)-xx(1); xpos = double(xmin+NAME_OFFSETX*xdiff);
        yy = ylim; ymax = yy(2); ydiff = yy(2)-yy(1); ypos = double(ymax-NAME_OFFSETY*ydiff);
        t=text(xpos, ypos,titl);
        axis off;
    end
    
    set(gca,'Box','off');
    set(gca,'Fontsize',TICKFONT);
    set(gca,'color',BACKCOLOR);
    if Erpflag == NO && NoTimeflag == NO
        if exist('NoTimesPassed')~=1
            l=xlabel('Time (ms)');
        else
            l=xlabel('Frames');
        end
        set(l,'Fontsize',LABELFONT);
    end
end

%
%% %%%%%%%%%%%%%%%%%% Overplot sortvar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if strcmpi(NoShow, 'no')
    
    if NoShowVar == YES
        fprintf('Not overplotting sorted sortvar on data.\n');
        
    elseif isnan(aligntime) % plot sortvar on un-aligned data
        
        if Nosort == NO;
            fprintf('Overplotting sorted sortvar on data.\n');
        end
        hold on;
        if TIMEX      % overplot sortvar
            plot(outsort,outtrials,'k','LineWidth',SORTWIDTH);
        else
            plot(outtrials,outsort,'k','LineWidth',SORTWIDTH);
        end
        drawnow
    else % plot re-aligned zeros on sortvar-aligned data
        if Nosort == NO;
            fprintf('Overplotting sorted sortvar on data.\n');
        end
        hold on;
        if TIMEX      % overplot re-aligned 0 time on image
            plot([aligntime aligntime],[min(outtrials) max(outtrials)],...
                'k','LineWidth',SORTWIDTH);
        else
            plot([[min(outtrials) max(outtrials)],aligntime aligntime],...
                'k','LineWidth',SORTWIDTH);
        end
        fprintf('Overplotting realigned times-zero on data.\n');
        hold on;
        
        if TIMEX      % overplot realigned sortvar on image
            plot(0+aligntime-outsort,outtrials,'k','LineWidth',ZEROWIDTH);
        else
            plot(0+outtrials,aligntime-outsort,'k','LineWidth',ZEROWIDTH);
        end
        drawnow
    end
end

%
%% %%%%%%%%%%%%%%%%%% Overplot auxvar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if strcmpi(NoShow, 'no')
    if ~isempty(auxvar)
        fprintf('Overplotting auxvar(s) on data.\n');
        hold on;
        auxtrials = outtrials(:)' ; % make row vector
        
        if exist('auxcolors')~=1 % If no auxcolors specified
            auxcolors = cell(1,size(auxvar,1));
            for c=1:size(auxvar,1)
                auxcolors(c) = {'k'};       % plot auxvars as black trace(s)
            end
        end
        if length(auxcolors) < size(auxvar,1)
            nauxColors = length(auxcolors);
            for k=nauxColors+1:size(auxvar,1)
                auxcolors = { auxcolors{:} auxcolors{1+rem(k-1,nauxColors)}};
            end
        end
        for c=1:size(auxvar,1)
            auxcolor = auxcolors{c};
            if ~isempty(auxcolor)
                if isnan(aligntime) % plot auxvar on un-aligned data
                    if TIMEX      % overplot auxvar
                        plot(auxvar(c,:)',auxtrials',auxcolor,'LineWidth',SORTWIDTH);
                    else
                        plot(auxtrials',auxvar(c,:)',auxcolor,'LineWidth',SORTWIDTH);
                    end
                    drawnow
                else % plot re-aligned zeros on sortvar-aligned data
                    if TIMEX      % overplot realigned 0-time on image
                        plot(auxvar(c,:)',auxtrials',auxcolor,'LineWidth',ZEROWIDTH);
                    else
                        plot(0+auxtrials',aligntime-auxvar(c,:)',auxcolor,'LineWidth',ZEROWIDTH);
                    end
                    drawnow
                end % aligntime
            end % if auxcolor
        end % c
    end % auxvar
    if exist('outpercent')
        for index = 1:length(outpercent)
            if isnan(aligntime) % plot auxvar on un-aligned data
                plot(outpercent{index},outtrials,'k','LineWidth',SORTWIDTH);
            else
                plot(aligntime-outpercent{index},outtrials,'k','LineWidth',SORTWIDTH);
            end
        end
    end
end
%
%% %%%%%%%%%%%%%%%%%%%%%% Plot colorbar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if strcmpi(NoShow, 'no')
    if Colorbar == YES
        pos=get(ax1,'Position');
        axcb=axes('Position',...
            [pos(1)+pos(3)+0.02 pos(2) ...
            0.03 pos(4)]);
        cbar(axcb,0,[mindat,maxdat]); % plot colorbar to right of image
        title(cbar_title);
        set(axcb,'fontsize',TICKFONT,'xtick',[]);
        % drawnow
        axes(ax1); % reset current axes to the erpimage
    end
end

%
%% %%%%%%%%%%%%%%%%%%%%% Compute ERP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
erp = [];

if Erpflag == YES
    if exist('erpalpha')
        if erp_ptiles>1,
            fprintf(['\nOnly plotting one ERP (i.e., not plotting ERPs from %d percentile splits) ' ...
                'because ''erpalpha'' option was chosen.  You can''t plot both.\n\n'],erp_ptiles)
        end
        [erp erpsig] = nan_mean(fastif(length(tsurdata) > 0, tsurdata',urdata'), ...
            erpalpha);
        fprintf('   Mean ERP (p<%g) significance threshold: +/-%g\n', ...
            erpalpha,mean(erpsig));
    else
        %potentially make ERPs of 50%, 33%, or 25% split of trials
        n_trials=size(urdata,2);
        trials_step=round(n_trials/erp_ptiles);
        erp=zeros(erp_ptiles,size(urdata,1));
        for ploop=1:erp_ptiles,
            ptile_trials=[1:trials_step]+(ploop-1)*trials_step;
            if max(ptile_trials)>n_trials,
                ptile_trials=ptile_trials(1):n_trials;
            end
            if length(tsurdata) > 0
                erp(ploop,:) = nan_mean(tsurdata(:,ptile_trials)');
            else
                erp(ploop,:) = nan_mean(urdata(:,ptile_trials)');
            end
        end
        % else
        %orig line
        %[erp] = nan_mean(fastif(length(tsurdata) > 0, tsurdata', urdata'));
        %end
    end % compute average ERP, ignoring nan's
end

if Erpflag == YES && strcmpi(NoShow, 'no')
    axes(ax1); % reset current axes to the erpimage
    xtick = get(ax1,'Xtick');               % remember x-axis tick locations
    xticklabel = get(ax1,'Xticklabel');     % remember x-axis tick locations
    set(ax1, 'xticklabel', []);
    widthxticklabel = size(xticklabel,2);
    xticklabel = cellstr(xticklabel);
    for tmpindex = 1:length(xticklabel)
        if length(xticklabel{tmpindex}) < widthxticklabel
            spaces = char(ones(1,ceil((widthxticklabel-length(xticklabel{tmpindex}))/2) )*32);
            xticklabel{tmpindex} = [spaces xticklabel{tmpindex}];
        end
    end
    xticklabel = strvcat(xticklabel);
    if Erpstdflag == YES
        stdev = nan_std(urdata');
    end
    %
    %%%%%% Plot ERP time series below image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    if ~isempty(erpsig)
        erpsig = [erpsig;-1*erpsig];
    end
    if isnan(maxerp)
        fac = 10;
        maxerp = 0;
        while maxerp == 0
            maxerp = round(fac*YEXPAND*max(erp))/fac; % minimal decimal places
            fac = 10*fac;
        end
        if Erpstdflag == YES
            fac = fac/10;
            maxerp = max(maxerp, round(fac*YEXPAND*max(erp+stdev))/fac);
        end
        if ~isempty(erpsig)
            maxerp = max(maxerp, round(fac*YEXPAND*max(erpsig))/fac);
        end
        maxerp=max(maxerp);
    end
    if isnan(minerp)
        fac = 1;
        minerp = 0;
        while minerp == 0
            minerp = round(fac*YEXPAND*min(erp))/fac; % minimal decimal places
            fac = 10*fac;
        end
        if Erpstdflag == YES
            fac = fac/10;
            minerp = min(minerp, round(fac*YEXPAND*min(erp-stdev))/fac);
        end
        if ~isempty(erpsig)
            minerp = min(minerp, round(fac*YEXPAND*min(erpsig))/fac);
        end
        minerp=min(minerp);
    end
    limit = [timelimits(1:2) minerp maxerp];
    
    if ~isnan(coherfreq)
        set(ax1,'Xticklabel',[]);     % remove tick labels from bottom of image
        ax2=axes('Position',...
            [gcapos(1) gcapos(2)+2/3*image_loy*gcapos(4) ...
            gcapos(3) (image_loy/3-YGAP)*gcapos(4)]);
    else
        ax2=axes('Position',...
            [gcapos(1) gcapos(2) ...
            gcapos(3) image_loy*gcapos(4)]);
    end
    fprintf('Plotting the ERP trace below the ERP image\n');
    if Erpstdflag == YES
        if Showwin
            tmph = plot1trace(ax2,times,erp,limit,[],stdev,times(winlocs),erp_grid,erp_vltg_ticks); % plot ERP +/-stdev
        else
            tmph = plot1trace(ax2,times,erp,limit, [], stdev,[],erp_grid,erp_vltg_ticks); % plot ERP +/-stdev
        end
    elseif ~isempty('erpsig')
        if Showwin
            tmph = plot1trace(ax2,times,erp,limit,erpsig,[],times(winlocs),erp_grid,erp_vltg_ticks); % plot ERP and 0+/-alpha threshold
        else
            tmph = plot1trace(ax2,times,erp,limit,erpsig,[],[],erp_grid,erp_vltg_ticks); % plot ERP and 0+/-alpha threshold
        end
    else % plot ERP alone - no significance or std dev plotted
        if Showwin
            tmph = plot1trace(ax2,times,erp,limit,[],[],times(winlocs),erp_grid,erp_vltg_ticks); % plot ERP alone
        else
            tmph = plot1trace(ax2,times,erp,limit,[],[],[],erp_grid,erp_vltg_ticks); % plot ERP alone
        end
    end
    
    if ~isnan(aligntime)
        line([aligntime aligntime],[limit(3:4)*1.1],'Color','k','LineWidth',ZEROWIDTH); % x=median sort value
        line([0 0],[limit(3:4)*1.1],'Color','k','LineWidth',ZEROWIDTH); % x=median sort value
        % remove y axis
        if length(tmph) > 1
            delete(tmph(end));
        end
    end
    
    set(ax2,'Xtick',xtick);        % use same Xticks as erpimage above
    if ~isnan(coherfreq)
        set(ax2,'Xticklabel',[]);    % remove tick labels from ERP x-axis
    else % bottom axis
        set(ax2,'Xticklabel',xticklabel); % add ticklabels to ERP x-axis
    end
    
    if isnan(coherfreq)            % if no amp and coher plots below . . .
        if TIMEX && NoTimeflag == NO
            if exist('NoTimesPassed')~=1
                l=xlabel('Time (ms)');
            else
                l=xlabel('Frames');
            end
            set(l,'FontSize',LABELFONT);
            % $$$       else
            % $$$         if exist('NoTimesPassed')~=1
            % $$$           l=ylabel('Time (ms)');
            % $$$         else
            % $$$           l=ylabel('Frames');
            % $$$         end
            % $$$         set(l,'FontSize',LABELFONT);
        end
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot vert lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if ~isempty(verttimes)
        if size(verttimes,1) == ntrials
            vts=sort(verttimes);
            vts = vts(ceil(ntrials/2),:); % plot median verttimes values if a matrix
        else
            vts = verttimes(:)';  % make verttimes a row vector
        end
        for vt = vts
            if isnan(aligntime)
                if TIMEX      % overplot vt on ERP
                    mydotstyle = DOTSTYLE;
                    if exist('auxcolors') && ...
                            length(verttimes) == length(verttimesColors)
                        mydotstyle = verttimesColors{find(verttimes == vt)};
                    end
                    plot([vt vt],[limit(3:4)],mydotstyle,'Linewidth',VERTWIDTH);
                else
                    plot([min(outtrials) max(outtrials)],[limit(3:4)],DOTSTYLE,...
                        'Linewidth',VERTWIDTH);
                end
            else
                if TIMEX      % overplot realigned vt on ERP
                    plot(repmat(median(aligntime+vt-outsort),1,2),[limit(3),limit(4)],...
                        DOTSTYLE,'LineWidth',VERTWIDTH);
                else
                    plot([limit(3),limit(4)],repmat(median(aligntime+vt-outsort),1,2),...
                        DOTSTYLE,'LineWidth',VERTWIDTH);
                end
            end
        end
    end
    
    limit = double(limit);
    ydelta = double(1/10*(limit(2)-limit(1)));
    ytextoffset = double(limit(1)-1.1*ydelta);
    ynumoffset  = double(limit(1)-0.3*ydelta); % double for Matlab 7
    
    %Far left axis max and min labels not needed now that there are tick
    %marks
    %t=text(ynumoffset,0.7*limit(3), num2str(limit(3)));
    %set(t,'HorizontalAlignment','right','FontSize',TICKFONT)
    %t=text(ynumoffset,0.7*limit(4), num2str(limit(4)));
    %set(t,'HorizontalAlignment','right','FontSize',TICKFONT)
    
    ynum = 0.7*(limit(3)+limit(4))/2;
    t=text(ytextoffset,ynum,yerplabel,'Rotation',90);
    set(t,'HorizontalAlignment','center','FontSize',LABELFONT)
    
    if ~exist('YDIR')
        error('\nerpimage(): Default YDIR not read from ''icadefs.m''');
    end
    if YDIR == 1
        set(ax2,'ydir','normal')
    else
        set(ax2,'ydir','reverse')
    end
    
    set(ax2,'Fontsize',TICKFONT);
    set(ax2,'Box','off','color',BACKCOLOR);
    drawnow
end

%
%% %%%%%%%%%%%%%%%%%%% Plot amp, coher time series %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if ~isnan(coherfreq)
    if freq > 0
        coherfreq = mean(freq); % use phase-sort frequency
    end
    %
    %%%%%% Plot amp axis below ERP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    if ~Allampsflag %%%% don't repeat computation if already done for 'allamps'
        
        fprintf('Computing and plotting amplitude at %g Hz.\n',coherfreq);
        
        if ~isnan(signifs) || Cohsigflag==NO % don't compute or plot signif. levels
            [amps,cohers] = phasecoher(urdata,size(times,2),srate,coherfreq,cycles);
            if ~isnan(signifs)
                ampsig = signifs([1 2]);
                fprintf('Using supplied amplitude significance levels: [%g,%g]\n',...
                    ampsig(1),ampsig(2));
                cohsig = signifs(3);
                fprintf('Using supplied coherence significance level: %g\n',cohsig);
            end
        else % compute amps, cohers with significance
            fprintf(...
                'Computing and plotting %g coherence significance level at %g Hz...\n',...
                alpha,                          coherfreq);
            [amps,cohers,cohsig,ampsig] = ...
                phasecoher(urdata,size(times,2),srate,coherfreq,cycles,alpha);
            fprintf('Coherence significance level: %g\n',cohsig);
            ampsig = 20*log10(ampsig); % convert to dB
        end
        amps = 20*log10(amps); % convert to dB
        
        if isnan(baseamp) % if baseamp not specified in 'limits'
            if ~isempty(baselinedb) && isnan(baselinedb(1))
                disp('Not removing amplitude baseline');
            else
                if isempty(baselinedb)
                    baselinedb = [times(1) DEFAULT_BASELINE_END];
                end
                base = find(times >= baselinedb(1) & times<=baselinedb(end)); % use default baseline end point (ms)
                if length(base)<2
                    base = 1:floor(length(times)/4); % default first quarter-epoch
                    fprintf('Using %g to %g ms as amplitude baseline.\n',...
                        times(1),times(base(end)));
                end
                [amps,baseamp] = rmbase(amps,length(times),base); % remove dB baseline
                fprintf('Removed baseline amplitude of %d dB for plotting.\n',baseamp);
            end
        else % if 'basedB' specified in 'limits' (in dB)
            fprintf('Removing specified baseline amplitude of %d dB for plotting.\n',...
                baseamp);
            amps = amps-baseamp; % remove specified dB baseline
        end
        
        fprintf('Data amplitude levels: [%g %g] dB\n',min(amps),max(amps));
        
        if alpha>0 % if computed significance levels
            ampsig = ampsig - baseamp;
            fprintf('Data amplitude significance levels: [%g %g] dB\n',ampsig(1),ampsig(2));
        end
        
    end % ~Allampsflag
    
    if strcmpi(NoShow, 'no')
        axis('off') % rm ERP axes axis and labels
        
        v=axis;
        minampERP=v(3);
        maxampERP=v(4);
        
        if ~exist('ynumoffset')
            limit = [timelimits(1:2) -max(abs([minampERP maxampERP])) max(abs([minampERP maxampERP]))];
            limit = double(limit);
            ydelta = double(1/10*(limit(2)-limit(1)));
            ytextoffset = double(limit(1)-1.1*ydelta);
            ynumoffset  = double(limit(1)-0.3*ydelta); % double for Matlab 7
        end
        
        t=text(double(ynumoffset),double(maxampERP),num2str(maxampERP,3));
        set(t,'HorizontalAlignment','right','FontSize',TICKFONT);
        
        t=text(double(ynumoffset),double(minampERP), num2str(minampERP,3));
        set(t,'HorizontalAlignment','right','FontSize',TICKFONT);
        
        ax3=axes('Position',...
            [gcapos(1) gcapos(2)+1/3*image_loy*gcapos(4) ...
            gcapos(3) (image_loy/3-YGAP)*gcapos(4)]);
        
        if isnan(maxamp) % if not specified
            fac = 1;
            maxamp = 0;
            while maxamp == 0
                maxamp = floor(YEXPAND*fac*max(amps))/fac; % minimal decimal place
                fac = 10*fac;
            end
            maxamp = maxamp + 10/fac;
            if Cohsigflag
                if ampsig(2)>maxamp
                    if ampsig(2)>0
                        maxamp = 1.01*(ampsig(2));
                    else
                        maxamp = 0.99*(ampsig(2));
                    end
                end
            end
        end
        if isnan(maxamp), maxamp = 0; end   % In case the above iteration went on
        % until fac = Inf and maxamp = NaN again.
        
        if isnan(minamp) % if not specified
            fac = 1;
            minamp = 0;
            while minamp == 0
                minamp = floor(YEXPAND*fac*max(-amps))/fac; % minimal decimal place
                fac = 10*fac;
            end
            minamp = minamp + 10/fac;
            minamp = -minamp;
            if Cohsigflag
                if ampsig(1)< minamp
                    if ampsig(1)<0
                        minamp = 1.01*(ampsig(1));
                    else
                        minamp = 0.99*(ampsig(1));
                    end
                end
            end
        end
        if isnan(minamp), minamp = 0; end   % In case the above iteration went on
        % until fac = Inf and minamp = NaN again.
        
        fprintf('Plotting the ERSP amplitude trace below the ERP\n');
        fprintf('Min, max plotting amplitudes: [%g, %g] dB\n',minamp,maxamp);
        fprintf('     relative to baseamp: %g dB\n',baseamp);
        if Cohsigflag
            ampsiglims = [repmat(ampsig(1)-mean(ampsig),1,length(times))];
            ampsiglims = [ampsiglims;-1*ampsiglims];
            plot1trace(ax3,times,amps,[timelimits minamp(1) maxamp(1)],ampsiglims,[],[],0); % plot AMP
        else
            plot1trace(ax3,times,amps,[timelimits minamp(1) maxamp(1)],[],[],[],0); % plot AMP
        end
        
        if ~isnan(aligntime)
            line([aligntime aligntime],[minamp(1) maxamp(1)]*1.1,'Color','k');
            % x=median sort value
        end
        if exist('xtick') % Added -JH
            set(ax3,'Xtick',xtick);
        end
        set(ax3,'Xticklabel',[]);   % remove tick labels from bottom of image
        set(ax3,'Yticklabel',[]);   % remove tick labels from left of image
        set(ax3,'YColor',BACKCOLOR);
        axis('off');
        set(ax3,'Box','off','color',BACKCOLOR);
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot vert marks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        if ~isempty(verttimes)
            if size(verttimes,1) == ntrials
                vts=sort(verttimes);
                vts = vts(ceil(ntrials/2),:); % plot median values if a matrix
            else
                vts=verttimes(:)';
            end
            for vt = vts
                if isnan(aligntime)
                    if TIMEX      % overplot vt on amp
                        mydotstyle = DOTSTYLE;
                        if exist('auxcolors') && ...
                                length(verttimes) == length(verttimesColors)
                            mydotstyle = verttimesColors{find(verttimes == vt)};
                        end
                        
                        plot([vt vt],[minamp(1) maxamp(1)],mydotstyle,...
                            'Linewidth',VERTWIDTH);
                    else
                        plot([min(outtrials) max(outtrials)],[minamp(1) maxamp(1)], ...
                            DOTSTYLE,...
                            'Linewidth',VERTWIDTH);
                    end
                else
                    if TIMEX      % overplot realigned vt on amp
                        plot(repmat(median(aligntime+vt-outsort),1,2), ...
                            [minamp(1),maxamp(1)],DOTSTYLE,...
                            'LineWidth',VERTWIDTH);
                    else
                        plot([minamp,maxamp],repmat(median(aligntime+vt-outsort),1,2), ...
                            DOTSTYLE,...
                            'LineWidth',VERTWIDTH);
                    end
                end
            end
        end
        
        if 0 % Cohsigflag % plot amplitude significance levels
            hold on
            plot([timelimits(1) timelimits(2)],[ampsig(1) ampsig(1)] - mean(ampsig),'r',...
                'linewidth',SIGNIFWIDTH);
            plot([timelimits(1) timelimits(2)],[ampsig(2) ampsig(2)] - mean(ampsig),'r',...
                'linewidth',SIGNIFWIDTH);
        end
        
        if ~exist('ynumoffset')
            limit = [timelimits(1:2) -max(abs([minamp maxamp])) max(abs([minamp maxamp]))];
            limit = double(limit);
            ydelta = double(1/10*(limit(2)-limit(1)));
            ytextoffset = double(limit(1)-1.1*ydelta);
            ynumoffset  = double(limit(1)-0.3*ydelta); % double for Matlab 7
        end
        
        t=text(ynumoffset,maxamp, num2str(maxamp,3));
        set(t,'HorizontalAlignment','right','FontSize',TICKFONT);
        
        t=text(ynumoffset,0, num2str(0));
        set(t,'HorizontalAlignment','right','FontSize',TICKFONT);
        
        t=text(ynumoffset,minamp, num2str(minamp,3));
        set(t,'HorizontalAlignment','right','FontSize',TICKFONT);
        
        t=text(ytextoffset,(maxamp+minamp)/2,'ERSP','Rotation',90);
        set(t,'HorizontalAlignment','center','FontSize',LABELFONT);
        
        axtmp = axis;
        dbtxt= text(1/13*(axtmp(2)-axtmp(1))+axtmp(1), ...
            11/13*(axtmp(4)-axtmp(3))+axtmp(3), ...
            [num2str(baseamp,4) ' dB']);
        set(dbtxt,'fontsize',TICKFONT);
        drawnow;
        set(ax3, 'xlim', timelimits);
        set(ax3, 'ylim', [minamp(1) maxamp(1)]);
        
        %
        %%%%%% Make coher axis below amp %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %
        ax4=axes('Position',...
            [gcapos(1) gcapos(2) ...
            gcapos(3) (image_loy/3-YGAP)*gcapos(4)]);
        if isnan(maxcoh)
            fac = 1;
            maxcoh = 0;
            while maxcoh == 0
                maxcoh = floor(YEXPAND*fac*max(cohers))/fac; % minimal decimal place
                fac = 10*fac;
            end
            maxcoh = maxcoh + 10/fac;
            if maxcoh>1
                maxcoh=1; % absolute limit
            end
        end
        if isnan(mincoh)
            mincoh = 0;
        end
        fprintf('Plotting the ITC trace below the ERSP\n');
        if Cohsigflag % plot coherence significance level
            cohsiglims = [repmat(cohsig,1,length(times));zeros(1,length(times))];
            coh_handle = plot1trace(ax4,times,cohers,[timelimits mincoh maxcoh],cohsiglims,[],[],0);
            % plot COHER, fill sorting window
        else
            coh_handle = plot1trace(ax4,times,cohers,[timelimits mincoh maxcoh],[],[],[],0); % plot COHER
        end
        if ~isnan(aligntime)
            line([aligntime aligntime],[[mincoh maxcoh]*1.1],'Color','k');
            % x=median sort value
        end
        % set(ax4,'Xticklabel',[]);    % remove tick labels from bottom of
        % image
        if exist('xtick')
            set(ax4,'Xtick',xtick);
            set(ax4,'Xticklabel',xticklabel);
        end
        set(ax4,'Ytick',[]);
        set(ax4,'Yticklabel',[]);      % remove tick labels from left of image
        set(ax4,'YColor',BACKCOLOR);
        
        if ~isempty(verttimes)
            if size(verttimes,1) == ntrials
                vts=sort(verttimes);
                vts = vts(ceil(ntrials/2),:); % plot median values if a matrix
            else
                vts = verttimes(:)';  % make verttimes a row vector
            end
            for vt = vts
                if isnan(aligntime)
                    if TIMEX      % overplot vt on coher
                        mydotstyle = DOTSTYLE;
                        if exist('auxcolors') && ...
                                length(verttimes) == length(verttimesColors)
                            mydotstyle = verttimesColors{find(verttimes == vt)};
                        end
                        
                        plot([vt vt],[mincoh maxcoh],mydotstyle,'Linewidth',VERTWIDTH);
                    else
                        plot([min(outtrials) max(outtrials)],...
                            [mincoh maxcoh],DOTSTYLE,'Linewidth',VERTWIDTH);
                    end
                else
                    if TIMEX      % overplot realigned vt on coher
                        plot(repmat(median(aligntime+vt-outsort),1,2),...
                            [mincoh,maxcoh],DOTSTYLE,'LineWidth',VERTWIDTH);
                    else
                        plot([mincoh,maxcoh],repmat(median(aligntime+vt-outsort),1,2),...
                            DOTSTYLE,'LineWidth',VERTWIDTH);
                    end
                end
            end
        end
        
        t=text(ynumoffset,0, num2str(0));
        set(t,'HorizontalAlignment','right','FontSize',TICKFONT);
        
        t=text(ynumoffset,maxcoh, num2str(maxcoh));
        set(t,'HorizontalAlignment','right','FontSize',TICKFONT);
        
        t=text(ytextoffset,maxcoh/2,'ITC','Rotation',90);
        set(t,'HorizontalAlignment','center','FontSize',LABELFONT);
        drawnow
        
        %if Cohsigflag % plot coherence significance level
        %hold on
        %plot([timelimits(1) timelimits(2)],[cohsig cohsig],'r',...
        %'linewidth',SIGNIFWIDTH);
        %end
        
        set(ax4,'Box','off','color',BACKCOLOR);
        set(ax4,'Fontsize',TICKFONT);
        if NoTimeflag==NO
            if exist('NoTimesPassed')~=1
                l=xlabel('Time (ms)');
            else
                l=xlabel('Frames');
            end
            set(l,'Fontsize',LABELFONT);
        end
        axtmp = axis;
        hztxt=text(10/13*(axtmp(2)-axtmp(1))+axtmp(1), ...
            8/13*(axtmp(4)-axtmp(3))+axtmp(3), ...
            [num2str(coherfreq,4) ' Hz']);
        set(hztxt,'fontsize',TICKFONT);
    end;% NoShow
else
    amps   = [];    % null outputs unless coherfreq specified
    cohers = [];
end

if VERS >= 8.04
    axhndls = {ax1 axcb ax2 ax3 ax4};
else
    axhndls = [ax1 axcb ax2 ax3 ax4];
end

if exist('ur_outsort')
    outsort = ur_outsort; % restore outsort clipped values, if any
end
if nargout<1
    data = []; % don't spew out data if no args out and no ;
end

%
%% %%%%%%%%%%%%% Plot a TOPOPLOT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if (~isempty(topomap)) && strcmpi(NoShow, 'no')
    h(12)=axes('Position',...
        [gcapos(1)+0.10*gcapos(3) gcapos(2)+0.92*gcapos(4),...
        0.20*gcapos(3) 0.14*gcapos(4)]);
    % h(12) = subplot('Position',[.10 .86 .20 .14]);
    fprintf('Plotting a topo map in upper left.\n');
    eloc_info.plotrad = [];
    if length(topomap) == 1
        try
            topoplot(topomap,eloc_file,'electrodes','off', ...
                'style', 'blank', 'emarkersize1chan', 10, 'chaninfo', eloc_info);
        catch
            fprintf('topoplot() plotting failed.\n');
        end
    else
        try
            topoplot(topomap,eloc_file,'electrodes','off', 'chaninfo', eloc_info);
        catch
            fprintf('topoplot() plotting failed.\n');
        end
    end
    axis('square')
    if VERS >= 8.04
        axhndls = {axhndls h(12)};
    else
        axhndls = [axhndls h(12)];
    end
end

%
%% %%%%%%%%%%%%% Plot a spectrum %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
SPECFONT = 10;
if (~isempty(lospecHz)) && strcmpi(NoShow, 'no')
    h(13)=axes('Position',...
        [gcapos(1)+0.82*gcapos(3) ...
        gcapos(2)+0.96*gcapos(4),...
        0.15*gcapos(3)*(0.8/gcapos(3))^0.5 ...
        0.10*gcapos(4)*(0.8/gcapos(4))^0.5]);
    
    % h(13) = subplot('Position',[.75 .88 .15 .10]);
    fprintf('Plotting the data spectrum in upper right.\n');
    winlength = frames;
    if winlength > 512
        for k=2:5
            if rem(winlength,k) == 0
                break
            end
        end
        winlength = winlength/k;
    end
    
    % [Pxx, Pxxc, F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP,P)
    if exist('pwelch') == 2
        [Pxx,F] = pwelch(reshape(urdata,1,size(urdata,1)*size(urdata,2)),...
            frames,0,max(1024,pow2(ceil(log2(frames)))),srate);
        % [Pxx,F] = psd(reshape(urdata,1,size(urdata,1)*size(urdata,2)),512,srate,winlength,0,0.05);
    else
        [Pxx,F] = spec(reshape(urdata,1,size(urdata,1)*size(urdata,2)),...
            max(1024,pow2(ceil(log2(frames)))),srate,frames,0);
        % [Pxx,F] = spec(reshape(urdata,1,size(urdata,1)*size(urdata,2)),512,srate,winlength,0);
    end
    
    plot(F,10*log10(Pxx));
    goodfs = find(F>= lospecHz & F <= hispecHz);
    maxgfs = max(10*log10(Pxx(goodfs)));
    mingfs = min(10*log10(Pxx(goodfs)));
    axis('square')
    axis([lospecHz hispecHz mingfs-1 maxgfs+1]);
    set(h(13),'Box','off','color',BACKCOLOR);
    set(h(13),'Fontsize',SPECFONT);
    set(h(13),'Xscale',SpecAxis); % added 'log' or 'linear' freq axis scaling -SM 5/31/12
    l=ylabel('dB');
    set(l,'Fontsize',SPECFONT);
    if ~isnan(coherfreq)
        hold on; plot([coherfreq,coherfreq],[mingfs maxgfs],'r');
    end
    if VERS >= 8.04
        axhndls = {axhndls h(13)};
    else
        axhndls = [axhndls h(13)];
    end
end

%
%% %%%%%%%%%%%%% save plotting limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
limits = [min(times) max(times) minerp maxerp minamp maxamp mincoh maxcoh];
limits = [limits baseamp coherfreq];  % add coherfreq to output limits array

%
%% %%%%%%%%%%%%% turn on AXCOPY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if strcmpi(NoShow, 'no')
    axcopy(gcf);
    % eegstr = 'img=get(gca,''''children''''); if (strcmp(img(end),''''type''''),''''image''''), img=get(img(end),''''CData''''); times=get(img(end),''''Xdata''''); clf; args = [''''limits'''' '''','''' times(1) '''','''' times(end)]; if exist(''''EEG'''')==1, args = [args '''','''' ''''srate'''' '''','''' EEG.srate]; end eegplot(img,args); end';
    % axcopy(gcf,eegstr);
end


% returning outsort
if exist('outpercent')
    outsort = { outsort outpercent };
end

fprintf('Done.\n\n');

%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%  End ERPIMAGE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if strcmpi(NoShow, 'no')
    axes('position',gcapos);
    axis off
end
warning on;

return
%
%% %%%%%%%%%%%%%%%%% function PLOT1TRACE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
function [plot_handle] = plot1trace(ax,times,trace,axlimits,signif,stdev,winlocs,erp_grid,erp_vltg_ticks)
%function [plot_handle] = plot1trace(ax,times,trace,axlimits,signif,stdev,winlocs,erp_grid,erp_vltg_ticks)
%                           If signif is [], plot trace +/- stdev
%                           Else if signif, plot trace and signif(1,:)&signif(2,:) fill.
%                           Else, plot trace alone.
%                           If winlocs not [], plot grey back image(s) in sort window
%                                       winlocs(1,1)-> winlocs(1,end) (ms)
%                                        ...
%                                       winlocs(end,1)-> winlocs(end,end) (ms)

FILLCOLOR    = [.66 .76 1];
WINFILLCOLOR    = [.88 .92 1];
ERPDATAWIDTH = 2;
ERPZEROWIDTH = 2;
axes(ax);

if nargin<9,
    erp_vltg_ticks=[]; %if erp_vltg_ticks is not empty, those will be the ticks used on the y-axis (voltage axis for ERPs)
end

if ~isempty(winlocs)
    for k=1:size(winlocs,1)
        winloc = winlocs(k,:);
        fillwinx = [winloc winloc(end:-1:1)];
        hannwin = makehanning(length(winloc));
        hannwin = hannwin./max(hannwin); % make max = 1
        hannwin = hannwin(:)'; % make row vector
        if ~isempty(axlimits) && sum(isnan(axlimits))==0
            % fillwiny = [repmat(axlimits(3),1,length(winloc)) repmat(axlimits(4),1,length(winloc))];
            fillwiny = [hannwin*axlimits(3) hannwin*axlimits(4)];
        else
            % fillwiny = [repmat(min(trace)*1.1,1,length(winloc)) repmat(max(trace)*1.1,1,length(winloc))];
            fillwiny = [hannwin*2*min(trace) hannwin*2*max(trace)];
        end
        fillwh = fill(fillwinx,fillwiny, WINFILLCOLOR); hold on    % plot 0+alpha
        set(fillwh,'edgecolor',WINFILLCOLOR-[.00 .00 0]); % make edges NOT highlighted
    end
end
if ~isempty(signif);% (2,times) array giving upper and lower signif limits
    filltimes = [times times(end:-1:1)];
    if size(signif,1) ~=2 || size(signif,2) ~= length(times)
        fprintf('plot1trace(): signif array must be size (2,frames)\n')
        return
    end
    fillsignif = [signif(1,:) signif(2,end:-1:1)];
    fillh = fill(filltimes,fillsignif, FILLCOLOR); hold on    % plot 0+alpha
    set(fillh,'edgecolor',FILLCOLOR-[.02 .02 0]); % make edges slightly highlighted
    % [plot_handle] = plot(times,signif, 'r','LineWidth',1); hold on    % plot 0+alpha
    % [plot_handle] = plot(times,-1*signif, 'r','LineWidth',1); hold on % plot 0-alpha
end
if ~isempty(stdev)
    [st1] = plot(times,trace+stdev, 'r--','LineWidth',1); hold on % plot trace+stdev
    [st2] = plot(times,trace-stdev, 'r--','LineWidth',1); hold on % plot trace-stdev
end
%linestyles={'r','m','c','b'};
% 'LineStyle',linestyles{traceloop}); hold on
linecolor={[0 0 1],[.25 0 .75],[.75 0 .25],[1 0 0]};
plot_handle=zeros(1,size(trace,1));
for traceloop=1:size(trace,1),
    [plot_handle(traceloop)] = plot(times,trace(traceloop,:),'LineWidth',ERPDATAWIDTH, ...
        'color',linecolor{traceloop}); hold on
end

%Assume that multiple traces are equally sized divisions of data
switch size(trace,1),
    case 2
        legend('Lower 50%','Higher 50%');
    case 3
        legend('Lowest 33%','Middle 33%','Highest 33%');
    case 4
        legend('Lowest 25%','2nd Lowest 25%','3rd Lowest 25%','Highest 25%');
end

if ~isempty(axlimits) && sum(isnan(axlimits))==0
    if axlimits(2)>axlimits(1) && axlimits(4)>axlimits(3)
        axis([axlimits(1:2) 1.1*axlimits(3:4)])
    end
    l1=line([axlimits(1:2)],[0 0],    'Color','k',...
        'linewidth',ERPZEROWIDTH); % y=zero-line
    timebar=0;
    l2=line([1 1]*timebar,axlimits(3:4)*1.1,'Color','k',...
        'linewidth',ERPZEROWIDTH); % x=zero-line
    
    %y-ticks
    if isempty(erp_vltg_ticks),
        shrunk_ylimits=axlimits(3:4)*.8;
        ystep=(shrunk_ylimits(2)-shrunk_ylimits(1))/4;
        if ystep>1,
            ystep=round(ystep);
        else
            ord=orderofmag(ystep);
            ystep=round(ystep/ord)*ord;
        end
        if (sign(shrunk_ylimits(2))*sign(shrunk_ylimits(1)))==1, %y shrunk_ylimits don't include 0
            erp_yticks=shrunk_ylimits(1):ystep:shrunk_ylimits(2);
        else %y limits include 0
            % erp_yticks=0:ystep:shrunk_ylimits(2); %ensures 0 is a tick point
            % erp_yticks=[erp_yticks [0:-ystep:shrunk_ylimits(1)]];
            erp_yticks=0:ystep:axlimits(4); %ensures 0 is a tick point
            erp_yticks=[erp_yticks [0:-ystep:axlimits(3)]];
        end
        if erp_grid,
            set(gca,'ytick',unique(erp_yticks),'ygrid','on');
        else
            set(gca,'ytick',unique(erp_yticks));
        end
    else
        set(gca,'ytick',erp_vltg_ticks);
    end
end

%make ERP traces on top
kids=get(gca,'children')';
for hndl_loop=plot_handle,
    id=find(kids==hndl_loop);
    kids(id)=[];
    kids=[hndl_loop kids];
end
set(gca,'children',kids');
plot_handle=[plot_handle l1 l2];

%
%% %%%%%%%%%%%%%%%%% function PHASEDET %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PHASEDET - function used in erpimage.m
%              Constructs a complex filter at frequency freq
%
function [ang,amp,win] = phasedet(data,frames,srate,nwin,freq)
% Typical values:
%   frames = 768;
%   srate = 256; % Hz
%   nwin = [200:300];
%   freq = 10; % Hz

data = reshape(data,[frames prod(size(data))/frames]);
% number of cycles depends on window size
% number of cycles automatically reduced if smaller window
% Note: as the number of cycles changes, the frequency shifts
%       a little -- this should be fixed

win = exp(2i*pi*freq(:)*[1:length(nwin)]/srate);
win = win .* repmat(makehanning(length(nwin))',length(freq),1);

%tmp =gcf; figure; plot(real(win)); figure(tmp);
%fprintf('ANY NAN ************************* %d\n', any(any(isnan( data(nwin,:)))));

tmpdata = data(nwin,:) - repmat(mean(data(nwin,:), 1), [size(data,1) 1]);
resp = win * tmpdata;
ang = angle(resp);
amp = abs(resp);

%
%% %%%%%%%%%%%% function PRCTLE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
function prctl = prctle(data,pc); % return percentile of a distribution
[prows pcols] = size(pc);
if prows ~= 1 && pcols ~= 1
    error('\nerpimage(): pc must be a scalar or a vector.');
end
if any(pc > 100) || any(pc < 0)
    error('\nerpimage(): pc must be between 0 and 100');
end
[i,j] = size(data);
sortdata = sort(data);
if i==1 || j==1 % if data is a vector
    i = max(i,j); j = 1;
    if i == 1,
        fprintf('  prctle() note: input data is a single scalar!\n')
        y = data*ones(length(pc),1); % if data is scalar, return it
        return;
    end
    sortdata = sortdata(:);
end
pt = [0 100*((1:i)-0.5)./i 100];
sortdata = [min(data); sortdata; max(data)];
prctl = interp1(pt,sortdata,pc);
%
%% %%%%%%%%%%%%%%%%%%%%% function NAN_MEAN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% NAN_MEAN - Take the column means of a matrix, ignoring NaN values.
%              Return significance bounds if alpha (0 < alpha< <1) is given.
%
function [out, outalpha]  = nan_mean(in,alpha)
NPERM = 500;

intrials = size(in,1);
inframes = size(in,2);
nans = find(isnan(in));
in(nans) = 0;
sums = sum(in);
nonnans = ones(size(in));
nonnans(nans) = 0;
nonnans = sum(nonnans);
nononnans = find(nonnans==0);
nonnans(nononnans) = 1;
out = sum(in)./nonnans;
outalpha = [];
if nargin>1
    if NPERM < round(3/alpha)
        NPERM = round(3/alpha);
    end
    fprintf('Performing a permuration test using %d permutations to determine ERP significance thresholds... ',NPERM);
    permerps = zeros(NPERM,inframes);
    for n=1:NPERM
        signs = sign(randn(1,intrials)'-0.5);
        permerps(n,:) = sum(repmat(signs,1,inframes).*in)./nonnans;
        if ~rem(n,50)
            fprintf('%d ',n);
        end
    end
    fprintf('\n');
    permerps = sort(abs(permerps));
    alpha_bnd = floor(2*alpha*NPERM); % two-sided probability threshold
    outalpha = permerps(end-alpha_bnd,:);
end
out(nononnans) = NaN;
%
%% %%%%%%%%%%%%%%%%%%%%% function NAN_STD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
function out = nan_std(in)

nans = find(isnan(in));
in(nans) = 0;

nonnans = ones(size(in));
nonnans(nans) = 0;
nonnans = sum(nonnans);
nononnans = find(nonnans==0);
nonnans(nononnans) = 1;

out = sqrt((sum(in.^2)-sum(in).^2./nonnans)./(nonnans-1));
out(nononnans) = NaN;

% symmetric hanning function
function w = makehanning(n)
if ~rem(n,2)
    w = 0.5*(1 - cos(2*pi*(1:n/2)'/(n+1)));
    w = [w; w(end:-1:1)];
else
    w = 0.5*(1 - cos(2*pi*(1:(n+1)/2)'/(n+1)));
    w = [w; w(end-1:-1:1)];
end
%
%% %%%%%%%%%%%%%%%%%%%%% function COMPUTE_PERCENTILE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
function outpercent = compute_percentile(sortvar, percent, outtrials, winsize);
ntrials = length(sortvar);
outtrials=round(outtrials);
sortvar = [ sortvar sortvar sortvar ];
winvals = [round(-winsize/2):round(winsize/2)];
outpercent = zeros(size(outtrials));
for index = 1:length(outtrials)
    sortvarval = sortvar(outtrials(index)+ntrials+winvals);
    sortvarval = sort(sortvarval);
    outpercent(index) = sortvarval(round((length(winvals)-1)*percent)+1);
end

%
%% %%%%%%%%%%%%%%%%%%%%% function ORDEROFMAG %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
function ord=orderofmag(val)
%function ord=orderofmag(val)
%
% Returns the order of magnitude of the value of 'val' in multiples of 10
% (e.g., 10^-1, 10^0, 10^1, 10^2, etc ...)
% used for computing erpimage trial axis tick labels as an alternative for
% plotting sorting variable

val=abs(val);
if val>=1
    ord=1;
    val=floor(val/10);
    while val>=1,
        ord=ord*10;
        val=floor(val/10);
    end
    return;
else
    ord=1/10;
    val=val*10;
    while val<1,
        ord=ord/10;
        val=val*10;
    end
    return;
end

%
%% %%%%%%%%%%%%%%%%%%%%% function FIND_CRSPND_PT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
function y_pt=find_crspnd_pt(targ,vals,outtrials)
%function id=find_crspnd_pt(targ,vals,outtrials)
%
% Inputs:
%   targ      - desired value of sorting variable
%   vals      - a vector of observed sorting variables (possibly smoothed)
%   outtrials - a vector of y-axis values corresponding to trials in the
%               ERPimage (this will just be 1:n_trials if there's no
%               smoothing)
%
% Output:
%   y_pt  - y-axis value (in units of trials) corresponding to "targ".
%          If "targ" matches more than one y-axis pt, the median point is
%          returned.  If "targ" falls between two points, y_pt is linearly
%          interpolated.
%
% Note: targ and vals should be in the same units (e.g., milliseconds)

%find closest point above
abv=find(vals>=targ);
if isempty(abv),
    %point lies outside of vals range, can't interpolate
    y_pt=[];
    return
end
abv=abv(1);

%find closest point below
blw=find(vals<=targ);
if isempty(blw),
    %point lies outside of vals range, can't interpolate
    y_pt=[];
    return
end
blw=blw(end);

if (vals(abv)==vals(blw)),
    %exact match
    ids=find(vals==targ);
    y_pt=median(outtrials(ids));
else
    %interpolate point
    
    %lst squares inear regression
    B=regress([outtrials(abv) outtrials(blw)]',[ones(2,1) [vals(abv) vals(blw)]']);
    
    %predict outtrial point from target value
    y_pt=[1 targ]*B;
    
end