erpimage

 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 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' (below) > (default) 3 for better frequency specificity,
              {default: plot potential, not amplitudes}

 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.
   '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}
 '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: erpimages(), phasecoher(), rmbase(), cbar(), movav()

0001 % erpimage() - Plot a colored image of a collection of single-trial data epochs, optionally
0002 %              sorted on and/or aligned to an input sorting variable and smoothed across
0003 %              trials with a Gaussian weighted moving-average. (To return event-aligned data
0004 %              without plotting, use eegalign()).  Optionally sort trials on value, amplitude
0005 %              or phase within a specified latency window. Optionally plot the ERP mean
0006 %              and std. dev.and moving-window spectral amplitude and inter-trial coherence
0007 %              at aselected or peak frequency. Optionally 'time warp' the single trial
0008 %              time-domain (potential) or power data to align the plotted data to a series
0009 %              of events with varying latencies that occur in each trial. Click on
0010 %              individual figures parts to examine them separately and zoom (using axcopy()).
0011 % Usage:
0012 %            >> figure; erpimage(data,[],times); % image trials in input order
0013 %
0014 %            >> figure; [outdata,outvar,outtrials,limits,axhndls, ...
0015 %                        erp,amps,cohers,cohsig,ampsig,outamps,...
0016 %                        phsangls,phsamp,sortidx,erpsig] ...
0017 %                            = erpimage(data,sortvar,times,'title',avewidth,decimate,...
0018 %                                  'key', 'val', ...); % use options
0019 % Required input:
0020 %   data     = [vector or matrix] Single-channel input data to image.
0021 %               Formats (1,frames*trials) or (frames,trials)
0022 %
0023 % Optional ordered inputs {with defaults}:
0024 %
0025 %   sortvar  = [vector | []] Variable to sort epochs on (length(sortvar) = nepochs)
0026 %              Example: sortvar may by subject response time in each epoch (in ms)
0027 %              {default|[]: plot in input order}
0028 %   times    = [vector | []] vector of latencies (in ms) for each epoch time point.
0029 %               Else [startms ntimes srate] = [start latency (ms), time points
0030 %               (=frames) per epoch, sampling rate (Hz)]. Else [] -> 0:nframes-1
0031 %               {default: []}
0032 %  'title'   = ['string'] Plot title {default: none}
0033 %   avewidth = [positive scalar (may be non-integer)]. If avg_type is set to 'boxcar'
0034 %               (the default), this is the number of trials used to smooth
0035 %               (vertically) with a moving-average. If avg_type is set to
0036 %               'Gaussian,' this is the standard deviation (in units of
0037 %               trials) of the Gaussian window used to smooth (vertically)
0038 %               with a moving-average.  Gaussian window extends three
0039 %               standard deviations below and three standard deviations above window
0040 %               center (trials beyond window are not incorporated into average). {default: no
0041 %               smoothing}
0042 %   decimate = Factor to decimate|interpolate ntrials by (may be non-integer)
0043 %               Else, if this is large (> sqrt(ntrials)), output this many epochs.
0044 %               {default|0->1}
0045 %
0046 % Optional unordered 'keyword',argument pairs:
0047 %
0048 % Re-align data epochs:
0049 %   'align'  = [latency] Time-lock data to sortvar. Plot sortvar at given latency
0050 %               (in ms). Else Inf -> plot sortvar at median sortvar latency
0051 %               {default: do not align}
0052 %   'timewarp' = {[events], [warpms], {colors}} Time warp ERP, amplitude and phase
0053 %               time-courses before smoothing. 'events' is a matrix whose columns
0054 %               specify the latencies (in ms) at which a series of successive events occur
0055 %               in each trial. 'warpms' is an optional vector of latencies (in ms) to which
0056 %               the series of events should be time locked. (Note: Epoch start and end
0057 %               should not be declared as events or warpms}. If 'warpms' is absent or [],
0058 %               the median of each 'events' column will be used. {colors} contains a
0059 %               list of Matlab linestyles to use for vertical lines marking the occurrence
0060 %               of the time warped events. If '', no line will be drawn for this event
0061 %               column. If fewer colors than event columns, cycles through the given color
0062 %               labels.  Note: Not compatible with 'vert' (below).
0063 %   'renorm' = ['yes'|'no'| formula] Normalize sorting variable to epoch
0064 %               latency range and plot. 'yes'= autoscale. formula must be a linear
0065 %               transformation in the format 'a*x+b'
0066 %               Example of formula: '3*x+2'. {default: 'no'}
0067 %               If sorting by string values like event type, suggested formulas for:
0068 %                 letter string: '1000*x', number string: '30000*x-1500'
0069 %   'noplot' = ['on'|'off'] Do not plot sortvar {default: Plot sortvar if in times range}
0070 %   'noshow' = ['on'|'off'] Do not plot erpimage, simply return outputs {default: 'off'}
0071 %
0072 % Sort data epochs:
0073 % 'nosort'       = ['on'|'off'] Do not sort data epochs. {default: Sort data epochs by
0074 %                  sortvar (see sortvar input above)}
0075 % 'replace_ties' = ['yes'|'no'] Replace trials with the same value of
0076 %                  sortvar with the mean of those trials.  Only works if sorting trials
0077 %                  by sortvar. {default: 'no'}
0078 % 'valsort'      = [startms endms direction] Sort data on (mean) value
0079 %                  between startms and (optional) endms. Direction is 1 or -1.
0080 %                  If -1, plot max-value epoch at bottom {default: sort on sortvar}
0081 % 'phasesort'    = [ms_center prct freq maxfreq topphase] Sort epochs by phase in
0082 %                  a 3-cycle window centered at latency ms_center (ms).
0083 %                  Percentile (prct) in range [0,100] gives percent of trials
0084 %                  to reject for (too) low amplitude. Else, if in range [-100,0],
0085 %                  percent of trials to reject for (too) high amplitude;
0086 %                  freq (Hz) is the phase-sorting frequency. With optional
0087 %                  maxfreq, sort by phase at freq of max power in the data in
0088 %                  the range [freq,maxfreq] (Note: 'phasesort' arg freq overrides
0089 %                  the frequency specified in 'coher'). With optional topphase,
0090 %                  sort by phase, putting topphase (degrees, in range [-180,180])
0091 %                  at the top of the image. Note: 'phasesort' now uses circular
0092 %                  smoothing. Use 'cycles' (below) for wavelet length.
0093 %                  {default: [0 25 8 13 180]}
0094 %  'ampsort'     = [center_ms prcnt freq maxfreq]  Sort epochs by amplitude.
0095 %                  (See 'phasesort' above). If ms_center is 'Inf', then sorting
0096 %                  is by mean power across the time window specified by 'sortwin'
0097 %                  below. If third arg, freq, is < 0, sort by mean power in the range
0098 %                  [ abs(freq)   maxfreq ].
0099 %  'sortwin'     = [start_ms end_ms] If center_ms == Inf in 'ampsort' arg (above), sorts
0100 %                  by mean amplitude across window centers shifted from start_ms
0101 %                  to end_ms by 10 ms.
0102 %  'showwin'     = ['on'|'off'] Show sorting window behind ERP trace. {default: 'off'}
0103 %
0104 % Plot time-varying spectral amplitude instead of potential:
0105 % 'plotamps' = ['on'|'off'] Image amplitudes at each trial and latency instead of potential
0106 %              values. Note: Currently requires 'coher' (below) with alpha signif.
0107 %              Use 'cycles' (below) > (default) 3 for better frequency specificity,
0108 %              {default: plot potential, not amplitudes}
0109 %
0110 % Specify plot parameters:
0111 %   'limits'         = [lotime hitime minerp maxerp lodB hidB locoher hicoher basedB]
0112 %                      Plot axes limits. Can use NaN (or nan, but not Nan) for missing items
0113 %                      and omit late items. Use last input, basedB, to set the
0114 %                      baseline dB amplitude in 'plotamps' plots {default: from data}
0115 %   'sortvar_limits' = [min max] minimum and maximum sorting variable
0116 %                      values to image. This only affects visualization of
0117 %                      ERPimage and ERPs (not smoothing).  Cannot be used
0118 %                      if sorting by any factor besides sortvar (e.g.,
0119 %                      phase).
0120 %   'signif'         = [lo_dB, hi_dB, coher_signif_level] Use precomputed significance
0121 %                      thresholds (as from outputs ampsig, cohsig) to save time. {default: none}
0122 %   'caxis'          = [lo hi] Set color axis limits. Else [fraction] Set caxis limits at
0123 %                      (+/-)fraction*max(abs(data)) {default: symmetrical in dB, based on data limits}
0124 %
0125 % Add epoch-mean ERP to plot:
0126 %   'erp'      = ['on'|'off'|1|2|3|4] Plot ERP time average of the trials below the
0127 %                image.  If 'on' or 1, a single ERP (the mean of all trials) is shown.  If 2,
0128 %                two ERPs (super and sub median trials) are shown.  If 3, the trials are split into
0129 %                tertiles and their three ERPs are shown.  If 4, the trials are split into quartiles
0130 %                and four ERPs are shown. Note, if you want negative voltage plotted up, change YDIR
0131 %                to -1 in icadefs.m.  If 'erpalpha' option is used, any values of 'erp' greater than
0132 %                1 will be reset to 1. {default no ERP plotted}
0133 %   'erpalpha' = [alpha] Visualizes two-sided significance threshold (i.e., a two-tailed test) for the
0134 %                null hypothesis of a zero mean, symmetric distribution (range: [.001 0.1]). Thresholds
0135 %                are determined via a permutation test. Requires 'erp' to be a value other than 'off'.
0136 %                If 'erp' is set  to a value greater than 1, it is reset to 1 to increase plot readability.
0137 %                {default: no alpha significance thresholds plotted}
0138 %   'erpstd'   = ['on'|'off'] Plot ERP +/- stdev. Requires 'erp' {default: no std. dev. plotted}
0139 %   'erp_grid' = If 'erp_grid' is added as an option voltage axis dashed grid lines will be
0140 %                 added to the ERP plot to facilitate judging ERP amplitude
0141 %   'rmerp'    = ['on'|'off'] Subtract the average ERP from each trial before processing {default: no}
0142 %
0143 % Add time/frequency information:
0144 %   'coher'  = [freq] Plot ERP average plus mean amplitude & coherence at freq (Hz)
0145 %               Else [minfrq maxfrq] = same, but select frequency with max power in
0146 %               given range. (Note: the 'phasesort' freq (above) overwrites these
0147 %               parameters). Else [minfrq maxfrq alpha] = plot coher. signif. level line
0148 %               at probability alpha (range: [0,0.1]) {default: no coher, no alpha level}
0149 %   'srate'  = [freq] Specify the data sampling rate in Hz for amp/coher (if not
0150 %               implicit in third arg, times) {default: as defined in icadefs.m}
0151 %   'cycles' = [float] Number of cycles in the wavelet time/frequency decomposition {default: 3}
0152 %
0153 % Add plot features:
0154 %   'cbar'           = ['on'|'off'] Plot color bar to right of ERP-image {default no}
0155 %   'cbar_title'     = [string] The title for the color bar (e.g., '\muV' for
0156 %                       microvolts).
0157 %   'topo'           = {map,chan_locs,eloc_info} Plot a 2-D scalp map at upper left of image.
0158 %                       map may be a single integer, representing the plotted data channel,
0159 %                       or a vector of scalp map channel values. chan_locs may be a channel locations
0160 %                       file or a chanlocs structure (EEG.chanlocs). See '>> topoplot example'
0161 %                       eloc_info (EEG.chaninfo), if empty ([]) or absent, implies the 'X' direction
0162 %                       points towards the nose and all channels are plotted {default: no scalp map}
0163 %   'spec'           = [loHz,hiHz] Plot the mean data spectrum at upper right of image.
0164 %   'horz'           = [epochs_vector] Plot horizontal lines at specified epoch numbers.
0165 %   'vert'           = [times_vector] Plot vertical dashed lines at specified latencies
0166 %   'auxvar'         = [size(nvars,ntrials) matrix] Plot auxiliary variable(s) for each trial
0167 %                       as separate traces. Else, 'auxvar',{[matrix],{colorstrings}}
0168 %                       to specify N trace colors.  Ex: colorstrings = {'r','bo-','','k:'}
0169 %                       (See also: 'vert' and 'timewarp' above). {default: none}
0170 %   'sortvarpercent' = [float vector] Plot percentiles for the sorting variable
0171 %                       for instance, [0.1 0.5 0.9] plots the 10th percentile, the median
0172 %                       and the 90th percentile.
0173 % Plot options:
0174 % 'noxlabel'          = ['on'|'off'] Do not plot "Time (ms)" on the bottom x-axis
0175 % 'yerplabel'         = ['string'] ERP ordinate axis label (default is ERP). Print uV with '\muV'
0176 % 'avg_type'          = ['boxcar'|'Gaussian'] The type of moving average used to smooth
0177 %                        the data. 'Boxcar' smoothes the data by simply taking the mean of
0178 %                        a certain number of trials above and below each trial.
0179 %                        'Gaussian' does the same but first weights the trials
0180 %                        according to a Gaussian distribution (e.g., nearby trials
0181 %                        receive greater weight).  The Gaussian is better than the
0182 %                        boxcar in that it rather evenly filters out high frequency
0183 %                        vertical components in the ERPimage. See 'avewidth' argument
0184 %                        description for more information. {default: boxcar}
0185 % 'img_trialax_label' = ['string'] The label of the axis corresponding to trials in the ERPimage
0186 %                        (e.g., 'Reaction Time').  Note, if img_trialax_label is set to something
0187 %                        besides 'Trials' or [], the tick marks on this axis will be set in units
0188 %                        of the sorting variable.  This is a useful alternative to plotting the
0189 %                        sorting variable when the sorting variable is not in milliseconds. This
0190 %                        option is not effective if sorting by amplitude, phase, or EEG value. {default: 'Trials'}
0191 % 'img_trialax_ticks' = Vector of sorting variable values at which tick marks (e.g., [300 350 400 450]
0192 %                        for reaction time in msec) will appear on the trial axis of the erpimage. Tick mark
0193 %                        values should be given in units img_trialax_label (e.g., 'Trials' or msec).
0194 %                        This option is not effective if sorting by amplitude, phase, or EEG value.
0195 %                        {default: automatic}
0196 % 'baseline'          = [low_boundary high_boundary] a time window (in msec) whose mean amplitude in
0197 %                        each trial will be removed from each trial (e.g., [-100 0]) after filtering.
0198 %                        Useful in conjunction with 'filt' option to re-basline trials after they have been
0199 %                        filtered. Not necessary if data have already been baselined and erpimage
0200 %                        processing does not affect baseline amplitude {default: no further baselining
0201 %                        of data}
0202 % 'filt'              = [low_boundary high_boundary] a two element vector indicating the frequency
0203 %                        cut-offs for a 3rd order Butterworth filter that will be applied to each
0204 %                        trial of data.  If low_boundary=0, the filter is a low pass filter.  If
0205 %                        high_boundary=srate/2, then the filter is a high pass filter.  If both
0206 %                        boundaries are between 0 and srate/2, then the filter is a bandpass filter.
0207 %                        If both boundaries are between 0 and -srate/2, then the filter is a bandstop
0208 %                        filter (with boundaries equal to the absolute values of low_boundary and
0209 %                        high_boundary).  Note, using this option requires the 'srate' option to be
0210 %                        specified and the signal processing toolbox function butter.m.  You should
0211 %                        probably use the 'baseline' option as well since the mean prestimulus baseline
0212 %                        may no longer be 0 after the filter is applied {default: no filtering}
0213 %
0214 % Optional outputs:
0215 %    outdata  = (times,epochsout) data matrix (after smoothing)
0216 %     outvar  = (1,epochsout) actual values trials are sorted on (after smoothing).
0217 %               if 'sortvarpercent' is used, this variable contains a cell array with
0218 %               { sorted_values { sorted_percent1 ... sorted_percentN } }
0219 %   outtrials = (1,epochsout)  smoothed trial numbers
0220 %     limits  = (1,10) array, 1-9 as in 'limits' above, then analysis frequency (Hz)
0221 %    axhndls  = vector of 1-7 plot axes handles (img,cbar,erp,amp,coh,topo,spec)
0222 %        erp  = plotted ERP average
0223 %       amps  = mean amplitude time course
0224 %      coher  = mean inter-trial phase coherence time course
0225 %     cohsig  = coherence significance level
0226 %     ampsig  = amplitude significance levels [lo high]
0227 %    outamps  = matrix of imaged amplitudes (from option 'plotamps')
0228 %   phsangls  = vector of sorted trial phases at the phase-sorting frequency
0229 %     phsamp  = vector of sorted trial amplitudes at the phase-sorting frequency
0230 %    sortidx  = indices of input data epochs in the sorting order
0231 %     erpsig  = trial average significance levels [2,frames]
0232 %
0233 % Example:  >> figure;
0234 %              erpimage(data,RTs,[-400 256 256],'Test',1,1,...
0235 %                            'erp','cbar','vert',-350);
0236 %    Plots an ERP-image of 1-s data epochs sampled at 256 Hz, sorted by RTs, with
0237 %    title ('Test'), and sorted epochs not smoothed or decimated (1,1). Overplots
0238 %    the (unsmoothed) RT latencies on the colored ERP-image. Also plots the
0239 %    epoch-mean (ERP), a color bar, and a dashed vertical line at -350 ms.
0240 %
0241 % Authors: Scott Makeig, Tzyy-Ping Jung & Arnaud Delorme,
0242 %          CNL/Salk Institute, La Jolla, 3-2-1998 -
0243 %
0244 % See also: erpimages(), phasecoher(), rmbase(), cbar(), movav()
0245 
0246 %123456789012345678901234567890123456789012345678901234567890123456789012
0247 
0248 % Copyright (C) Scott Makeig & Tzyy-Ping Jung, CNL / Salk Institute, La Jolla 3-2-98
0249 %
0250 % This program is free software; you can redistribute it and/or modify
0251 % it under the terms of the GNU General Public License as published by
0252 % the Free Software Foundation; either version 2 of the License, or
0253 % (at your option) any later version.
0254 %
0255 % This program is distributed in the hope that it will be useful,
0256 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0257 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0258 % GNU General Public License for more details.
0259 %
0260 % You should have received a copy of the GNU General Public License
0261 % along with this program; if not, write to the Free Software
0262 % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
0263 
0264 % Uses external toolbox functions: phasecoher(), rmbase(), cbar(), movav()
0265 % Uses included functions:         plot1trace(), phasedet()
0266 
0267 % UNIMPLEMENTED - 'allcohers',[data2] -> image the coherences at each latency & epoch.
0268 %                 Requires arg 'coher' with alpha significance.
0269 %                 Shows projection on grand mean coherence vector at each latency
0270 %                 and trial. {default: no}
0271 
0272 %% LOG COMMENTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0273 % $Log: erpimage.m,v $
0274 %
0275 % 2010/01/01 David Groppe, Kutaslab
0276 % -Fixed error in help description of 'filt' option and orderofmag bug
0277 %
0278 % 2009/09/25 David Groppe, Kutaslab
0279 % Added:
0280 % -option avg_type (thanks to Jeffrey S. Johnson for the idea)
0281 % -option for plotting median/tertile/quartile split ERPs
0282 % -option erp_grid
0283 % -erps are plot in front of voltage=0 and time=0 lines
0284 % -Added y-axis ticks for ERP
0285 % -option cbar_title
0286 % -Removed xticks for colorbar
0287 % -option img_trialax_label
0288 % -option img_trialax_ticks
0289 % -option filt
0290 % -option baseline
0291 % -max and min ERP voltage limits can now be independently specified
0292 %  (before max=max(abs([minerp maxerp])); and min=-max)
0293 % -sortvar limits can be specified (affects visualization, not smoothing)
0294 % -erpalpha documentation is fixed (It used to say significance thresholds
0295 % were based on bootstrapping. They are actually based on a permutation
0296 % test.)
0297 % -erp_vltg_ticks
0298 %
0299 % Revision 1.283  2009/07/08 00:04:13  arno
0300 % remove warnings
0301 %
0302 % Revision 1.282  2009/07/03 18:00:58  arno
0303 % typo for octave
0304 %
0305 % Revision 1.281  2009/04/28 01:52:09  arno
0306 % header
0307 %
0308 % Revision 1.280  2008/10/16 17:06:32  julie
0309 % was missing a % in line 3387.
0310 %
0311 % Revision 1.279  2008/06/24 01:06:29  scott
0312 % in help specified 'limits'  loamp & hiamp --> lodB & hidB
0313 % changed 'bamp' to basedB for clarity in help and code
0314 %
0315 % Revision 1.278  2008/06/24 00:37:43  scott
0316 % documented the default bamp in limits - should be in dB
0317 %
0318 % Revision 1.277  2008/06/24 00:34:09  scott
0319 % corrected use of baseamp for allamps plotting -sm & rsh
0320 %
0321 % Revision 1.276  2008/06/24 00:19:27  scott
0322 % nothing
0323 %
0324 % Revision 1.275  2008/04/08 15:00:32  arno
0325 % fix phasedet baseline removal
0326 %
0327 % Revision 1.274  2008/01/11 00:28:24  arno
0328 % make all parameters 'key', 'val' while preserving backward compatibility
0329 %
0330 % Revision 1.273  2007/07/25 23:32:04  scott
0331 % restored default cycles = 3 !!?!!
0332 % the change from 3 -> 5 was never documented!?
0333 %
0334 % Revision 1.272  2007/05/10 03:44:26  toby
0335 % doc edit related to eeg_getepochevent.m changes
0336 %
0337 % Revision 1.271  2007/03/05 21:38:14  toby
0338 % Changed reference to deprecated function eventlock.m to eegalign.m
0339 %
0340 % Revision 1.270  2007/02/09 02:58:35  toby
0341 % bug fix: auto-calculation of yaxis limits occasionally fails
0342 %
0343 % Revision 1.269  2007/01/26 17:58:49  arno
0344 % Add option for compatibility with metaplottopo
0345 %
0346 % Revision 1.268  2007/01/06 13:48:54  scott
0347 % help msg corrections
0348 %
0349 % Revision 1.267  2006/09/23 23:03:15  scott
0350 % nothing
0351 %
0352 % Revision 1.266  2006/09/23 22:04:43  scott
0353 % fixed auxcolor bug, edited help message
0354 %
0355 % Revision 1.265  2006/09/22 16:25:42  scott
0356 % text edit
0357 %
0358 % Revision 1.264  2006/09/22 16:18:52  scott
0359 % added test that timewarp latencies are in ascending order in each trial
0360 %
0361 % Revision 1.263  2006/09/22 15:36:12  scott
0362 % made auxvar more flexible
0363 %
0364 % Revision 1.262  2006/09/22 02:10:19  scott
0365 % debugged 'timewarp' feature using test script /home/scott/matlab/testerpimagewarp.m
0366 %
0367 % Revision 1.261  2006/09/22 00:07:56  scott
0368 % fixed 'timewarp' call as per Diane's suggestion. -sm
0369 %
0370 % Revision 1.260  2006/07/26 03:35:13  toby
0371 % Changed so that smoothing does not bias the trial numbers
0372 %
0373 % Revision 1.259  2006/07/15 02:00:46  toby
0374 % help edit, minor code streamlining.
0375 %
0376 % Revision 1.258  2006/07/15 01:05:26  toby
0377 % Smart indent, created programming cells
0378 %
0379 % Revision 1.257  2006/05/14 20:10:29  toby
0380 % bug fix
0381 %
0382 % Revision 1.256  2006/04/12 04:12:02  scott
0383 % timeWarp -> timewarp()
0384 %
0385 % Revision 1.255  2006/04/12 03:31:32  scott
0386 % myphasecoher() -> phasecoher()
0387 %
0388 % Revision 1.254  2006/04/12 03:26:37  scott
0389 % help msg
0390 %
0391 % Revision 1.253  2006/04/12 03:05:34  scott
0392 % replaced 'timestretch' with 'timewarp', arguments in ms
0393 %
0394 % Revision 1.248  2006/01/12 01:33:04  toby
0395 % Comment edit
0396 %
0397 % Revision 1.245  2005/09/29 15:18:04  scott
0398 % documented the 'topo' flag args; set eloc_info.plotrad = [] to avoid lower-head topoplot() problem
0399 % when plotrad < radius of channel plotted -sm
0400 %
0401 % Revision 1.244  2005/03/10 16:54:32  arno
0402 % returning percentiles
0403 %
0404 % Revision 1.243  2005/03/10 16:47:52  arno
0405 % chaning variable name
0406 %
0407 % Revision 1.242  2005/03/10 02:33:38  arno
0408 % implementing percentile
0409 %
0410 % Revision 1.241  2005/03/07 22:55:21  arno
0411 % fixing error
0412 %
0413 % Revision 1.240  2005/03/07 21:19:40  arno
0414 % chaninfo
0415 %
0416 % Revision 1.239  2005/02/15 02:38:39  arno
0417 % same
0418 %
0419 % Revision 1.238  2005/02/15 02:35:39  arno
0420 % return amplitude image if 'plotamps' is used
0421 %
0422 % Revision 1.237  2005/02/14 01:35:49  arno
0423 % more debugging
0424 %
0425 % Revision 1.236  2005/02/14 01:32:47  arno
0426 % debuging ampsort
0427 %
0428 % Revision 1.235  2005/02/11 03:27:38  arno
0429 % help msg
0430 %
0431 % Revision 1.234  2005/02/10 22:19:01  arno
0432 % debug auxvar re-alignment
0433 %
0434 % Revision 1.233  2005/02/10 02:38:50  arno
0435 % auxvar re-alignement
0436 %
0437 % Revision 1.232  2005/01/21 16:45:27  scott
0438 % help msg
0439 %
0440 % Revision 1.231  2005/01/04 18:38:41  scott
0441 % function phasedet() is FAILING! just touching up phase sorting code
0442 %
0443 % Revision 1.230  2004/11/30 22:39:40  hilit
0444 % fixing 'erpstd' option
0445 %
0446 % Revision 1.229  2004/11/30 06:24:19  scott
0447 % help msg
0448 %
0449 % Revision 1.228  2004/11/22 17:57:37  scott
0450 % read YDIR from icadefs.m to set the direction of the ERP trace axes
0451 %
0452 % Revision 1.227  2004/11/15 17:22:46  scott
0453 % compute spectrum on urdata instead of data, for consistency
0454 %
0455 % Revision 1.226  2004/11/12 23:02:10  scott
0456 % making spectral peak estimates more consistent
0457 %
0458 % Revision 1.225  2004/11/05 18:36:54  arno
0459 % debug && removal
0460 %
0461 % Revision 1.224  2004/09/21 16:52:41  hilit
0462 % change && -> &
0463 %
0464 % Revision 1.223  2004/09/06 18:15:52  scott
0465 % final? time/freq window sorting debug -sm
0466 %
0467 % Revision 1.222  2004/09/03 19:02:22  scott
0468 % fixed fprintf ...s -sm
0469 %
0470 % Revision 1.221  2004/09/02 16:35:07  scott
0471 % modified 'showwin' for smoothing over amp in multiple windows
0472 % improved amp-sorting printout messages. -sm
0473 %
0474 % Revision 1.220  2004/09/02 00:04:00  scott
0475 % added 'sortwin' and ampsort by mean freq in a range. -sm
0476 %
0477 % Revision 1.219  2004/09/01 21:56:28  scott
0478 % adding 'sortwin' for 'ampsort' -sm
0479 %
0480 % Revision 1.218  2004/08/30 18:32:06  scott
0481 % added figure(curfig) before plotting functions to keep
0482 % version 7.0.0 from reverting plotting to the EEGLAB menu window
0483 % when plotting from the gui (only).  -sm
0484 %
0485 % Revision 1.217  2004/08/30 17:08:49  scott
0486 % debug last
0487 %
0488 % Revision 1.216  2004/08/30 17:05:17  scott
0489 % reshape verttimes if necessary
0490 %
0491 % Revision 1.215  2004/08/20 00:39:04  arno
0492 % debut phasearg if only using 1 frequency
0493 %
0494 % Revision 1.214  2004/08/13 20:27:12  scott
0495 % help msg
0496 %
0497 % Revision 1.213  2004/07/29 23:18:23  arno
0498 % same
0499 %
0500 % Revision 1.212  2004/07/29 23:16:50  arno
0501 % convert to double for Matlab 7
0502 %
0503 % Revision 1.211  2004/06/09 01:48:43  arno
0504 % make limits and image data consistent
0505 %
0506 % Revision 1.210  2004/06/05 01:52:45  arno
0507 % aligntime as string to avoid really realigning
0508 %
0509 % Revision 1.209  2004/06/05 01:33:07  arno
0510 % spetial option to preserve backward compatibility
0511 %
0512 % Revision 1.208  2004/05/07 04:47:17  scott
0513 % made sotvar = [] work
0514 %
0515 % Revision 1.207  2004/03/26 00:21:06  arno
0516 % same
0517 %
0518 % Revision 1.206  2004/03/26 00:19:48  arno
0519 % minimum number of trials
0520 %
0521 % Revision 1.205  2004/03/26 00:18:11  arno
0522 % plot aligntime
0523 %
0524 % Revision 1.204  2004/02/24 23:04:40  arno
0525 % fixed vertical lines in ERP when RT-aligned
0526 %
0527 % Revision 1.203  2004/01/24 22:01:23  scott
0528 % *** empty log message ***
0529 %
0530 % Revision 1.202  2004/01/24 21:58:33  scott
0531 % same
0532 %
0533 % Revision 1.201  2004/01/24 21:53:38  scott
0534 % same
0535 %
0536 % Revision 1.200  2004/01/24 21:43:59  scott
0537 % *** empty log message ***
0538 %
0539 % Revision 1.199  2004/01/24 21:36:01  scott
0540 % *** empty log message ***
0541 %
0542 % Revision 1.198  2004/01/24 21:26:50  scott
0543 % same
0544 %
0545 % Revision 1.197  2004/01/24 21:25:04  scott
0546 % same
0547 %
0548 % Revision 1.196  2004/01/24 21:22:59  scott
0549 % same
0550 %
0551 % Revision 1.195  2004/01/24 21:17:00  scott
0552 % same
0553 %
0554 % Revision 1.194  2004/01/24 21:16:32  scott
0555 % same
0556 %
0557 % Revision 1.193  2004/01/24 21:10:31  scott
0558 % same
0559 %
0560 % Revision 1.192  2004/01/24 21:08:29  scott
0561 % same
0562 %
0563 % Revision 1.191  2004/01/24 21:04:09  scott
0564 % same
0565 %
0566 % Revision 1.190  2004/01/24 21:01:18  scott
0567 % same
0568 %
0569 % Revision 1.189  2004/01/24 20:59:53  scott
0570 % same
0571 %
0572 % Revision 1.188  2004/01/24 20:51:57  scott
0573 % same
0574 %
0575 % Revision 1.187  2004/01/24 20:45:11  scott
0576 % same
0577 %
0578 % Revision 1.186  2004/01/24 20:42:05  scott
0579 % same
0580 %
0581 % Revision 1.185  2004/01/24 20:40:23  scott
0582 % plotting sorting window
0583 %
0584 % Revision 1.184  2003/12/17 21:42:13  scott
0585 % adjust same
0586 %
0587 % Revision 1.183  2003/12/17 21:40:12  scott
0588 % change y-labels on traces
0589 %
0590 % Revision 1.182  2003/12/04 17:53:52  arno
0591 % debug spec()
0592 %
0593 % Revision 1.181  2003/12/03 02:18:48  arno
0594 % use spec if psd is absent
0595 %
0596 % Revision 1.180  2003/11/26 18:16:23  scott
0597 % help msg
0598 %
0599 % Revision 1.179  2003/11/19 01:06:22  arno
0600 % including makehanning function
0601 %
0602 % Revision 1.178  2003/11/16 17:48:00  scott
0603 % plot1erp() -> plot1trace(); printf "Done."
0604 %
0605 % Revision 1.177  2003/11/14 17:01:59  scott
0606 % bootstrap msg
0607 %
0608 % Revision 1.176  2003/11/14 16:58:02  scott
0609 % debug last
0610 %
0611 % Revision 1.175  2003/11/14 16:56:16  scott
0612 % refining erpsig, sig fills
0613 %
0614 % Revision 1.174  2003/11/14 16:44:06  scott
0615 % changed signif fill color
0616 %
0617 % Revision 1.173  2003/11/14 16:34:31  scott
0618 % debug same
0619 %
0620 % Revision 1.172  2003/11/14 16:32:55  scott
0621 % debug same
0622 %
0623 % Revision 1.171  2003/11/14 16:27:27  scott
0624 % fill coher signif limits
0625 %
0626 % Revision 1.170  2003/11/13 02:32:42  scott
0627 % fill the dB signif limits
0628 %
0629 % Revision 1.169  2003/11/13 02:29:51  scott
0630 % debug
0631 %
0632 % Revision 1.168  2003/11/13 02:28:12  scott
0633 % debug
0634 %
0635 % Revision 1.167  2003/11/13 02:26:25  scott
0636 % debug
0637 %
0638 % Revision 1.166  2003/11/13 02:25:29  scott
0639 % debug
0640 %
0641 % Revision 1.165  2003/11/13 02:24:19  scott
0642 % fill amp signif
0643 %
0644 % Revision 1.164  2003/11/13 02:14:16  scott
0645 % same
0646 %
0647 % Revision 1.163  2003/11/13 01:58:10  scott
0648 % same
0649 %
0650 % Revision 1.162  2003/11/13 01:47:13  scott
0651 % make erpalpha fill less saturated
0652 %
0653 % Revision 1.161  2003/11/10 23:40:32  scott
0654 % fill the bootstrap limits behind the erp if erpalpha
0655 %
0656 % Revision 1.160  2003/11/10 16:57:46  arno
0657 % msg typo
0658 %
0659 % Revision 1.159  2003/10/29 22:15:52  scott
0660 % adjust erp alpha NBOOT to (low) alpha
0661 %
0662 % Revision 1.158  2003/10/29 22:07:10  arno
0663 % nothing
0664 %
0665 % Revision 1.157  2003/09/24 19:36:19  scott
0666 % fixed auxvar bug
0667 %
0668 % Revision 1.156  2003/09/24 18:56:22  scott
0669 % debug
0670 %
0671 % Revision 1.155  2003/09/24 00:45:30  scott
0672 % debug same
0673 %
0674 % Revision 1.154  2003/09/24 00:43:15  scott
0675 % debug same
0676 %
0677 % Revision 1.153  2003/09/24 00:42:10  scott
0678 % debug same
0679 %
0680 % Revision 1.152  2003/09/24 00:39:05  scott
0681 % adding 'horz' -> horizontal line plotting
0682 %
0683 % Revision 1.151  2003/09/21 21:12:27  scott
0684 % edited comments
0685 %
0686 % Revision 1.150  2003/09/11 22:23:25  scott
0687 % debug same
0688 %
0689 % Revision 1.149  2003/09/11 22:19:49  scott
0690 % made 'plotamps' and 'auxvar' work together
0691 %
0692 % Revision 1.148  2003/09/09 23:26:48  arno
0693 % change && to &
0694 %
0695 % Revision 1.147  2003/09/07 00:41:51  arno
0696 % fixing stdev, do not know why it crashed
0697 %
0698 % Revision 1.146  2003/09/06 22:45:04  scott
0699 % same
0700 %
0701 % Revision 1.145  2003/09/06 22:43:57  scott
0702 % add erpsig output
0703 %
0704 % Revision 1.144  2003/09/06 22:24:55  scott
0705 % debug last, add \n before first printed line "Plotting input...
0706 %
0707 % Revision 1.143  2003/09/06 22:15:23  scott
0708 % adjust auxvar if phase sort or amp sort
0709 %
0710 % Revision 1.142  2003/08/27 17:45:07  scott
0711 % header
0712 %
0713 % Revision 1.141  2003/08/25 22:38:25  scott
0714 % auxvars adjust tests -sm
0715 %
0716 % Revision 1.140  2003/08/24 04:49:01  scott
0717 % help msg
0718 %
0719 % Revision 1.139  2003/08/24 04:38:18  scott
0720 % same
0721 %
0722 % Revision 1.138  2003/08/24 04:37:47  scott
0723 % debug last
0724 %
0725 % Revision 1.137  2003/08/24 04:36:29  scott
0726 % fprintf
0727 %
0728 % Revision 1.136  2003/08/24 04:35:19  scott
0729 % added help for 'erpalpha' -sm
0730 %
0731 % Revision 1.135  2003/08/24 04:27:41  scott
0732 % fprintf adjust
0733 %
0734 % Revision 1.134  2003/08/24 04:25:05  scott
0735 % adjust erpalpha line type
0736 %
0737 % Revision 1.133  2003/08/24 04:23:40  scott
0738 % debug
0739 %
0740 % Revision 1.132  2003/08/24 04:20:26  scott
0741 % same
0742 %
0743 % Revision 1.131  2003/08/24 04:19:58  scott
0744 % added fprintf
0745 %
0746 % Revision 1.130  2003/08/24 04:17:42  scott
0747 % same
0748 %
0749 % Revision 1.129  2003/08/24 04:17:17  scott
0750 % same
0751 %
0752 % Revision 1.128  2003/08/24 04:16:39  scott
0753 % debug same
0754 %
0755 % Revision 1.127  2003/08/24 04:12:33  scott
0756 % added (erpalpha) bootstrap ERP estimation -sm
0757 %
0758 % Revision 1.126  2003/07/26 17:18:51  scott
0759 % debug same
0760 %
0761 % Revision 1.125  2003/07/26 17:16:50  scott
0762 % debug same
0763 %
0764 % Revision 1.124  2003/07/26 17:14:12  scott
0765 % debug same
0766 %
0767 % Revision 1.123  2003/07/26 17:05:26  scott
0768 % debug auxvars warning, hide vert matrix option
0769 %
0770 % Revision 1.122  2003/07/24 23:41:05  arno
0771 % correcting output
0772 %
0773 % Revision 1.121  2003/07/24 18:17:33  scott
0774 % *** empty log message ***
0775 %
0776 % Revision 1.120  2003/07/22 16:03:06  scott
0777 % phasesort help message
0778 %
0779 % Revision 1.119  2003/07/22 15:40:26  scott
0780 % adding circular smoothing to phase-sorted, allamps plots -sm
0781 %
0782 % Revision 1.118  2003/07/22 00:52:02  scott
0783 % debug
0784 %
0785 % Revision 1.117  2003/07/22 00:48:36  scott
0786 % debug
0787 %
0788 % Revision 1.116  2003/07/22 00:44:21  scott
0789 % debug
0790 %
0791 % Revision 1.115  2003/07/21 21:38:32  scott
0792 % debug
0793 %
0794 % Revision 1.114  2003/07/21 21:36:36  scott
0795 % debug
0796 %
0797 % Revision 1.113  2003/07/21 21:33:48  scott
0798 % debug
0799 %
0800 % Revision 1.112  2003/07/21 21:25:44  scott
0801 % debug
0802 %
0803 % Revision 1.111  2003/07/21 21:19:15  scott
0804 % debug
0805 %
0806 % Revision 1.110  2003/07/21 21:17:52  scott
0807 % debug
0808 %
0809 % Revision 1.109  2003/07/21 21:13:09  scott
0810 % debug
0811 %
0812 % Revision 1.108  2003/07/21 21:09:59  scott
0813 % debug
0814 %
0815 % Revision 1.107  2003/07/21 20:46:19  scott
0816 % debug
0817 %
0818 % Revision 1.106  2003/07/21 20:45:01  scott
0819 % fg
0820 %
0821 % debug last
0822 %
0823 % Revision 1.105  2003/07/21 20:43:15  scott
0824 % using wraparound smoothing for phase-sorted epochs
0825 %
0826 % Revision 1.104  2003/07/15 18:01:05  scott
0827 % trial number -> trials throughout
0828 %
0829 % Revision 1.103  2003/05/06 17:31:29  arno
0830 % debug Rmerp
0831 %
0832 % Revision 1.102  2003/05/06 15:58:51  scott
0833 % header edits
0834 %
0835 % Revision 1.101  2003/05/06 15:54:15  scott
0836 % edit header
0837 %
0838 % Revision 1.100  2003/05/06 00:49:52  arno
0839 % debug last
0840 %
0841 % Revision 1.99  2003/05/06 00:45:43  arno
0842 % implementing new option rmerp
0843 %
0844 % Revision 1.98  2003/04/26 01:06:38  arno
0845 % debuging ampsort for phase sorting
0846 %
0847 % Revision 1.97  2003/04/25 22:32:26  arno
0848 % doing the same for ampsort
0849 %
0850 % Revision 1.96  2003/04/25 22:30:40  arno
0851 % interpolating phase value
0852 %
0853 % Revision 1.95  2003/04/25 18:04:31  arno
0854 % nothing
0855 %
0856 % Revision 1.94  2003/04/24 21:54:24  arno
0857 % removing debug message
0858 %
0859 % Revision 1.93  2003/04/23 23:59:48  arno
0860 % remove debug msg
0861 %
0862 % Revision 1.92  2003/04/23 23:58:02  arno
0863 % debuging phasedet in erpimage
0864 %
0865 % Revision 1.91  2003/04/23 23:41:19  arno
0866 % restoring cycles to a default of 3, adding cycle paramete
0867 % r
0868 %
0869 % Revision 1.90  2003/04/23 22:09:21  arno
0870 % adding cycles to the phasedet script
0871 %
0872 % Revision 1.89  2003/04/23 01:24:15  arno
0873 % chaning default to 3 cycles at 5 Hz
0874 %
0875 % Revision 1.88  2003/04/23 01:18:08  arno
0876 % typo
0877 %
0878 % Revision 1.87  2003/04/23 01:15:51  arno
0879 % debuging cycles
0880 %
0881 % Revision 1.86  2003/04/23 01:10:39  arno
0882 % remove 3-cycle from help
0883 %
0884 % Revision 1.85  2003/04/23 01:05:55  arno
0885 % changing default cycle for coherence
0886 %
0887 % Revision 1.84  2003/04/18 17:20:13  arno
0888 % noshow option
0889 %
0890 % Revision 1.83  2003/03/23 22:46:39  scott
0891 % *** empty log message ***
0892 %
0893 % Revision 1.82  2003/03/15 18:01:34  scott
0894 % help msg
0895 %
0896 % Revision 1.81  2003/03/13 03:20:33  scott
0897 % restoring
0898 %
0899 % Revision 1.80  2003/03/07 22:21:46  scott
0900 % same
0901 %
0902 % Revision 1.79  2003/03/07 22:15:19  scott
0903 % same
0904 %
0905 % Revision 1.78  2003/03/07 22:13:51  scott
0906 % same
0907 %
0908 % Revision 1.77  2003/03/07 22:11:58  scott
0909 % debugging axcopy call
0910 %
0911 % Revision 1.76  2003/03/07 22:09:47  scott
0912 % typo
0913 %
0914 % Revision 1.75  2003/03/07 22:08:30  scott
0915 % axcopy - adding ''''
0916 %
0917 % Revision 1.74  2003/03/07 21:18:13  scott
0918 % testing axcopy
0919 %
0920 % Revision 1.73  2003/03/07 21:15:25  scott
0921 % adding eegplot option if click on erpimage -sm
0922 %
0923 % Revision 1.72  2003/03/04 20:29:18  arno
0924 % header typo
0925 %
0926 % Revision 1.71  2003/01/31 23:23:34  arno
0927 % adding erpstd for ploting standard deviation of ERP
0928 %
0929 % Revision 1.70  2003/01/02 16:55:02  scott
0930 % allowed sortvars to be negative -sm
0931 %
0932 % Revision 1.69  2002/11/19 23:27:55  arno
0933 % debugging spectrum for very long epochs
0934 %
0935 % Revision 1.68  2002/11/09 21:19:57  scott
0936 % added example
0937 %
0938 % Revision 1.67  2002/11/09 21:07:29  scott
0939 % reorganized help message optional argument list
0940 %
0941 % Revision 1.66  2002/10/15 17:52:00  scott
0942 % help msg - only two args required
0943 %
0944 % Revision 1.65  2002/10/14 16:12:38  scott
0945 % fprintf ms
0946 %
0947 % Revision 1.64  2002/10/14 16:09:53  scott
0948 % valsort fprintf
0949 %
0950 % Revision 1.63  2002/10/14 16:07:29  scott
0951 % removed moreargs - consolidated help message -sm
0952 %
0953 % Revision 1.62  2002/10/14 14:56:45  scott
0954 % working on ampsort
0955 %
0956 % Revision 1.61  2002/10/14 00:42:58  scott
0957 % added valsort direction and help msg -sm
0958 %
0959 % Revision 1.60  2002/10/13 23:56:54  scott
0960 % *** empty log message ***
0961 %
0962 % Revision 1.59  2002/10/13 23:55:47  scott
0963 % debugging valsort
0964 %
0965 % Revision 1.58  2002/10/13 23:51:22  scott
0966 % edit valsort fprint
0967 %
0968 % Revision 1.57  2002/10/13 23:49:43  scott
0969 % *** empty log message ***
0970 %
0971 % Revision 1.56  2002/10/13 23:48:30  scott
0972 % *** empty log message ***
0973 %
0974 % Revision 1.55  2002/10/13 23:48:04  scott
0975 % *** empty log message ***
0976 %
0977 % Revision 1.54  2002/10/13 23:46:37  scott
0978 % *** empty log message ***
0979 %
0980 % Revision 1.53  2002/10/13 23:43:47  scott
0981 % debugging
0982 %
0983 % Revision 1.52  2002/10/13 23:41:57  scott
0984 % debug ampargs
0985 %
0986 % Revision 1.51  2002/10/13 23:41:11  scott
0987 % debug ampargs
0988 %
0989 % Revision 1.50  2002/10/13 23:37:01  scott
0990 % debug valsort
0991 %
0992 % Revision 1.49  2002/10/13 23:35:51  scott
0993 % debug
0994 %
0995 % Revision 1.48  2002/10/13 23:33:10  scott
0996 % valsort debug
0997 %
0998 % Revision 1.47  2002/10/13 23:28:18  scott
0999 % added Ampflag, Valflag defaults -sm
1000 %
1001 % Revision 1.46  2002/10/13 23:24:48  scott
1002 % typo
1003 %
1004 % Revision 1.45  2002/10/13 23:23:32  scott
1005 % added ampsort, valsort args -sm
1006 %
1007 % Revision 1.44  2002/09/03 21:35:27  arno
1008 % removing oridata completelly
1009 %
1010 % Revision 1.43  2002/09/03 21:15:23  scott
1011 % go back to averaging urdata instead of oridata -sm
1012 %
1013 % Revision 1.42  2002/08/31 17:00:50  arno
1014 % add yerplabel option
1015 %
1016 % Revision 1.41  2002/08/30 18:18:09  arno
1017 % same
1018 %
1019 % Revision 1.40  2002/08/30 18:14:08  arno
1020 % same
1021 %
1022 % Revision 1.39  2002/08/30 18:12:12  arno
1023 % same
1024 %
1025 % Revision 1.38  2002/08/30 18:09:15  arno
1026 % same
1027 %
1028 % Revision 1.37  2002/08/30 18:07:20  arno
1029 % same.
1030 %
1031 % Revision 1.36  2002/08/30 18:01:11  arno
1032 % same
1033 %
1034 % Revision 1.35  2002/08/30 18:00:22  arno
1035 % debug erp axis and average
1036 %
1037 % Revision 1.34  2002/08/21 18:36:13  arno
1038 % adding error message
1039 %
1040 % Revision 1.33  2002/08/19 19:48:25  arno
1041 % commenting crosscoher for Mac compatibility
1042 %
1043 % Revision 1.32  2002/08/12 01:37:00  arno
1044 % color
1045 %
1046 % Revision 1.31  2002/08/11 22:34:04  arno
1047 % color
1048 %
1049 % Revision 1.30  2002/08/09 16:28:07  arno
1050 % debugging allamps
1051 %
1052 % Revision 1.29  2002/08/05 18:04:55  arno
1053 % performing moving average on allamps amplitude (and not log)
1054 %
1055 % Revision 1.28  2002/07/27 01:25:45  arno
1056 % debugging vert
1057 %
1058 % Revision 1.27  2002/07/26 16:18:32  arno
1059 % removing debugging messages
1060 %
1061 % Revision 1.26  2002/07/26 16:14:03  arno
1062 % removing trials with Nan values for sortvar, debugging 'vert'
1063 %
1064 % Revision 1.25  2002/07/15 02:00:48  arno
1065 % same
1066 %
1067 % Revision 1.24  2002/07/15 01:55:45  arno
1068 % debugging minamp maxamp
1069 %
1070 % Revision 1.23  2002/07/15 01:48:20  arno
1071 % force yscale on amp ploterp
1072 %
1073 % Revision 1.22  2002/07/14 03:17:57  arno
1074 % making allamps moving average of log power
1075 %
1076 % Revision 1.21  2002/07/14 02:39:16  arno
1077 % same
1078 %
1079 % Revision 1.20  2002/07/14 02:23:16  arno
1080 % same
1081 %
1082 % Revision 1.19  2002/07/14 01:58:00  arno
1083 % same
1084 %
1085 % Revision 1.18  2002/07/14 01:47:23  arno
1086 % testing amps limits
1087 %
1088 % Revision 1.17  2002/05/23 16:59:09  scott
1089 % replaced nanmean with nan_mean() -sm
1090 %
1091 % Revision 1.16  2002/05/22 06:06:28  marissa
1092 % changed line 1047 to remove nonexistent variable 'baseall'
1093 %
1094 % Revision 1.15  2002/05/20 17:58:53  scott
1095 % adding fprintf info about ampsig plotting -sm
1096 %
1097 % Revision 1.14  2002/04/25 17:52:51  arno
1098 % removing debug message
1099 %
1100 % Revision 1.13  2002/04/25 17:51:50  arno
1101 % including renorm inside function
1102 %
1103 % Revision 1.12  2002/04/25 17:08:03  arno
1104 % correcting error message
1105 %
1106 % Revision 1.11  2002/04/24 19:04:45  arno
1107 % further check for coherfreq
1108 %
1109 % Revision 1.10  2002/04/24 18:29:21  arno
1110 % two inputs for coher and time centering for phase
1111 %
1112 % Revision 1.9  2002/04/24 17:35:50  arno
1113 % rechanged log
1114 %
1115 % 3/5/98 added nosort option -sm
1116 % 3/22/98 added colorbar ylabel, sym. range finding -sm
1117 % 5/08/98 added noplot option -sm
1118 % 6/09/98 added align, erp, coher options -sm
1119 % 6/10/98 added limits options -sm
1120 % 6/26/98 made number of variables output 8, as above -sm
1121 % 9/16/98 plot out-of-bounds sortvars at nearest times boundary -sm
1122 % 10/27/98 added cohsig, alpha -sm
1123 % 10/28/98 adjust maxamp, maxcoh computation -sm
1124 % 05/03/99 added horizontal ticks beneath coher trace, fixed vert.
1125 %          scale printing -t-pj
1126 % 05/07/99 converted amps plot to log scaling -sm
1127 % 05/14/99 added sort by alpha phase -se
1128 % 07/23/99 made "Phase-sorted" axis label -sm
1129 % 07/24/99 added 'allamps' option -sm
1130 % 08/04/99 added new times spec., 'srate' arg, made 'phase' and 'allamps'
1131 %          work together, plot re-aligned time zeros  -sm
1132 % 06/26/99 debugged cbar; added vert lines at aligntime to plot1erp() axes -sm
1133 % 09/29/99 fixed srate computation from times -sm & se
1134 % 01/18/00 output outsort without clipping -sm
1135 % 02/29/00 added 'vert' arg, fixed xticklabels, added ampsig -sm
1136 % 03/03/00 added 'vert' arg lines to erp/amp/coher axes -sm
1137 % 03/17/00 added axcopy -sm & tpj
1138 % 03/29/00 added 'topo' option -sm
1139 % 05/05/00 fixed y-axis label bug when time limits given -sm
1140 % 06/01/00 added topphase arg to 'phase' option for phasemovie.m -sm
1141 % 07/12/00 adjusted prctle()
1142 % 07/12/00 added 'spec' option -sm
1143 % 08/22/00 added coherfreq to limits output -sm
1144 % 09/13/00 added hard limit (1) to maxcoh -sm
1145 % 09/14/00 made topoplot() and psd() plots relative to gca, not gcf -sm
1146 % 10/10/00 added NoTimeflag -sm
1147 % 11/03/00 changed method of rejecting small amplitude trials for phase sort -sm
1148 % 11/03/00 added number_of_trials_out option for decfactor -sm
1149 % 11/16/00 added ampoffset to center sig lines around baseline mean amp (0) -sm
1150 % 01/06/01 edited help message; adjusted ampsig plot limits;initialized outputs -sm
1151 %          rm'd 'allcoher' from help message - not fully implemented -sm
1152 % 01/09/01 documented limits arg 'bamp' (baseline amplitude) -sm
1153 % 02/13/01 debugged use of stored baseamp, ampsig parameters -sm
1154 % 03/28/01 made erpimage(data) possible. Debugged ampsig change limits -sm
1155 % 08/31/01 fixed allamps bug -sm
1156 % 09/03/01 added 'auxvar' plotting -sm
1157 % 01-25-02 reformated help & license, added links -ad
1158 % 02-16-02 added matrix option to arg 'vert' -sm
1159 % 04-05-02 corrected zero alignment problem (display only) -ad
1160 %
1161 % Known Bugs:
1162 % 'limits', [lotime hitime] may not work with 'erp'
1163 % 'limits', [... loerp hierp] (still??) may leave "ghost" grey numbers
1164 %       on the coher axis when printed (-djpeg or -depsc)
1165 % 'allcohers' - not fully implemented, and has been omitted from the help msg
1166 
1167 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)
1168 
1169 %
1170 %% %%%%%%%%%%%%%%%%% Define defaults %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1171 %
1172 % Initialize optional output variables:
1173 warning off;
1174 erp = []; amps = []; cohers = []; cohsig = []; ampsig = [];
1175 allamps = []; phaseangles = []; phsamp = []; sortidx = [];
1176 auxvar = []; erpsig = []; winloc = [];winlocs = [];
1177 timeStretchColors = {};
1178 curfig = gcf;   % note current figure - to avoid v7.0.0 bug that draws
1179 % some elements on the EEGLAB window -sm 8-30-04
1180 
1181 YES = 1;  % logical variables
1182 NO  = 0;
1183 
1184 DEFAULT_BASELINE_END = 0; % ms
1185 TIMEX = 1;          % 1 -> plot time on x-axis;
1186 % 0 -> plot trials on x-axis
1187 
1188 BACKCOLOR = [0.8 0.8 0.8]; % grey background
1189 try, icadefs; catch, end;
1190 % read BACKCOLOR for plot from defs file (edit this)
1191 % read DEFAULT_SRATE for coher,phase,allamps, etc.
1192 % read YDIR for plotting ERP
1193 % Fix plotting text and line style parameters
1194 SORTWIDTH = 2.5;    % Linewidth of plotted sortvar
1195 ZEROWIDTH = 3.0;    % Linewidth of vertical 0 line
1196 VERTWIDTH = 2.5;    % Linewidth of optional vertical lines
1197 HORZWIDTH = 2.1;    % Linewidth of optional vertical lines
1198 SIGNIFWIDTH = 1.9;  % Linewidth of red significance lines for amp, coher
1199 DOTSTYLE   = 'k--'; % line style to use for vetical dotted/dashed lines
1200 LINESTYLE = '-';    % solid line
1201 LABELFONT = 10;     % font sizes for axis labels, tick labels
1202 TICKFONT  = 10;
1203 
1204 PLOT_HEIGHT = 0.2;  % fraction of y dim taken up by each time series axes
1205 YGAP = 0.03;        % fraction gap between time axes
1206 YEXPAND = 1.3;      % expansion factor for y-axis about erp, amp data limits
1207 
1208 DEFAULT_SDEV  = 1/7; % smooth trials with this window size by default if Gaussian window
1209 DEFAULT_AVEWIDTH  = 1; % smooth trials with this window size by default
1210 DEFAULT_DECFACTOR = 1; % decimate by this factor by default
1211 DEFAULT_CYCLES    = 3; % use this many cycles in amp,coher computation window
1212 cycles = DEFAULT_CYCLES;
1213 DEFAULT_CBAR      = NO;% do not plot color bar by default
1214 DEFAULT_PHARGS = [0 25 8 13]; % Default arguments for phase sorting
1215 DEFAULT_ALPHA     = 0.01;
1216 alpha     = 0;      % default alpha level for coherence significance
1217 
1218 MIN_ERPALPHA = 0.001; % significance bounds for ERP
1219 MAX_ERPALPHA = 0.1;
1220 
1221 Noshow    = NO;     % show sortvar by default
1222 Nosort    = NO;     % sort on sortvar by default
1223 Caxflag   = NO;     % use default caxis by default
1224 
1225 timestretchflag = NO; % Added -JH
1226 
1227 mvavg_type='boxcar'; % use the original rectangular moving average -DG
1228 erp_grid = NO; % add y-tick grids to ERP plot -DG
1229 cbar_title = []; % title to add above ERPimage color bar (e.g., '\muV') -DG
1230 img_ylab = 'Trials'; % make the ERPimage y-axis in units of the sorting variable -DG
1231 img_ytick_lab = []; % the values at which tick marks will appear on the trial axis of the ERPimage (the y-axis by default).
1232 %Note, this is in units of the sorting variable if img_ylab~='Trials', otherwise it is in units of trials -DG
1233 baseline = []; %time window of each trial whose mean amp will be used to baseline the trial -DG
1234 flt=[]; %frequency domain filter parameters -DG
1235 sortvar_limits=[]; %plotting limits for sorting variable/trials; limits only affect visualization, not smoothing -DG
1236 replace_ties = NO; %if YES, trials with the exact same value of a sorting variable will be replaced by their average -DG
1237 erp_vltg_ticks=[]; %if non-empty, these are the voltage axis ticks for the ERP plot
1238 
1239 Caxis     = [];
1240 caxfraction = [];
1241 Coherflag = NO;     % don't compute or show amp,coher by default
1242 Cohsigflag= NO;     % default: do not compute coherence significance
1243 Allampsflag=NO;     % don't image the amplitudes by default
1244 Allcohersflag=NO;   % don't image the coherence amplitudes by default
1245 Topoflag  = NO;     % don't plot a topoplot in upper left
1246 Specflag  = NO;     % don't plot a spectrum in upper right
1247 Erpflag   = NO;     % don't show erp average by default
1248 Erpstdflag= NO;
1249 Erpalphaflag= NO;
1250 Alignflag = NO;     % don't align data to sortvar by default
1251 Colorbar  = NO;     % if YES, plot a colorbar to right of erp image
1252 Limitflag = NO;     % plot whole times range by default
1253 Phaseflag = NO;     % don't sort by phase
1254 Ampflag   = NO;     % don't sort by amplitude
1255 Sortwinflag = NO;   % sort by amplitude over a window
1256 Valflag   = NO;     % don't sort by value
1257 Srateflag = NO;     % srate not given
1258 Vertflag  = NO;
1259 Horzflag  = NO;
1260 titleflag = NO;
1261 Noshowflag  = NO;
1262 Renormflag = NO;
1263 Showwin = NO;
1264 yerplabel = 'ERP';
1265 yerplabelflag = NO;
1266 verttimes = [];
1267 horzepochs = [];
1268 NoTimeflag= NO;     % by default DO print "Time (ms)" below bottom axis
1269 Signifflag= NO;     % compute significance instead of receiving it
1270 Auxvarflag= NO;
1271 plotmodeflag= NO;
1272 plotmode = 'normal';
1273 Cycleflag = NO;
1274 signifs   = NaN;
1275 coherfreq = nan;    % amp/coher-calculating frequency
1276 freq = 0;           % phase-sorting frequency
1277 srate = DEFAULT_SRATE; % from icadefs.m
1278 aligntime = nan;
1279 timelimits= nan;
1280 topomap   = [];     % topo map vector
1281 lospecHz  = [];     % spec lo frequency
1282 topphase = 180;     % default top phase for 'phase' option
1283 renorm    = 'no';
1284 noshow    = 'no';
1285 Rmerp     = 'no';
1286 percentiles = [];
1287 percentileflag = NO;
1288 
1289 minerp = NaN; % default limits
1290 maxerp = NaN;
1291 minamp = NaN;
1292 maxamp = NaN;
1293 mincoh = NaN;
1294 maxcoh = NaN;
1295 baseamp =NaN;
1296 allamps = []; % default return
1297 
1298 ax1    = NaN; % default axes handles
1299 axcb   = NaN;
1300 ax2    = NaN;
1301 ax3    = NaN;
1302 ax4    = NaN;
1303 
1304 
1305 timeStretchRef = [];
1306 timeStretchMarks = [];
1307 tsurdata = []; % time-stretched data before smoothing (for timestretched-erp computation)
1308 % If time-stretching is off, this variable remains empty
1309 %
1310 %% %%%%%%%%%%%%%%%%% Test, fill in commandline args %%%%%%%%%%%%%%%%%%%%%%%%%%%%
1311 %
1312 if nargin < 1
1313     help erpimage
1314     return
1315 end
1316 
1317 if nargin < 3 | isempty(times)
1318     if size(data,1)==1 | size(data,2)==1
1319         fprintf('erpimage(): either input a times vector or make data size = (frames,trials).\n')
1320         return
1321     end
1322     times = 1:size(data,1);
1323     NoTimesPassed= 1;
1324 end
1325 
1326 if nargin < 2 | isempty(sortvar)
1327     sortvar = 1:size(data,2);
1328     Noshow = 1; % don't plot the dummy sortvar
1329 end
1330 
1331 framestot = size(data,1)*size(data,2);
1332 ntrials = length(sortvar);
1333 if ntrials < 2
1334     help erpimage
1335     fprintf('\nerpimage(): too few trials.\n');
1336     return
1337 end
1338 
1339 frames = floor(framestot/ntrials);
1340 if frames*ntrials ~= framestot
1341     help erpimage
1342     fprintf(...
1343         '\nerpimage(); length of sortvar doesn''t divide number of data elements??\n')
1344     return
1345 end
1346 
1347 if nargin < 6
1348     decfactor = 0;
1349 end
1350 if nargin < 5
1351     avewidth = DEFAULT_AVEWIDTH;
1352 end
1353 if nargin<4
1354     titl = ''; % default no title
1355 end
1356 if nargin<3
1357     times = NO;
1358 end
1359 if (length(times) == 1) | (times == NO),  % make default times
1360     times = 0:frames-1;
1361     srate = 1000*(length(times)-1)/(times(length(times))-times(1));
1362     fprintf('Using sampling rate %g Hz.\n',srate);
1363 elseif length(times) == 3
1364     mintime = times(1);
1365     frames = times(2);
1366     srate = times(3);
1367     times = mintime:1000/srate:mintime+(frames-1)*1000/srate;
1368     fprintf('Using sampling rate %g Hz.\n',srate);
1369 else
1370     % Note: might use default srate read from icadefs here...
1371     srate = 1000*(length(times)-1)/(times(end)-times(1));
1372 end
1373 if length(times) ~= frames
1374     fprintf(...
1375         'erpimage(): length(data)(%d) ~= length(sortvar)(%d) * length(times)(%d).\n\n',...
1376         framestot,              length(sortvar),   length(times));
1377     return
1378 end
1379 
1380 if decfactor == 0,
1381     decfactor = DEFAULT_DECFACTOR;
1382 elseif decfactor > ntrials
1383     fprintf('Setting variable decfactor to max %d.\n',ntrials)
1384     decfactor = ntrials;
1385 end
1386 %
1387 %% %%%%%%%%%%%%%%% Collect optional args %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1388 %
1389 if nargin > 6
1390     flagargs = [];
1391 
1392     a = 6;
1393     while a < nargin % for each remaining Arg
1394         a = a + 1;
1395         
1396         Arg = eval(['arg' int2str(a-6)]);
1397         if Caxflag == YES
1398             if size(Arg,1) ~= 1 || size(Arg,2) > 2
1399                 help erpimage
1400                 fprintf('\nerpimage(): caxis arg must be a scalar or (1,2) vector.\n');
1401                 return
1402             end
1403             if size(Arg,2) == 2
1404                 Caxis = Arg;
1405             else
1406                 caxfraction = Arg;
1407             end
1408             Caxflag = NO;
1409 
1410         elseif timestretchflag == YES % Added -JH
1411             timeStretchMarks = Arg{1};
1412             timeStretchMarks = round(1+(timeStretchMarks-times(1))*srate/1000); % convert from ms to frames -sm
1413             [smc smr] = find(diff(timeStretchMarks') < 0);
1414             if ~isempty(smr)
1415                  fprintf('erpimage(): Timewarp event latencies not in ascending order in trial %d.\n',smr)
1416                  return
1417             end
1418 
1419             timeStretchMarks = [ ...
1420                 repmat(1, [size(timeStretchMarks,1), 1]), ...% Epoch begins
1421                 timeStretchMarks, ...
1422                 repmat(length(times), [size(timeStretchMarks,1), 1])]; % Epoch ends
1423             if length(Arg) < 2 || isempty(Arg{2})
1424                 timeStretchRef = median(timeStretchMarks);
1425             else
1426                 timeStretchRef = Arg{2};
1427                 timeStretchRef = round(1+(timeStretchRef-times(1))*srate/1000); % convert from ms to frames -sm
1428                 timeStretchRef = [1 timeStretchRef length(times)];  % add epoch beginning, end
1429             end
1430             if length(Arg) < 3 || isempty(Arg{3})
1431                 timeStretchColors = {};
1432             else
1433                 timeStretchColors = Arg{3};
1434             end
1435             fprintf('The %d events specified in each trial will be time warped to latencies:',length(timeStretchRef)-2);
1436             fprintf(' %.0f', times(1)+1000*(timeStretchRef(2:end-1)-1)/srate); % converted from frames to ms -sm
1437             fprintf(' ms\n');
1438             timestretchflag = NO;
1439         elseif Coherflag == YES
1440             if length(Arg) > 3 || length(Arg) < 1
1441                 help erpimage
1442                 fprintf('\nerpimage(): coher arg must be size <= 3.\n');
1443                 return
1444             end
1445             coherfreq = Arg(1);
1446             if size(Arg,2) == 1
1447                 coherfreq = Arg(1);
1448             else
1449                 coherfreq = Arg(1:2);
1450             end;
1451             if size(Arg,2) == 3
1452                 Cohsigflag = YES;
1453                 alpha  = Arg(3);
1454                 if alpha < 0 || alpha > 0.1
1455                     fprintf('erpimage(): alpha value %g out of bounds.\n',alpha);
1456                     return
1457                 end
1458             end
1459             Coherflag = NO;
1460             Erpflag = YES;  % plot amp, coher below erp time series
1461         elseif Topoflag == YES;
1462             if length(Arg) < 2
1463                 help erpimage
1464                 fprintf('\nerpimage(): topo arg must be a list of length 2 or 3.\n');
1465                 return
1466             end
1467             topomap = Arg{1};
1468             eloc_file = Arg{2};
1469             if length(Arg) > 2, eloc_info = Arg{3};
1470             else                eloc_info = [];
1471             end;
1472             Topoflag = NO;
1473         elseif Specflag == YES;
1474             if length(Arg) ~= 2
1475                 help erpimage
1476                 fprintf('\nerpimage(): spec arg must be a list of length 2.\n');
1477                 return
1478             end
1479             lospecHz = Arg(1);
1480             hispecHz = Arg(2);
1481             Specflag = NO;
1482         elseif Renormflag == YES
1483             renorm = Arg;
1484             Renormflag = NO;
1485         elseif Noshowflag == YES
1486             noshow = Arg;
1487             Noshowflag = NO;
1488         elseif Alignflag == YES
1489             aligntime = Arg;
1490             Alignflag = NO;
1491         elseif percentileflag == YES
1492             percentiles = Arg;
1493             percentileflag = NO;
1494         elseif Limitflag == YES
1495             %  [lotime hitime loerp hierp loamp hiamp locoher hicoher]
1496             if size(Arg,1) ~= 1 || size(Arg,2) < 2 || size(Arg,2) > 9
1497                     help erpimage
1498                 fprintf('\nerpimage(): limits arg must be a vector sized (1,2<->9).\n');
1499                 return
1500             end
1501             if  ~isnan(Arg(1)) & (Arg(2) <= Arg(1))
1502                 help erpimage
1503                 fprintf('\nerpimage(): time limits out of order or out of range.\n');
1504                 return
1505             end
1506             if Arg(1) < min(times)
1507                 Arg(1) = min(times); 
1508                 fprintf('Adjusting mintime limit to first data value %g\n',min(times));
1509             end
1510             if Arg(2) > max(times)
1511                 Arg(2) = max(times); 
1512                 fprintf('Adjusting maxtime limit to last data value %g\n',max(times));
1513             end
1514             timelimits = Arg(1:2);
1515             if length(Arg)> 2
1516                 minerp = Arg(3);
1517             end
1518             if length(Arg)> 3
1519                 maxerp = Arg(4);
1520             end
1521             if ~isnan(maxerp) & maxerp <= minerp
1522                 help erpimage
1523                 fprintf('\nerpimage(): erp limits args out of order.\n');
1524                 return
1525             end
1526             if length(Arg)> 4
1527                 minamp = Arg(5);
1528             end
1529             if length(Arg)> 5
1530                 maxamp = Arg(6);
1531             end
1532             if maxamp <= minamp
1533                 help erpimage
1534                 fprintf('\nerpimage(): amp limits args out of order.\n');
1535                 return
1536             end
1537             if length(Arg)> 6
1538                 mincoh = Arg(7);
1539             end
1540             if length(Arg)> 7
1541                 maxcoh = Arg(8);
1542             end
1543             if maxcoh <= mincoh
1544                 help erpimage
1545                 fprintf('\nerpimage(): coh limits args out of order.\n');
1546                 return
1547             end
1548             if length(Arg)>8
1549                 baseamp = Arg(9);    % for 'allamps'
1550             end
1551             Limitflag = NO;
1552 
1553         elseif Srateflag == YES
1554             srate = Arg(1);
1555             Srateflag = NO;
1556         elseif Cycleflag == YES
1557             cycles = Arg;
1558             Cycleflag = NO;
1559         elseif Auxvarflag == YES;
1560             if isa(Arg,'cell')==YES && length(Arg)==2
1561                 auxvar = Arg{1};
1562                 auxcolors = Arg{2};
1563             elseif isa(Arg,'cell')==YES
1564                 fprintf('erpimage(): auxvars argument must be a matrix or length-2 cell array.\n');
1565                 return
1566             else
1567                 auxvar = Arg; % no auxcolors specified
1568             end
1569             [xr,xc] = size(auxvar);
1570             lns = length(sortvar);
1571             if xr ~= lns && xc ~= lns
1572                error('auxvar columns different from the number of epochs in data');
1573             elseif xr == lns && xc ~= lns
1574                auxvar = auxvar';   % exchange rows/cols
1575             end
1576             Auxvarflag = NO;
1577         elseif Vertflag == YES
1578             verttimes = Arg;
1579             Vertflag = NO;
1580         elseif Horzflag == YES
1581             horzepochs = Arg;
1582             Horzflag = NO;
1583         elseif yerplabelflag == YES
1584             yerplabel = Arg;
1585             yerplabelflag = NO;
1586         elseif Signifflag == YES
1587             signifs = Arg; % [low_amp hi_amp coher]
1588             if length(signifs) ~= 3
1589                 fprintf('\nerpimage(): signif arg [%g] must have 3 values\n',Arg);
1590                 return
1591             end
1592             Signifflag = NO;
1593         elseif Allcohersflag == YES
1594             data2=Arg;
1595             if size(data2) ~= size(data)
1596                 fprintf('erpimage(): allcohers data matrix must be the same size as data.\n');
1597                 return
1598             end
1599             Allcohersflag = NO;
1600         elseif Phaseflag == YES
1601             n = length(Arg);
1602             if n > 5
1603                 error('erpimage(): Too many arguments for keyword ''phasesort''');
1604             end
1605             phargs = Arg;
1606 
1607             if phargs(3) < 0
1608                 error('erpimage(): Invalid negative frequency argument for keyword ''phasesort''');
1609             end
1610             if n>=4
1611                 if phargs(4) < 0
1612                     error('erpimage(): Invalid negative argument for keyword ''phasesort''');
1613                 end
1614             end
1615             if min(phargs(1)) < times(1) | max(phargs(1)) > times(end)
1616                 error('erpimage(): time for phase sorting filter out of bound.');
1617             end
1618             if phargs(2) >= 100 | phargs(2) < -100
1619                 error('%-argument for keyword ''phasesort'' must be (-100;100)');
1620             end
1621             if length(phargs) >= 4 & phargs(3) > phargs(4)
1622                 error('erpimage(): Phase sorting frequency range must be increasing.');
1623             end
1624             if length(phargs) == 5
1625                 topphase = phargs(5);
1626             end
1627             Phaseflag = NO;
1628         elseif Sortwinflag == YES % 'ampsort' mean amplitude over a time window
1629             n = length(Arg);
1630             sortwinarg = Arg;
1631             if n > 2
1632                 error('erpimage(): Too many arguments for keyword ''sortwin''');
1633             end
1634             if min(sortwinarg(1)) < times(1) | max(sortwinarg(1)) > times(end)
1635                 error('erpimage(): start time for value sorting out of bounds.');
1636             end
1637             if n > 1
1638                 if min(sortwinarg(2)) < times(1) | max(sortwinarg(2)) > times(end)
1639                     error('erpimage(): end time for value sorting out of bounds.');
1640                 end
1641             end
1642             if n > 1 & sortwinarg(1) > sortwinarg(2)
1643                 error('erpimage(): Value sorting time range must be increasing.');
1644             end
1645             Sortwinflag = NO;
1646         elseif Ampflag == YES % 'ampsort',[center_time,prcnt_reject,minfreq,maxfreq]
1647             n = length(Arg);
1648             if n > 4
1649                 error('erpimage(): Too many arguments for keyword ''ampsort''');
1650             end
1651             ampargs = Arg;
1652 
1653             % if ampargs(3) < 0
1654             %    error('erpimage(): Invalid negative argument for keyword ''ampsort''');
1655             % end
1656             if n>=4
1657                 if ampargs(4) < 0
1658                     error('erpimage(): Invalid negative argument for keyword ''ampsort''');
1659                 end
1660             end
1661 
1662             if ~isinf(ampargs(1))
1663                 if min(ampargs(1)) < times(1) | max(ampargs(1)) > times(end)
1664                     error('erpimage(): time for amplitude sorting filter out of bounds.');
1665                 end
1666             end
1667 
1668             if ampargs(2) >= 100 | ampargs(2) < -100
1669                 error('percentile argument for keyword ''ampsort'' must be (-100;100)');
1670             end
1671 
1672             if length(ampargs) == 4 & abs(ampargs(3)) > abs(ampargs(4))
1673                 error('erpimage(): Amplitude sorting frequency range must be increasing.');
1674             end
1675             Ampflag = NO;
1676 
1677         elseif Valflag == YES % sort by potential value in a given window
1678             % Usage: 'valsort',[mintime,maxtime,direction]
1679             n = length(Arg);
1680             if n > 3
1681                 error('erpimage(): Too many arguments for keyword ''valsort''');
1682             end
1683             valargs = Arg;
1684 
1685             if min(valargs(1)) < times(1) | max(valargs(1)) > times(end)
1686                 error('erpimage(): start time for value sorting out of bounds.');
1687             end
1688             if n > 1
1689                 if min(valargs(2)) < times(1) | max(valargs(2)) > times(end)
1690                     error('erpimage(): end time for value sorting out of bounds.');
1691                 end
1692             end
1693             if n > 1 & valargs(1) > valargs(2)
1694                 error('erpimage(): Value sorting time range must be increasing.');
1695             end
1696             if n==3 & (~isnumeric(valargs(3)) | valargs(3)==0)
1697                 error('erpimage(): Value sorting direction must be +1 or -1.');
1698             end
1699             Valflag = NO;
1700         elseif plotmodeflag == YES
1701             plotmode = Arg; plotmodeflag = NO;
1702         elseif titleflag == YES
1703             titl = Arg; titleflag = NO;
1704         elseif Erpalphaflag == YES
1705             erpalpha = Arg(1);
1706             if erpalpha < MIN_ERPALPHA | erpalpha > MAX_ERPALPHA
1707                 fprintf('erpimage(): erpalpha value is out of bounds [%g, %g]\n',...
1708                     MIN_ERPALPHA,MAX_ERPALPHA);
1709                 return
1710             end
1711             Erpalphaflag = NO;
1712         % -----------------------------------------------------------------------
1713         % -----------------------------------------------------------------------
1714         % -----------------------------------------------------------------------
1715        elseif strcmpi(Arg,'avg_type')
1716             if a < nargin, 
1717                 a=a+1;
1718                 Arg = eval(['arg' int2str(a-6)]); 
1719                 if strcmpi(Arg, 'Gaussian'), mvavg_type='gaussian'; 
1720                 elseif strcmpi(Arg, 'Boxcar'), mvavg_type='boxcar'; 
1721                 else error('Invalid value for optional argument ''avg_type''.');
1722                 end;
1723             else
1724                 error('Optional argument ''avg_type'' needs to be assigned a value.');
1725             end
1726         elseif strcmp(Arg,'nosort')
1727             Nosort = YES;
1728             if a < nargin, 
1729                 Arg = eval(['arg' int2str(a+1-6)]); 
1730                 if strcmpi(Arg, 'on'),     Nosort = YES; a = a+1;
1731                 elseif strcmpi(Arg, 'off') Nosort = NO;  a = a+1;
1732                 end;
1733             end;
1734         elseif strcmp(Arg,'showwin')
1735             Showwin = YES;
1736             if a < nargin, 
1737                 Arg = eval(['arg' int2str(a+1-6)]); 
1738                 if strcmpi(Arg, 'on'),     Showwin = YES; a = a+1;
1739                 elseif strcmpi(Arg, 'off') Showwin = NO;  a = a+1;
1740                 end;
1741             end;
1742         elseif strcmp(Arg,'noplot')
1743             Noshow = YES;
1744             if a < nargin, 
1745                 Arg = eval(['arg' int2str(a+1-6)]); 
1746                 if strcmpi(Arg, 'on'),     Noshow = YES; a = a+1;
1747                 elseif strcmpi(Arg, 'off') Noshow = NO;  a = a+1;
1748                 end;
1749             end;
1750         elseif strcmpi(Arg,'replace_ties')
1751             if a < nargin,
1752                 a = a+1;
1753                 temp = eval(['arg' int2str(a-6)]);
1754                 if strcmpi(temp,'on'),
1755                     replace_ties = YES;
1756                 elseif strcmpi(temp,'off') replace_ties = NO;
1757                 else
1758                     error('Argument ''replace_ties'' needs to be followed by the string ''on'' or ''off''.');
1759                 end
1760             else
1761                 error('Argument ''replace_ties'' needs to be followed by the string ''on'' or ''off''.');
1762             end
1763         elseif strcmpi(Arg,'sortvar_limits')
1764             if a < nargin,
1765                 a = a+1;
1766                 sortvar_limits = eval(['arg' int2str(a-6)]);
1767                 if ischar(sortvar_limits) || length(sortvar_limits)~=2
1768                   error('Argument ''sortvar_limits'' needs to be followed by a two element vector.');  
1769                 end
1770             else
1771                 error('Argument ''sortvar_limits'' needs to be followed by a two element vector.');
1772             end
1773         elseif strcmpi(Arg,'erp')
1774             Erpflag = YES;
1775             erp_ptiles=1;
1776             if a < nargin, 
1777                 Arg = eval(['arg' int2str(a+1-6)]); 
1778                 if strcmpi(Arg, 'on'),     Erpflag = YES; erp_ptiles=1; a = a+1;
1779                 elseif strcmpi(Arg, 'off') Erpflag = NO;  a = a+1;
1780                 elseif strcmpi(Arg,'1') | (Arg==1) Erplag = YES; erp_ptiles=1; a=a+1;
1781                 elseif strcmpi(Arg,'2') | (Arg==2) Erplag = YES; erp_ptiles=2; a=a+1;
1782                 elseif strcmpi(Arg,'3') | (Arg==3) Erplag = YES; erp_ptiles=3; a=a+1;
1783                 elseif strcmpi(Arg,'4') | (Arg==4) Erplag = YES; erp_ptiles=4; a=a+1;
1784                 end;
1785             end;            
1786         elseif strcmpi(Arg,'rmerp')
1787             Rmerp = 'yes';
1788             if a < nargin, 
1789                 Arg = eval(['arg' int2str(a+1-6)]); 
1790                 if strcmpi(Arg, 'on'),     Rmerp = 'yes'; a = a+1;
1791                 elseif strcmpi(Arg, 'off') Rmerp = 'no';  a = a+1;
1792                 end;
1793             end;            
1794         elseif strcmp(Arg,'cbar') | strcmp(Arg,'colorbar')
1795             Colorbar = YES;
1796             if a < nargin, 
1797                 Arg = eval(['arg' int2str(a+1-6)]); 
1798                 if strcmpi(Arg, 'on'),     Colorbar = YES; a = a+1;
1799                 elseif strcmpi(Arg, 'off') Colorbar = NO;  a = a+1;
1800                 end;
1801             end;            
1802         elseif (strcmp(Arg,'allamps') | strcmp(Arg,'plotamps'))
1803             Allampsflag = YES;
1804             if a < nargin, 
1805                 Arg = eval(['arg' int2str(a+1-6)]); 
1806                 if strcmpi(Arg, 'on'),     Allampsflag = YES; a = a+1;
1807                 elseif strcmpi(Arg, 'off') Allampsflag = NO;  a = a+1;
1808                 end;
1809             end;            
1810         elseif strcmpi(Arg,'erpstd')
1811             Erpstdflag = YES;
1812             if a < nargin, 
1813                 Arg = eval(['arg' int2str(a+1-6)]); 
1814                 if strcmpi(Arg, 'on'),     Erpstdflag = YES; a = a+1;
1815                 elseif strcmpi(Arg, 'off') Erpstdflag = NO;  a = a+1;
1816                 end;
1817             end;            
1818         elseif strcmp(Arg,'noxlabel') | strcmp(Arg,'noxlabels') | strcmp(Arg,'nox')
1819             NoTimeflag = YES;
1820             if a < nargin, 
1821                 Arg = eval(['arg' int2str(a+1-6)]); 
1822                 if strcmpi(Arg, 'on'),     NoTimeflag = YES; a = a+1;
1823                 elseif strcmpi(Arg, 'off') NoTimeflag = NO;  a = a+1;
1824                 end;
1825             end;            
1826         elseif strcmp(Arg,'plotmode')
1827             plotmodeflag = YES;
1828         elseif strcmp(Arg,'sortvarpercent')
1829             percentileflag = YES;
1830         elseif strcmp(Arg,'renorm')
1831             Renormflag = YES;
1832         elseif strcmp(Arg,'noshow')
1833             Noshowflag = YES;
1834         elseif strcmp(Arg,'caxis')
1835             Caxflag = YES;
1836         elseif strcmp(Arg,'title')
1837             titleflag = YES;
1838         elseif strcmp(Arg,'coher')
1839             Coherflag = YES;
1840         elseif strcmp(Arg,'timestretch') | strcmp(Arg,'timewarp') % Added -JH
1841             timestretchflag = YES;
1842         elseif strcmp(Arg,'allcohers')
1843             Allcohersflag = YES;
1844         elseif strcmp(Arg,'topo') | strcmp(Arg,'topoplot')
1845             Topoflag = YES;
1846         elseif strcmp(Arg,'spec') | strcmp(Arg,'spectrum')
1847             Specflag = YES;
1848         elseif strcmpi(Arg,'erpalpha')
1849             Erpalphaflag = YES;
1850         elseif strcmp(Arg,'align')
1851             Alignflag = YES;
1852         elseif strcmp(Arg,'limits')
1853             Limitflag = YES;
1854         elseif (strcmp(Arg,'phase') | strcmp(Arg,'phasesort'))
1855             Phaseflag = YES;
1856         elseif strcmp(Arg,'ampsort')
1857             Ampflag = YES;
1858         elseif strcmp(Arg,'sortwin')
1859             Sortwinflag = YES;
1860         elseif strcmp(Arg,'valsort')
1861             Valflag = YES;
1862         elseif strcmp(Arg,'auxvar')
1863             Auxvarflag = YES;
1864         elseif strcmp(Arg,'cycles')
1865             Cycleflag = YES;
1866         elseif strcmpi(Arg,'yerplabel')
1867             yerplabelflag = YES;
1868         elseif strcmpi(Arg,'srate')
1869             Srateflag = YES;
1870         elseif strcmpi(Arg,'erp_grid')
1871             erp_grid = YES;
1872         elseif strcmpi(Arg,'baseline')
1873             if a < nargin,
1874                 a = a+1;
1875                 baseline = eval(['arg' int2str(a-6)]);
1876             else
1877                 error('Argument ''baseline'' needs to be followed by a two element vector.');
1878             end
1879         elseif strcmpi(Arg,'filt')
1880             if a < nargin,
1881                 a = a+1;
1882                 flt = eval(['arg' int2str(a-6)]);
1883             else
1884                 error('Argument ''filt'' needs to be followed by a two element vector.');
1885             end
1886         elseif strcmpi(Arg,'erp_vltg_ticks')
1887             if a < nargin,
1888                 a = a+1;
1889                 erp_vltg_ticks=eval(['arg' int2str(a-6)]);
1890             else
1891                 error('Argument ''erp_vltg_ticks'' needs to be followed by a vector.');
1892             end          
1893         elseif strcmpi(Arg,'img_trialax_label')
1894             if a < nargin,
1895                 a = a+1;
1896                 img_ylab = eval(['arg' int2str(a-6)]);
1897             else
1898                 error('Argument ''img_trialax_label'' needs to be followed by a string.');
1899             end;
1900         elseif strcmpi(Arg,'img_trialax_ticks')
1901             if a < nargin,
1902                 a = a+1;
1903                 img_ytick_lab = eval(['arg' int2str(a-6)]);
1904             else
1905                 error('Argument ''img_trialax_ticks'' needs to be followed by a vector of values at which tick marks will appear.');
1906             end;
1907         elseif strcmpi(Arg,'cbar_title')
1908             if a < nargin, 
1909                 a = a+1;
1910                 cbar_title = eval(['arg' int2str(a-6)]); 
1911             else
1912                 error('Argument ''cbar_title'' needs to be followed by a string.');
1913             end;   
1914         elseif strcmp(Arg,'vert') ||  strcmp(Arg,'verttimes')
1915             Vertflag = YES;
1916         elseif strcmp(Arg,'horz') ||  strcmp(Arg,'horiz') || strcmp(Arg,'horizontal')
1917             Horzflag = YES;
1918         elseif strcmp(Arg,'signif') || strcmp(Arg,'signifs') || strcmp(Arg,'sig') || strcmp(Arg,'sigs')
1919             Signifflag = YES;
1920         else
1921             help erpimage
1922             if isstr(Arg)
1923                 fprintf('\nerpimage(): unknown arg %s\n',Arg);
1924             else
1925                 fprintf('\nerpimage(): unknown arg %d, size(%d,%d)\n',a,size(Arg,1),size(Arg,2));
1926             end
1927             return
1928         end
1929     end % Arg
1930 end
1931 
1932 if exist('img_ylab','var') || exist('img_ytick_lab','var'),
1933     oops=0;
1934     if exist('phargs','var'),
1935         fprintf('********* Warning *********\n');
1936         fprintf('Options ''img_ylab'' and ''img_ytick_lab'' have no effect when sorting by phase.\n');
1937         oops=0;
1938     elseif exist('valargs','var'),
1939         fprintf('********* Warning *********\n');
1940         fprintf('Options ''img_ylab'' and ''img_ytick_lab'' have no effect when sorting by EEG voltage.\n');
1941         oops=0;
1942     elseif exist('ampargs','var'),
1943         fprintf('********* Warning *********\n');
1944         fprintf('Options ''img_ylab'' and ''img_ytick_lab'' have no effect when sorting by frequency amplitude.\n');
1945         oops=0;
1946     end
1947     if oops
1948         img_ylab=[];
1949         img_ytick_lab=[];
1950     end
1951 end
1952 
1953 
1954 if strcmpi(noshow, 'off'), noshow = 'no'; end;
1955 
1956 if   Caxflag == YES ...
1957         |Coherflag == YES ...
1958         |Alignflag == YES ...
1959         |Limitflag == YES
1960     help erpimage
1961     fprintf('\nerpimage(): missing option arg.\n')
1962     return
1963 end
1964 if (Allampsflag | exist('data2')) & ( any(isnan(coherfreq)) | ~Cohsigflag )
1965     fprintf('\nerpimage(): allamps and allcohers flags require coher freq, srate, and cohsig.\n');
1966     return
1967 end
1968 if Allampsflag & exist('data2')
1969     fprintf('\nerpimage(): cannot image both allamps and allcohers.\n');
1970     return
1971 end
1972 if ~exist('srate') | srate <= 0
1973     fprintf('\nerpimage(): Data srate must be specified and > 0.\n');
1974     return
1975 end
1976 if ~isempty(auxvar)
1977     % whos auxvar
1978     if size(auxvar,1) ~= ntrials & size(auxvar,2) ~= ntrials
1979         fprintf('erpimage(): auxvar size should be (N,ntrials), e.g., (N,%d)\n',...
1980             ntrials);
1981         return
1982     end
1983     if size(auxvar,1) == ntrials & size(auxvar,2) ~= ntrials  % make (N,frames)
1984         auxvar = auxvar';
1985     end
1986     if size(auxvar,2) ~= ntrials
1987         fprintf('erpimage(): auxvar size should be (N,ntrials), e.g., (N,%d)\n',...
1988             ntrials);
1989         return
1990     end
1991     if exist('auxcolors')==YES % if specified
1992         if isa(auxcolors,'cell')==NO % if auxcolors is not a cell array
1993             fprintf(...
1994                 'erpimage(): auxcolors argument to auxvar flag must be a cell array.\n');
1995             return
1996         end
1997     end
1998 elseif exist('timeStretchRef') & ~isempty(timeStretchRef)
1999     if ~isnan(aligntime)
2000         fprintf(['erpimage(): options "align" and ' ...
2001             '"timewarp" are not compatiable.\n']);
2002         return;
2003     end
2004 
2005     if ~isempty(timeStretchColors)
2006         if length(timeStretchColors) < length(timeStretchRef)
2007             nColors = length(timeStretchColors);
2008             for k=nColors+1:length(timeStretchRef)-2
2009                 timeStretchColors = { timeStretchColors{:} timeStretchColors{1+rem(k-1,nColors)}};
2010             end
2011         end
2012         timeStretchColors = {'' timeStretchColors{:} ''};
2013     else
2014         timeStretchColors = { 'k--'};
2015         for k=2:length(timeStretchRef)
2016             timeStretchColors = { timeStretchColors{:} 'k--'};
2017         end
2018     end
2019 
2020 
2021     auxvarInd = 1-strcmp('',timeStretchColors); % indicate which lines to draw
2022     newauxvars = ((timeStretchRef(find(auxvarInd))-1)/srate+times(1)/1000) * 1000; % convert back to ms
2023     fprintf('Overwriting vert with auxvar\n');
2024     verttimes = [newauxvars'];
2025     verttimesColors = {timeStretchColors{find(auxvarInd)}};
2026     newauxvars = repmat(newauxvars, [1 ntrials]);
2027 
2028     if isempty(auxvar) % Initialize auxvar & auxcolors
2029         %  auxvar = newauxvars;
2030         auxcolors = {timeStretchColors{find(auxvarInd)}};
2031     else % Append auxvar & auxcolors
2032         if ~exist('auxcolors')
2033             % Fill with default color (k-- for now)
2034             auxcolors = {};
2035             for j=1:size(auxvar,1)
2036                 auxcolors{end+1} = 'k--';
2037             end
2038         end
2039         for j=find(auxvarInd)
2040             auxcolors{end+1} = timeStretchColors{j};
2041         end
2042         auxvar = [auxvar; newauxvars];
2043     end
2044 end
2045 
2046 % No need to turn off ERP anymore; we now time-stretch potentials
2047 % separately (see tsurdata below) so the (warped) ERP is actually meaningful
2048 
2049 if exist('phargs')
2050     if phargs(3) > srate/2
2051         fprintf(...
2052             'erpimage(): Phase-sorting frequency (%g Hz) must be less than Nyquist rate (%g Hz).',...
2053             phargs(3),srate/2);
2054     end
2055 
2056     if frames < cycles*srate/phargs(3)
2057         fprintf('\nerpimage(): phase-sorting freq. (%g) too low: epoch length < %d cycles.\n',...
2058             phargs(3),cycles);
2059         return
2060     end
2061     if length(phargs)==4 & phargs(4) > srate/2
2062         phargs(4) = srate/2;
2063     end
2064     if length(phargs)==5 & (phargs(5)>180 | phargs(5) < -180)
2065         fprintf('\nerpimage(): coher topphase (%g) out of range.\n',topphase);
2066         return
2067     end
2068 end
2069 if exist('ampargs')
2070     if abs(ampargs(3)) > srate/2
2071         fprintf(...
2072             'erpimage(): amplitude-sorting frequency (%g Hz) must be less than Nyquist rate (%g Hz).',...
2073             abs(ampargs(3)),srate/2);
2074     end
2075 
2076     if frames < cycles*srate/abs(ampargs(3))
2077         fprintf('\nerpimage(): amplitude-sorting freq. (%g) too low: epoch length < %d cycles.\n',...
2078             abs(ampargs(3)),cycles);
2079         return
2080     end
2081     if length(ampargs)==4 & abs(ampargs(4)) > srate/2
2082         ampargs(4) = srate/2;
2083         fprintf('> Reducing max ''ampsort'' frequency to Nyquist rate (%g Hz)\n',srate/2)
2084     end
2085 end
2086 if ~any(isnan(coherfreq))
2087     if coherfreq(1) <= 0 | srate <= 0
2088         fprintf('\nerpimage(): coher frequency (%g) out of range.\n',coherfreq(1));
2089         return
2090     end
2091     if coherfreq(end) > srate/2 | srate <= 0
2092         fprintf('\nerpimage(): coher frequency (%g) out of range.\n',coherfreq(end));
2093         return
2094     end
2095     if frames < cycles*srate/coherfreq(1)
2096         fprintf('\nerpimage(): coher freq. (%g) too low:  epoch length < %d cycles.\n',...
2097             coherfreq(1),cycles);
2098         return
2099     end
2100 end
2101 
2102 if isnan(timelimits)
2103     timelimits = [min(times) max(times)];
2104 end
2105 if ~isstr(aligntime) & ~isnan(aligntime)
2106     if ~isinf(aligntime) ...
2107             & (aligntime < timelimits(1) | aligntime > timelimits(2))
2108         help erpimage
2109         fprintf('\nerpimage(): requested align time outside of time limits.\n');
2110         return
2111     end
2112 end
2113 %
2114 %% %%%%%%%%%%%%%%  Replace nan's with 0s %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2115 %
2116 nans = find(isnan(data));
2117 if length(nans)
2118     fprintf('Replaced %d nan in data with 0s.\n');
2119     data(nans) = 0;
2120 end
2121 %
2122 %% %%%%%%%%%%%% Reshape data to (frames,ntrials) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2123 %
2124 if size(data,2) ~= ntrials
2125     if size(data,1)>1
2126         % fprintf('frames %d, ntrials %d length(data) %d\n',frames,ntrials,length(data));
2127         data=reshape(data,1,frames*ntrials);
2128     end
2129     data=reshape(data,frames,ntrials);
2130 end
2131 fprintf('Plotting input data as %d epochs of %d frames sampled at %3.1f Hz.\n',...
2132     ntrials,frames,srate);
2133 %
2134 %% %%%%%%%%%%%% Reshape data2 to (frames,ntrials) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2135 %
2136 if exist('data2') == 1
2137     if size(data2,2) ~= ntrials
2138         if size(data2,1)>1
2139             data2=reshape(data2,1,frames*ntrials);
2140         end
2141         data2=reshape(data2,frames,ntrials);
2142     end
2143 end
2144 %
2145 %% %%%%%%%%%%%%% if sortvar=NaN, remove lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ad
2146 %
2147 if any(isnan(sortvar))
2148     nanlocs = find(isnan(sortvar));
2149     fprintf('Removing %d trials with NaN sortvar values.\n', length(nanlocs));
2150     data(:,nanlocs) = [];
2151     sortvar(nanlocs) = [];
2152     if exist('data2') == 1
2153         data2(:,nanlocs) = [];
2154     end;
2155     if ~isempty(auxvar)
2156         auxvar(:,nanlocs) = [];
2157     end
2158     if ~isempty(verttimes)
2159         if size(verttimes,1) == ntrials
2160             verttimes(nanlocs,:) = [];
2161         end;
2162     end;
2163     ntrials = size(data,2);
2164     if ntrials <= 1, close(gcf); error('Too few trials'); end;
2165 end;
2166 
2167 %% Create moving average window %%
2168 if strcmpi(mvavg_type,'Gaussian'),
2169     %construct Gaussian window to weight trials
2170     if avewidth == 0,
2171         avewidth = DEFAULT_SDEV;
2172     elseif avewidth < 1,
2173         help erpimage
2174         fprintf('\nerpimage(): Variable avewidth cannot be < 1.\n')
2175         fprintf('\nerpimage(): avewidth needs to be a positive integer.\n')
2176         return
2177     end
2178     wt_wind=exp(-0.5*([-3*avewidth:3*avewidth]/avewidth).^2)';
2179     wt_wind=wt_wind/sum(wt_wind); %normalize to unit sum
2180     avewidth=length(wt_wind);
2181     
2182     if avewidth > ntrials
2183         avewidth=floor((ntrials-1)/6);
2184         if avewidth==0,
2185             avewidth=DEFAULT_SDEV; %should be a window with one time point (smallest possible)
2186         end
2187         wt_wind=exp(-0.5*([-3*avewidth:3*avewidth]/avewidth).^2)';
2188         wt_wind=wt_wind/sum(wt_wind);
2189         fprintf('avewidth is too big for this number of trials.\n');
2190         fprintf('Changing avewidth to maximum possible size: %d\n',avewidth);
2191         avewidth=length(wt_wind);
2192     end
2193 else
2194     %construct rectangular "boxcar" window to equally weight trials within
2195     %window
2196     if avewidth == 0,
2197         avewidth = DEFAULT_AVEWIDTH;
2198     elseif avewidth < 1
2199         help erpimage
2200         fprintf('\nerpimage(): Variable avewidth cannot be < 1.\n')
2201         return
2202     elseif avewidth > ntrials
2203         fprintf('Setting variable avewidth to max %d.\n',ntrials)
2204         avewidth = ntrials;
2205     end
2206     wt_wind=ones(1,avewidth)/avewidth; 
2207 end
2208 
2209 
2210 %% Filter data with Butterworth filter (if requested) %%%%%%%%%%%%
2211 %
2212 if ~isempty(flt)
2213     %error check
2214     if length(flt)~=2,
2215        error('''filt'' parameter argument should be a two element vector.');
2216     elseif max(flt)>(srate/2),
2217        error('''filt'' parameters need to be less than or equal to sampling rate/2 (i.e., %f).',srate/2);
2218     elseif (flt(2)==(srate/2)) && (flt(1)==0),
2219         error('If second element of ''filt'' parameter is srate/2, then the first element must be greater than 0.');
2220     elseif abs(flt(2))<=abs(flt(1)),
2221        error('Second element of ''filt'' parameters must be greater than first in absolute value.'); 
2222     elseif (flt(1)<0) || (flt(2)<0),
2223         if (flt(1)>=0) || (flt(2)>=0),
2224            error('BOTH parameters of ''filt'' need to be greater than or equal to zero OR need to be negative.');
2225         end
2226         if min(flt)<=(-srate/2),
2227             error('''filt'' parameters need to be greater than sampling rate/2 (i.e., -%f) when creating a stop band.',srate/2);
2228         end
2229     end
2230     
2231     fprintf('\nFiltering data with 3rd order Butterworth filter: ');
2232     if (flt(1)==0),
2233         %lowpass filter the data
2234         [B A]=butter(3,flt(2)*2/srate,'low');
2235         fprintf('lowpass at %.0f Hz\n',flt(2));
2236     elseif (flt(2)==(srate/2)),
2237         %highpass filter the data
2238         [B A]=butter(3,flt(1)*2/srate,'high');
2239         fprintf('highpass at %.0f Hz\n',flt(1));
2240     elseif (flt(1)<0)
2241         %bandstop filter the data
2242         flt=-flt;
2243         [B A]=butter(3,flt*2/srate,'stop');
2244         fprintf('stopband from %.0f to %.0f Hz\n',flt(1),flt(2));
2245     else
2246         %bandpass filter the data
2247         [B A]=butter(3,flt*2/srate);
2248         fprintf('bandpass from %.0f to %.0f Hz\n',flt(1),flt(2));
2249     end
2250     s=size(data);
2251     for trial=1:s(2),
2252         data(:,trial)=filtfilt(B,A,data(:,trial));
2253     end
2254     if isempty(baseline)
2255         fprintf('Note, you might want to re-baseline the data using the erpiamge ''baseline'' option.\n\n');
2256     end
2257 end
2258 
2259 %% Mean Baseline Each Trial (if requested) %%
2260 if ~isempty(baseline),
2261     %check argument values for errors
2262     if baseline(2)<baseline(1),
2263        error('First element of ''baseline'' argument needs to be less than or equal to second argument.'); 
2264     elseif baseline(2)<times(1),
2265         error('Second element of ''baseline'' argument needs to be greater than or equal to epoch start time %.1f.',times(1));
2266     elseif baseline(1)>times(end),
2267         error('First element of ''baseline'' argument needs to be less than or equal to epoch end time %.1f.',times(end));
2268     end
2269     
2270     %convert msec into time points
2271     if baseline(1)<times(1),
2272         strt_pt=1;
2273     else
2274         strt_pt=find_crspnd_pt(baseline(1),times,1:length(times));
2275         strt_pt=ceil(strt_pt);
2276     end
2277     if baseline(end)>times(end),
2278         end_pt=length(times);
2279     else
2280         end_pt=find_crspnd_pt(baseline(2),times,1:length(times));
2281         end_pt=floor(end_pt);
2282     end
2283     fprintf('\nRemoving pre-stimulus mean baseline from %.1f to %.1f msec.\n\n',times(strt_pt),times(end_pt));
2284     bsln_mn=mean(data(strt_pt:end_pt,:),1);
2285     data=data-repmat(bsln_mn,length(times),1);
2286 end
2287 
2288 
2289 
2290 %
2291 %% %%%%%%%%%%%%%%%%% Renormalize sortvar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2292 %
2293 switch lower(renorm)
2294     case 'yes',
2295         disp('erpimage warning: *** sorting variable renormalized ***');
2296         sortvar = (sortvar-min(sortvar)) / (max(sortvar) - min(sortvar)) * ...
2297             0.5 * (max(times) - min(times)) + min(times) + 0.4*(max(times) - min(times));
2298     case 'no',;
2299     otherwise,
2300         if ~isempty(renorm)
2301             locx = findstr('x', lower(renorm));
2302             if length(locx) ~= 1, error('erpimage: unrecognized renormalizing formula'); end;
2303             eval( [ 'sortvar =' renorm(1:locx-1) 'sortvar' renorm(locx+1:end) ';'] );
2304         end;
2305 end;
2306 %
2307 %% %%%%%%%%%%%%%%%%% Align data to sortvar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2308 %
2309 
2310 if isstr(aligntime) | ~isnan(aligntime)
2311     if ~isstr(aligntime) & isinf(aligntime)
2312         aligntime= median(sortvar);
2313         fprintf('Aligning data to median sortvar.\n');
2314         % Alternative below: trimmed median - ignore top/bottom 5%
2315         %   ssv = sort(sortvar); % ssv = 'sorted sortvar'
2316         %   aligntime= median(ssv(ceil(ntrials/20)):floor(19*ntrials/20));
2317     end
2318 
2319     if ~isstr(aligntime)
2320         fprintf('Realigned sortvar plotted at %g ms.\n',aligntime);
2321         aligndata=zeros(frames,ntrials); % begin with matrix of zeros()
2322         shifts = zeros(1,ntrials);
2323         for t=1:ntrials, %%%%%%%%% foreach trial %%%%%%%%%
2324             shft = round((aligntime-sortvar(t))*srate/1000);
2325             shifts(t) = shft;
2326             if shft>0, % shift right
2327                 if frames-shft > 0
2328                     aligndata(shft+1:frames,t)=data(1:frames-shft,t);
2329                 else
2330                     fprintf('No aligned data for epoch %d - shift (%d frames) too large.\n',t,shft);
2331                 end
2332             elseif shft < 0 % shift left
2333                 if frames+shft > 0
2334                     aligndata(1:frames+shft,t)=data(1-shft:frames,t);
2335                 else
2336                     fprintf('No aligned data for epoch %d - shift (%d frames) too large.\n',t,shft);
2337                 end
2338             else % shft == 0
2339                 aligndata(:,t) = data(:,t);
2340             end
2341         end % end trial
2342         if ~isempty(auxvar)
2343             auxvar = auxvar+shifts;
2344         end;
2345         fprintf('Shifted epochs by %d to %d frames.\n',min(shifts),max(shifts));
2346         data = aligndata;                       % now data is aligned to sortvar
2347     else
2348         aligntime = str2num(aligntime);
2349         if isinf(aligntime),  aligntime= median(sortvar); end;
2350     end;
2351 end
2352 
2353 %
2354 %% %%%%%%%%%%%%%%%%%%%% Remove the ERP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2355 %
2356 if strcmpi(Rmerp, 'yes')
2357     data = data - repmat(nan_mean(data')', [1 size(data,2)]);
2358 end;
2359 
2360 %
2361 %% %%%%%%%%%%%%% Sort the data trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2362 %
2363 if exist('phargs') == 1 % if phase-sort the data trials
2364     if length(phargs) >= 4 & phargs(3) ~= phargs(4) % find max frequency
2365         % in specified band
2366         if exist('psd') == 2 % requires Signal Processing Toolbox
2367             fprintf('Computing data spectrum using psd().\n');
2368             [pxx,freqs] = psd(data(:),max(1024, pow2(ceil(log2(frames)))),srate,frames,0);
2369         else % EEGLAB native work-around
2370             fprintf('Computing data spectrum using spec().\n');
2371             [pxx,freqs] = spec(data(:),max(1024, pow2(ceil(log2(frames)))),srate,frames,0);
2372         end;
2373         % gf = gcf; % figure;plot(freqs,pxx); %xx=axis; %axis([phargs(3) phargs(4) xx(3) xx(4)]); %figure(gf);
2374         pxx = 10*log10(pxx);
2375         n = find(freqs >= phargs(3) & freqs <= phargs(4));
2376         if ~length(n)
2377             freq = (phargs(3)+phargs(4))/2;
2378         end
2379         [dummy maxx] = max(pxx(n));
2380         freq = freqs(n(maxx));
2381     else
2382         freq = phargs(3); % else use specified frequency
2383     end
2384     fprintf('Sorting trials on phase at %.2g Hz.\n',freq);
2385 
2386     [amps, cohers, cohsig, ampsig, allamps, allphs] = ...
2387         phasecoher(data,length(times),srate,freq,cycles,0, ...
2388         [], [], timeStretchRef, timeStretchMarks);
2389 
2390     phwin = phargs(1);
2391     [dummy minx] = min(abs(times-phwin)); % closest time to requested
2392     winlen = floor(cycles*srate/freq);
2393     winloc = minx-linspace(floor(winlen/2), floor(-winlen/2), winlen+1);
2394     tmprange = find(winloc>0 & winloc<=frames);
2395     winloc = winloc(tmprange); % sorting window times
2396     winlocs = winloc;
2397     %
2398     %%%%%%%%%%%%%%%%%%%% Compute phsamp and phaseangles %%%%%%%%%%%%%%%%%%%%
2399     %
2400     % $$$     [phaseangles phsamp] = phasedet(data,frames,srate,winloc,freq);
2401     phaseangles = allphs(minx,:);
2402     phsamp = allamps(minx,:);
2403     % $$$     % hist(phaseangles,50);
2404     % $$$     % return
2405     % $$$
2406     % $$$     %
2407     % $$$     % Print facts on commandline %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2408     % $$$     %
2409     % $$$     if length(tmprange) ~=  winlen+1
2410     % $$$       filtersize = cycles * length(tmprange) / (winlen+1);
2411     % $$$       timecenter = median(winloc)/srate*1000+times(1); % center of window in ms
2412     % $$$       phaseangles = phaseangles + 2*pi*(timecenter-phargs(1))*freq;
2413     % $$$       fprintf('Sorting data epochs by phase at frequency %2.1f Hz: \n', freq);
2414     % $$$       fprintf(...
2415     % $$$           '    Data time limits reached -> now uses a %1.1f cycles (%1.0f ms) window centered at %1.0f ms\n', ...
2416     % $$$           filtersize, 1000/freq*filtersize, timecenter);
2417     % $$$       fprintf(...
2418     % $$$           '    Filter length is %d; Phase has been linearly interpolated to latency at %1.0f ms.\n', ...
2419     % $$$           length(winloc), phargs(1));
2420     % $$$     else
2421     % $$$       fprintf(...
2422     % $$$           'Sorting data epochs by phase at %2.1f Hz in a %1.1f-cycle (%1.0f ms) window centered at %1.0f ms.\n',...
2423     % $$$           freq,cycles,1000/freq*cycles,times(minx));
2424     % $$$       fprintf('Phase is computed using a wavelet of %d frames.\n',length(winloc));
2425     % $$$     end;
2426     %
2427     % Reject small (or large) phsamp trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2428     %
2429     phargs(2) = phargs(2)/100; % convert rejection rate from % to fraction
2430     amprej = phargs(2);
2431     [tmp ampsortidx] = sort(phsamp); % sort amplitudes
2432     if amprej>=0
2433         ampsortidx = ampsortidx(ceil(amprej*length(ampsortidx))+1:end); % if amprej==0, select all trials
2434         fprintf('Retaining %d epochs (%g percent) with largest power at the analysis frequency,\n',...
2435             length(ampsortidx),100*(1-amprej));
2436     else % amprej < 0
2437         amprej = 1+amprej; % subtract from end
2438         ampsortidx = ampsortidx(1:floor(amprej*length(ampsortidx)));
2439         fprintf('Retaining %d epochs (%g percent) with smallest power at the analysis frequency,\n',...
2440             length(ampsortidx),amprej*100);
2441     end
2442     %
2443     % Remove low|high-amplitude trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2444     %
2445     data = data(:,ampsortidx); % amp-sort the data, removing rejected-amp trials
2446     phsamp = phsamp(ampsortidx);           % amp-sort the amps
2447     phaseangles = phaseangles(ampsortidx); % amp-sort the phaseangles
2448     sortvar = sortvar(ampsortidx);         % amp-sort the trial indices
2449     ntrials = length(ampsortidx);          % number of trials retained
2450     if ~isempty(auxvar)
2451         auxvar = auxvar(:,ampsortidx);
2452     end
2453     if ~isempty(timeStretchMarks)
2454         timeStretchMarks =  timeStretchMarks(:,ampsortidx);
2455     end
2456 
2457     %
2458     % Sort remaining data by phase angle %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2459     %
2460     phaseangles = -phaseangles;
2461     topphase = (topphase/360)*2*pi; % convert from degrees to radians
2462     ip = find(phaseangles>topphase);
2463     phaseangles(ip) = phaseangles(ip)-2*pi; % rotate so topphase at top of plot
2464 
2465     [phaseangles sortidx] = sort(phaseangles); % sort trials on (rotated) phase
2466     data    =  data(:,sortidx);                % sort data by phase
2467     phsamp  =  phsamp(sortidx);                % sort amps by phase
2468     sortvar = sortvar(sortidx);                % sort input sortvar by phase
2469     if ~isempty(auxvar)
2470         auxvar = auxvar(:,sortidx);
2471     end
2472     if ~isempty(timeStretchMarks)
2473         timeStretchMarks =  timeStretchMarks(:,sortidx);
2474     end
2475     phaseangles = -phaseangles; % Note: phsangles now descend from pi
2476     % TEST auxvar = 360 + (1000/256)*(256/5)*phaseangles/(2*pi); % plot phase+360 in ms for test
2477 
2478     fprintf('Size of data = [%d,%d]\n',size(data,1),size(data,2));
2479     sortidx = ampsortidx(sortidx); % return original trial indices in final sorted order
2480     %
2481     % %%%%%%%%%%%%%%% Sort data by amplitude %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2482     %
2483 elseif exist('ampargs') == 1 % if amplitude-sort
2484     if length(ampargs) == 4 % find max frequency in specified band
2485         if exist('psd') == 2
2486             fprintf('Computing data spectrum using psd().\n');
2487             [pxx,freqs] = psd(data(:),max(1024, pow2(ceil(log2(frames)))),srate,frames,0);
2488         else
2489             fprintf('Computing data spectrum using spec().\n');
2490             [pxx,freqs] = spec(data(:),max(1024, pow2(ceil(log2(frames)))),srate,frames,0);
2491         end;
2492         pxx = 10*log10(pxx);
2493         if ampargs(3) == ampargs(4)
2494             [freq n] = min(abs(freqs - ampargs(3)));
2495         else
2496             n = find(freqs >= abs(ampargs(3)) & freqs <= abs(ampargs(4)));
2497         end;
2498         if ~length(n)
2499             freq = mean([abs(ampargs(3)),abs(ampargs(4))]);
2500         end
2501         if ampargs(3)>=0
2502             [dummy maxx] = max(pxx(n));
2503             freq = freqs(n(maxx)); % use the highest-power frequency
2504         else
2505             freq = freqs(n);  % use all frequencies in the specified range
2506         end
2507     else
2508         freq = abs(ampargs(3)); % else use specified frequency
2509     end
2510     if length(freq) == 1
2511         fprintf('Sorting data epochs by amplitude at frequency %2.1f Hz \n', freq);
2512     else
2513         fprintf('Sorting data epochs by amplitude at %d frequencies (%2.1f Hz to %.1f Hz) \n',...
2514             length(freq),freq(1),freq(end));
2515     end
2516     SPECWININCR = 10;    % make spectral sorting time windows increment by 10 ms
2517     if isinf(ampargs(1))
2518         ampwins = sortwinarg(1):SPECWININCR:sortwinarg(2);
2519     else
2520         ampwins = ampargs(1);
2521     end
2522     if ~isinf(ampargs(1)) % single time given
2523         if length(freq) == 1
2524             fprintf('   in a %1.1f-cycle (%1.0f ms) time window centered at %1.0f ms.\n',...
2525                 cycles,1000/freq(1)*cycles,ampargs(1));
2526         else
2527             fprintf('   in %1.1f-cycle (%1.0f-%1.0f ms) time windows centered at %1.0f ms.\n',...
2528                 cycles,1000/freq(1)*cycles,1000/freq(end)*cycles,ampargs(1));
2529         end
2530     else % range of times
2531         [dummy sortwin_st ] = min(abs(times-ampwins(1)));
2532         [dummy sortwin_end] = min(abs(times-ampwins(end)));
2533         if length(freq) == 1
2534             fprintf('   in %d %1.1f-cycle (%1.0f ms) time windows centered from %1.0f to  %1.0f ms.\n',...
2535                 length(ampwins),cycles,1000/freq(1)*cycles,times(sortwin_st),times(sortwin_end));
2536         else
2537             fprintf('   in %d %1.1f-cycle (%1.0f-%1.0f ms) time windows centered from %1.0f to %1.0f ms.\n',...
2538                 length(ampwins),cycles,1000/freq(1)*cycles,1000/freq(end)*cycles,times(sortwin_st),times(sortwin_end));
2539         end
2540     end
2541 
2542     phsamps = 0; %%%%%%%%%%%%%%%%%%%%%%%%%% sort by (mean) amplitude %%%%%%%%%%%%%%%%%%%%%%%%%%
2543     minxs = [];
2544     for f = 1:length(freq)  % use one or range of frequencies
2545         frq = freq(f);
2546         [amps, cohers, cohsig, ampsig, allamps, allphs] = ...
2547             phasecoher(data,length(times),srate,frq,cycles,0, ...
2548             [], [], timeStretchRef, timeStretchMarks);
2549 
2550         for ampwin = ampwins
2551             [dummy minx] = min(abs(times-ampwin)); % find nearest time point to requested
2552             minxs = [minxs minx];
2553             winlen = floor(cycles*srate/frq);
2554             % winloc = minx-[winlen:-1:0]; % ending time version
2555             winloc = minx-linspace(floor(winlen/2), floor(-winlen/2), winlen+1);
2556             tmprange = find(winloc>0 & winloc<=frames);
2557             winloc = winloc(tmprange); % sorting window frames
2558             if f==1
2559                 winlocs = [winlocs;winloc];  % store tme windows
2560             end
2561             % $$$         [phaseangles phsamp] = phasedet(data,frames,srate,winloc,frq);
2562             % $$$         phsamps = phsamps+phsamps;  % accumulate amplitudes across 'sortwin'
2563             phaseangles = allphs(minx,:);
2564             phsamps = phsamps+allamps(minx,:);  % accumulate amplitudes across 'sortwin'
2565         end
2566     end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2567     if length(tmprange) ~=  winlen+1 % ?????????
2568         filtersize = cycles * length(tmprange) / (winlen+1);
2569         timecenter = median(winloc)/srate*1000+times(1); % center of window in ms
2570         phaseangles = phaseangles + 2*pi*(timecenter-ampargs(1))*freq(end);
2571         fprintf(...
2572             '    Data time limits reached -> now uses a %1.1f cycles (%1.0f ms) window centered at %1.0f ms\n', ...
2573             filtersize, 1000/freq(1)*filtersize, timecenter);
2574         fprintf(...
2575             '    Wavelet length is %d; Phase has been linearly interpolated to latency et %1.0f ms.\n', ...
2576             length(winloc(1,:)), ampargs(1));
2577     end
2578     if length(freq) == 1
2579         fprintf('Amplitudes are computed using a wavelet of %d frames.\n',length(winloc(1,:)));
2580     end
2581 
2582     %
2583     % Reject small (or large) phsamp trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2584     %
2585     ampargs(2) = ampargs(2)/100; % convert rejection rate from % to fraction
2586     [tmp n] = sort(phsamps); % sort amplitudes
2587     if ampargs(2)>=0
2588         n = n(ceil(ampargs(2)*length(n))+1:end); % if rej 0, select all trials
2589         fprintf('Retaining %d epochs (%g percent) with largest power at the analysis frequency,\n',...
2590             length(n),100*(1-ampargs(2)));
2591     else % ampargs(2) < 0
2592         ampargs(2) = 1+ampargs(2); % subtract from end
2593         n = n(1:floor(ampargs(2)*length(n)));
2594         fprintf(...
2595             'Retaining %d epochs (%g percent) with smallest power at the analysis frequency,\n',...
2596             length(n),ampargs(2)*100);
2597     end
2598     %
2599     % Remove low|high-amplitude trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2600     %
2601 
2602     data = data(:,n); % amp-sort the data, removing rejected-amp trials
2603     phsamps = phsamps(n);           % amp-sort the amps
2604     phaseangles = phaseangles(n); % amp-sort the phaseangles
2605     sortvar = sortvar(n);         % amp-sort the trial indices
2606     ntrials = length(n);          % number of trials retained
2607     if ~isempty(auxvar)
2608         auxvar = auxvar(:,n);
2609     end
2610     %
2611     % Sort remaining data by amplitude %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2612     %
2613     [phsamps sortidx] = sort(phsamps); % sort trials on amplitude
2614     data    =  data(:,sortidx);                % sort data by amp
2615     phaseangles  =  phaseangles(sortidx);      % sort angles by amp
2616     sortvar = sortvar(sortidx);                % sort input sortvar by amp
2617     if ~isempty(auxvar)
2618         auxvar = auxvar(:,sortidx);
2619     end
2620     fprintf('Size of data = [%d,%d]\n',size(data,1),size(data,2));
2621 
2622     sortidx = n(sortidx); % return original trial indices in final sorted order
2623     %
2624     %%%%%%%%%%%%%%%%%%%%%% Don't Sort trials %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2625     %
2626 elseif Nosort == YES
2627     fprintf('Not sorting data on input sortvar.\n');
2628     sortidx = 1:ntrials;
2629     %
2630     %%%%%%%%%%%%%%%%%%%%%% Sort trials on (mean) value %%%%%%%%%%%%%%%%%%%%%%%%%%%%
2631     %
2632 elseif exist('valargs')
2633     [sttime stframe] = min(abs(times-valargs(1)));
2634     sttime = times(stframe);
2635     if length(valargs)>1
2636         [endtime endframe] = min(abs(times-valargs(2)));
2637         endtime = times(endframe);
2638     else
2639         endframe = stframe;
2640         endtime = times(endframe);
2641     end
2642     if length(valargs)==1 || sttime == endtime
2643         fprintf('Sorting data on value at time %4.0f ms.\n',sttime);
2644     elseif length(valargs)>1
2645         fprintf('Sorting data on mean value between %4.0f and %4.0f ms.\n',...
2646             sttime,endtime);
2647     end
2648     if endframe>stframe
2649         sortval = mean(data(stframe:endframe,:));
2650     else
2651         sortval = data(stframe,:);
2652     end
2653     [sortval,sortidx] = sort(sortval);
2654     if length(valargs)>2
2655         if valargs(3) <0
2656             sortidx = sortidx(end:-1:1); % plot largest values on top
2657             % if direction < 0
2658         end
2659     end
2660     data = data(:,sortidx);
2661     sortvar = sortvar(sortidx);                % sort input sortvar by amp
2662     if ~isempty(auxvar)
2663         auxvar = auxvar(:,sortidx);
2664     end
2665     if ~isempty(phaseangles)
2666         phaseangles  =  phaseangles(sortidx);      % sort angles by amp
2667     end
2668     winloc = [stframe,endframe];
2669     winlocs = winloc;
2670     %
2671     %%%%%%%%%%%%%%%%%%%%%% Sort trials on sortvar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2672     %
2673 else
2674     fprintf('Sorting data on input sortvar.\n');
2675     [sortvar,sortidx] = sort(sortvar);
2676     data = data(:,sortidx);
2677     if ~isempty(auxvar)
2678         auxvar = auxvar(:,sortidx);
2679     end
2680     uni_svar=unique(sortvar);
2681     n_ties=0;
2682     tie_dist=zeros(1,length(uni_svar));
2683     loop_ct=0;
2684     for tie_loop=uni_svar,
2685         ids=find(sortvar==tie_loop);
2686         n_ids=length(ids);
2687         if n_ids>1,
2688             if replace_ties==YES,
2689                 mn=mean(data(:,ids),2);
2690                 data(:,ids)=repmat(mn,1,n_ids);
2691                 if ~isempty(auxvar)
2692                     mn=mean(auxvar(:,ids),2);
2693                     auxvar(:,ids) = repmat(mn,1,n_ids);
2694                 end
2695             end
2696             n_ties=n_ties+n_ids;
2697             loop_ct=loop_ct+1;
2698             tie_dist(loop_ct)=n_ids;
2699         end
2700     end
2701     fprintf('%.2f%c of the trials (i.e., %d out of %d) have the same sortvar value as at least one other trial.\n', ...
2702         100*n_ties/length(sortvar),37,n_ties,length(sortvar));
2703     fprintf('Distribution of number ties per unique value of sortvar:\n');
2704     fprintf('Min: %d, 25th ptile: %d, Median: %d, 75th ptile: %d, Max: %d\n',min(tie_dist),round(prctile(tie_dist,25)), ...
2705         round(median(tie_dist)),round(prctile(tie_dist,75)),max(tie_dist));
2706     if replace_ties==YES,
2707         fprintf('Trials with tied sorting values will be replaced by their mean.\n');
2708     end
2709     fprintf('\n');
2710 end
2711 
2712 %if max(sortvar)<0
2713 %   fprintf('Changing the sign of sortvar: making it positive.\n');
2714 %   sortvar = -sortvar;
2715 %end
2716 %
2717 %% %%%%%%%%%%%%%%%%% Adjust decfactor if phargs or ampargs %%%%%%%%%%%%%%%%%%%%%
2718 %
2719 if decfactor < 0
2720     decfactor = -decfactor;
2721     invdec = 1;
2722 else
2723     invdec = 0;
2724 end;
2725 if decfactor > sqrt(ntrials) % if large, output this many trials
2726     n = 1:ntrials;
2727     if exist('phargs') & length(phargs)>1
2728         if phargs(2)>0
2729             n = n(ceil(phargs(2)*ntrials)+1:end); % trials after rejection
2730         elseif phargs(2)<0
2731             n = n(1:floor(phargs(2)*length(n)));  % trials after rejection
2732         end
2733     elseif exist('ampargs') & length(ampargs)>1
2734         if ampargs(2)>0
2735             n = n(ceil(ampargs(2)*ntrials)+1:end); % trials after rejection
2736         elseif ampargs(2)<0
2737             n = n(1:floor(ampargs(2)*length(n)));  % trials after rejection
2738         end
2739     end
2740     if invdec
2741         decfactor = (length(n)-avewidth)/decfactor;
2742     else
2743         decfactor = length(n)/decfactor;
2744     end;
2745 end
2746 
2747 %
2748 %% %%%%%%%%%%%%%%%% Smooth data using moving average %%%%%%%%%%%%%%%%%%%%%%%%%%%
2749 %
2750 urdata = data; % save data to compute amp, coher on unsmoothed data
2751 
2752 if ~Allampsflag & ~exist('data2') % if imaging potential,
2753     if length(timeStretchRef) > 0 & length(timeStretchMarks) > 0
2754         %
2755         % Perform time-stretching here -JH %%%%%%%%%%%%%%%%
2756         %
2757         for t=1:size(data,2)
2758             M = timewarp(timeStretchMarks(t,:)', timeStretchRef');
2759             data(:,t) = M*data(:,t);
2760         end
2761     end
2762     tsurdata = data;
2763     % Time-stretching ends here %%%%%%%%%%%%%%
2764 
2765     if avewidth > 1 || decfactor > 1
2766         if Nosort == YES
2767             fprintf('Smoothing the data using a window width of %g epochs ',avewidth);
2768         else
2769             fprintf('Smoothing the sorted epochs with a %g-epoch moving window.',...
2770                 avewidth);
2771         end
2772         fprintf('\n');
2773         fprintf('  and a decimation factor of %g\n',decfactor);
2774 
2775         if ~exist('phargs') % if not phase-sorted trials
2776             [data,outtrials] = movav(data,1:ntrials,avewidth,decfactor,[],[],wt_wind);
2777             % Note: movav() here sorts using square window
2778             [outsort,outtrials] = movav(sortvar,1:ntrials,avewidth,decfactor,[],[],wt_wind);
2779 
2780         else % if phase-sorted trials, use circular / wrap-around smoothing
2781             backhalf  = floor(avewidth/2);
2782             fronthalf = floor((avewidth-1)/2);
2783             if avewidth > 2
2784                 [data,outtrials] = movav([data(:,[(end-backhalf+1):end]),...
2785                     data,...
2786                     data(:,[1:fronthalf])],...
2787                     [1:(ntrials+backhalf+fronthalf)],avewidth,decfactor,[],[],wt_wind);
2788                 [outsort,outtrials] = movav([sortvar((end-backhalf+1):end),...
2789                     sortvar,...
2790                     sortvar(1:fronthalf)],...
2791                     1:(ntrials+backhalf+fronthalf),avewidth,decfactor,[],[],wt_wind);
2792                 % Shift elements of outtrials so the first element is 1
2793                 outtrials = outtrials - outtrials(1) + 1;
2794             else % avewidth==2
2795                 [data,outtrials] = movav([data(:,end),data],...
2796                     [1:(ntrials+1)],avewidth,decfactor,[],[],wt_wind);
2797                 % Note: movav() here sorts using square window
2798                 [outsort,outtrials] = movav([sortvar(end) sortvar],...
2799                     1:(ntrials+1),avewidth,decfactor,[],[],wt_wind);
2800                 % Shift elements of outtrials so the first element is 1
2801                 outtrials = outtrials - outtrials(1) + 1;
2802             end
2803         end
2804         for index=1:length(percentiles)
2805             outpercent{index} = compute_percentile( sortvar, percentiles(index), outtrials, avewidth);
2806         end;
2807         if ~isempty(auxvar)
2808             if ~exist('phargs') % if not phase-sorted trials
2809                 [auxvar,tmp] = movav(auxvar,1:ntrials,avewidth,decfactor,[],[],wt_wind);
2810             else % if phase-sorted trials
2811                 if avewidth>2
2812                     [auxvar,tmp] = movav([auxvar(:,[(end-backhalf+1):end]),...
2813                         auxvar,...
2814                         auxvar(:,[1:fronthalf])],...
2815                         [1:(ntrials+backhalf+fronthalf)],avewidth,decfactor,[],[],wt_wind);
2816                     % Shift elements of tmp so the first element is 1
2817                     tmp = tmp - tmp(1) + 1;
2818                 else % avewidth==2
2819                     [auxvar,tmp] = movav([auxvar(:,end),auxvar],[1:(ntrials+1)],avewidth,decfactor,[],[],wt_wind);
2820                     % Shift elements of tmp so the first element is 1
2821                     tmp = tmp - tmp(1) + 1;
2822                 end
2823             end
2824         end
2825       %  if ~isempty(sortvar_limits),
2826       %      fprintf('Output data will be %d frames by %d smoothed trials.\n',...
2827       %          frames,length(outtrials));
2828       %      fprintf('Outtrials: %3.2f to %4.2f\n',min(outtrials),max(outtrials));
2829       %  end
2830     else % don't smooth
2831         outtrials = 1:ntrials;
2832         outsort = sortvar;
2833     end
2834 
2835     %
2836     %%%%%%%%%%%%%%%%%%%%%%%%% Find color axis limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2837     %
2838     if ~isempty(Caxis)
2839         mindat = Caxis(1);
2840         maxdat = Caxis(2);
2841         fprintf('Using the specified caxis range of [%g,%g].\n', mindat, maxdat);
2842     else
2843         mindat = min(min(data));
2844         maxdat = max(max(data));
2845         maxdat =  max(abs([mindat maxdat])); % make symmetrical about 0
2846         mindat = -maxdat;
2847         if ~isempty(caxfraction)
2848             adjmax = (1-caxfraction)/2*(maxdat-mindat);
2849             mindat = mindat+adjmax;
2850             maxdat = maxdat-adjmax;
2851             fprintf(...
2852                 'The caxis range will be %g times the sym. abs. data range -> [%g,%g].\n',...
2853                 caxfraction,mindat,maxdat);
2854         else
2855             fprintf(...
2856                 'The caxis range will be the sym. abs. data range -> [%g,%g].\n',...
2857                 mindat,maxdat);
2858         end
2859     end
2860 end % if ~Allampsflag & ~exist('data2')
2861 
2862 %
2863 %% %%%%%%%%%%%%%%%%%%%%%%%% Set time limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2864 %
2865 if isnan(timelimits(1))
2866     timelimits = [min(times) max(times)];
2867 end
2868 fprintf('Data will be plotted between %g and %g ms.\n',timelimits(1),timelimits(2));
2869 
2870 %
2871 %% %%%%%%%%%%% Image the aligned/sorted/smoothed data %%%%%%%%%%%%%%%%%%%%%%%%%%
2872 %
2873 if strcmpi(noshow, 'no')
2874     if ~any(isnan(coherfreq))       % if plot three time axes
2875         image_loy = 3*PLOT_HEIGHT;
2876     elseif Erpflag == YES   % elseif if plot only one time axes
2877         image_loy = 1*PLOT_HEIGHT;
2878     else                    % else plot erp-image only
2879         image_loy = 0*PLOT_HEIGHT;
2880     end
2881     gcapos=get(gca,'Position');
2882     delete(gca)
2883     if isempty(topomap)
2884         image_top = 1;
2885     else
2886         image_top = 0.9;
2887     end
2888     ax1=axes('Position',...
2889         [gcapos(1) gcapos(2)+image_loy*gcapos(4) ...
2890         gcapos(3) (image_top-image_loy)*gcapos(4)]);
2891 end;
2892 ind = isnan(data);    % find nan's in data
2893 [i j]=find(ind==1);
2894 if ~isempty(i)
2895     data(i,j) = 0;      % plot shifted nan data as 0 (=green)
2896 end
2897 
2898 %
2899 %% %%%%%%%%%%% Determine coherence freqeuncy %%%%%%%%%%%%%%%%%%%%%%%%%%
2900 %
2901 if length(coherfreq) == 2 & coherfreq(1) ~= coherfreq(2) & freq <= 0
2902     % find max frequency in specified band - should use Matlab pwelch()?
2903     if exist('psd') == 2 % from Signal Processing Toolbox
2904         [pxx,tmpfreq] = psd(urdata(:),max(1024,pow2(ceil(log2(frames)))),srate,frames,0);
2905     else % substitute from EEGLAB
2906         [pxx,tmpfreq] = spec(urdata(:),max(1024,pow2(ceil(log2(frames)))),srate,frames,0);
2907     end;
2908     pxx = 10*log10(pxx);
2909     n = find(tmpfreq >= coherfreq(1) & tmpfreq <= coherfreq(2));
2910     % [tmpfreq(n) pxx(n)]
2911     % coherfreqs = coherfreq; % save for debugging spectrum plotting
2912 
2913     if ~length(n)
2914         coherfreq = coherfreq(1);
2915     end
2916     [dummy maxx] = max(pxx(n));
2917     coherfreq = tmpfreq(n(maxx));
2918 else
2919     coherfreq = coherfreq(1);
2920 end
2921 
2922 
2923 
2924 if ~Allampsflag & ~exist('data2') %%%%%%%% Plot ERP image %%%%%%%%%%
2925 
2926     %Stretch the data array
2927     % $$$     keyboard;
2928 
2929     if strcmpi(noshow, 'no')
2930         if TIMEX
2931             h_eim=imagesc(times,outtrials,data',[mindat,maxdat]);% plot time on x-axis
2932             set(gca,'Ydir','normal');
2933             axis([timelimits(1) timelimits(2) ...
2934                 min(outtrials) max(outtrials)]);
2935         else
2936             h_eim=imagesc(outtrials,times,data,[mindat,maxdat]); % plot trials on x-axis
2937             axis([min(outtrials) max(outtrials)...
2938                 timelimits(1) timelimits(2)]);
2939         end
2940         hold on
2941         drawnow
2942     end;
2943 
2944 elseif Allampsflag %%%%%%%%%%%%%%%% Plot allamps instead of data %%%%%%%%%%%%%%
2945 
2946     if freq > 0
2947         coherfreq = mean(freq); % use phase-sort frequency
2948     end
2949 
2950     if ~isnan(signifs) % plot received significance levels
2951         fprintf(['Computing and plotting received ERSP and ITC signif. ' ...
2952             'levels...\n']);
2953         [amps,cohers,cohsig,ampsig,allamps] = ...
2954             phasecoher(urdata,length(times),srate,coherfreq,cycles,0, ...
2955             [], [], timeStretchRef, timeStretchMarks);
2956         % need to receive cohsig and ampsig to get allamps <---
2957         ampsig = signifs([1 2]); % assume these already in dB
2958         cohsig = signifs(3);
2959 
2960     elseif alpha>0 % compute significance levels
2961         fprintf('Computing and plotting %g ERSP and ITC signif. level...\n',alpha);
2962         [amps,cohers,cohsig,ampsig,allamps] = ...
2963             phasecoher(urdata,length(times),srate,coherfreq, ...
2964             cycles, alpha, [], [], ...
2965             timeStretchRef, timeStretchMarks');
2966         % need to receive cohsig and ampsig to get allamps
2967         fprintf('Coherence significance level: %g\n',cohsig);
2968 
2969     else % no plotting of significance
2970         [amps,cohers,cohsig,ampsig,allamps] = ...
2971             phasecoher(urdata,length(times),srate,coherfreq, ...
2972             cycles,0,[], [], timeStretchRef, timeStretchMarks);
2973         % need to receive cohsig and ampsig to get allamps
2974     end
2975 
2976     % fprintf('#1 Size of allamps = [%d %d]\n',size(allamps,1),size(allamps,2));
2977 
2978     base = find(times<=DEFAULT_BASELINE_END);
2979     if length(base)<2
2980         base = 1:floor(length(times)/4); % default first quarter-epoch
2981         fprintf('Using %g to %g ms as amplitude baseline.\n',...
2982             times(1),times(base(end)));
2983     end
2984 
2985 
2986     % fprintf('#2 Size of allamps = [%d %d]\n',size(allamps,1),size(allamps,2));
2987 
2988     fprintf('Subtracting the mean baseline log amplitude \n');
2989 
2990     %fprintf('Subtracting the mean baseline log amplitude %g\n',baseall);
2991     % allamps = allamps./baseall;
2992     % fprintf('#3 Size of allamps = [%d %d]\n',size(allamps,1),size(allamps,2));
2993 
2994     if avewidth > 1 || decfactor > 1
2995         if Nosort == YES
2996             fprintf(...
2997                 'Smoothing the amplitude epochs using a window width of %g epochs ',...
2998                 avewidth);
2999         else % sort trials
3000             fprintf(...
3001                 'Smoothing the sorted amplitude epochs with a %g-epoch moving window.',...
3002                 avewidth);
3003         end
3004         fprintf('\n');
3005         fprintf('  and a decimation factor of %g\n',decfactor);
3006 
3007         %fprintf('4 Size of allamps = [%d %d]\n',size(allamps,1),size(allamps,2));
3008 
3009         if exist('phargs') % if phase-sorted trials, use circular/wrap-around smoothing
3010             backhalf  = floor(avewidth/2);
3011             fronthalf = floor((avewidth-1)/2);
3012             if avewidth > 2
3013                 [allamps,outtrials] = movav([allamps(:,[(end-backhalf+1):end]),...
3014                     allamps,...
3015                     allamps(:,[1:fronthalf])],...
3016                     [1:(ntrials+backhalf+fronthalf)],avewidth,decfactor,[],[],wt_wind);
3017                 % Note: sort using square window
3018                 [outsort,outtrials] = movav([sortvar((end-backhalf+1):end),...
3019                     sortvar,...
3020                     sortvar(1:fronthalf)],...
3021                     1:(ntrials+backhalf+fronthalf),avewidth,decfactor,[],[],wt_wind);
3022                 % Shift elements of outtrials so the first element is 1
3023                 outtrials = outtrials - outtrials(1) + 1;
3024                 if ~isempty(auxvar)
3025                     [auxvar,tmp] = movav([auxvar(:,[(end-backhalf+1):end]),...
3026                         auxvar,...
3027                         auxvar(:,[1:fronthalf])],...
3028                         [1:(ntrials+backhalf+fronthalf)],avewidth,decfactor,[],[],wt_wind);
3029                     % Shift elements of outtrials so the first element is 1
3030                     outtrials = outtrials - outtrials(1) + 1;
3031                 end
3032             else % avewidth==2
3033                 [allamps,outtrials] = movav([allamps(:,end),allamps],...
3034                     [1:(ntrials+1)],avewidth,decfactor,[],[],wt_wind);
3035                 % Note: sort using square window
3036                 [outsort,outtrials] = movav([sortvar(end) sortvar],...
3037                     1:(ntrials+1),avewidth,decfactor,[],[],wt_wind);
3038                 % Shift elements of outtrials so the first element is 1
3039                 outtrials = outtrials - outtrials(1) + 1;
3040                 [auxvar,tmp] = movav([auxvar(:,end),auxvar],[1:(ntrials+1)],avewidth,decfactor,[],[],wt_wind);
3041                 % Shift elements of tmp so the first element is 1
3042                 tmp = tmp - tmp(1) + 1;
3043             end
3044         else % if trials not phase sorted, no wrap-around
3045             [allamps,outtrials] = movav(allamps,1:ntrials,avewidth,decfactor,[],[],wt_wind);
3046             %fprintf('5 Size of allamps = [%d %d]\n',size(allamps,1),size(allamps,2));
3047             [outsort,outtrials] = movav(sortvar,1:ntrials,avewidth,decfactor,[],[],wt_wind);
3048             if ~isempty(auxvar)
3049                 [auxvar,tmp] = movav(auxvar,1:ntrials,avewidth,decfactor,[],[],wt_wind);
3050             end
3051         end
3052         for index=1:length(percentiles)
3053             outpercent{index} = compute_percentile( sortvar, percentiles(index), outtrials, avewidth);
3054         end;
3055         fprintf('Output allamps data will be %d frames by %d smoothed trials.\n',...
3056             frames,length(outtrials));
3057 
3058     else % if no smoothing
3059         outtrials = 1:ntrials;
3060         outsort = sortvar;
3061     end
3062 
3063     allamps = 20*log10(allamps); % convert allamps to dB
3064     amps = 20*log10(amps);       % convert latency mean amps to dB
3065     ampsig = 20*log10(ampsig);   % convert amplitude signif thresholds to dB
3066 
3067     if alpha>0
3068         fprintf('Amplitude significance levels: [%g %g] dB\n',ampsig(1),ampsig(2));
3069     end
3070 
3071     if isnan(baseamp) % if not specified in 'limits'
3072         [amps,baseamp] = rmbase(amps,length(times),base); % subtract the dB baseline (baseamp)
3073                                                           % amps are the means at each latency
3074         allamps = allamps - baseamp; % subtract dB baseline from allamps
3075         % amplitude
3076         ampsig = ampsig - baseamp; % subtract dB baseline from ampsig
3077 
3078     else % if baseamp specified in 'limits' (as last argument 'basedB', see help)
3079         amps = amps-baseamp; % use specified (log) baseamp
3080         allamps = allamps - baseamp; % subtract dB baseline
3081         if isnan(signifs);
3082             ampsig = ampsig-baseamp; % subtract dB baseline
3083         end
3084     end
3085 
3086 
3087     %
3088     %%%%%%%%%%%%%%%%%%%%%%%%% Find color axis limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3089     %
3090     if ~isempty(Caxis)
3091         mindat = Caxis(1);
3092         maxdat = Caxis(2);
3093         fprintf('Using the specified caxis range of [%g,%g].\n',...
3094             mindat,maxdat);
3095     else
3096 
3097         % Changed -JH
3098         % 0.8 is Scott's suggestion to make the erp image show small
3099         % variations better
3100         maxdat = 0.8 * max(max(abs(allamps)));
3101         mindat = -maxdat;
3102         % $$$         mindat = min(min(allamps));
3103         % $$$         maxdat = max(max(allamps));
3104         % $$$         maxdat =  max(abs([mindat maxdat])); % make symmetrical about 0
3105         % $$$         mindat = -maxdat;
3106 
3107         if ~isempty(caxfraction)
3108             adjmax = (1-caxfraction)/2*(maxdat-mindat);
3109             mindat = mindat+adjmax;
3110             maxdat = maxdat-adjmax;
3111             fprintf(...
3112                 'The caxis range will be %g times the sym. abs. data range -> [%g,%g].\n',...
3113                 caxfraction,mindat,maxdat);
3114         else
3115             fprintf(...
3116                 'The caxis range will be the sym. abs. data range -> [%g,%g].\n',...
3117                 mindat,maxdat);
3118         end
3119     end
3120     %
3121     %%%%%%%%%%%%%%%%%%%%% Image amplitudes at coherfreq %%%%%%%%%%%%%%%%%%%%%%%%%%
3122     %
3123 
3124     if strcmpi(noshow, 'no')
3125         fprintf('Plotting amplitudes at freq %g Hz instead of potentials.\n',coherfreq);
3126         if TIMEX
3127             imagesc(times,outtrials,allamps',[mindat,maxdat]);% plot time on x-axis
3128             set(gca,'Ydir','normal');
3129             axis([timelimits(1) timelimits(2) ...
3130                 min(outtrials) max(outtrials)]);
3131         else
3132             imagesc(outtrials,times,allamps,[mindat,maxdat]); % plot trials on x-axis
3133             axis([min(outtrials) max(outtrials)...
3134                 timelimits(1) timelimits(2)]);
3135         end
3136         drawnow
3137         hold on
3138     end;
3139     data = allamps;
3140 
3141 elseif exist('data2') %%%%%% Plot allcohers instead of data %%%%%%%%%%%%%%%%%%%
3142     %%%%%%%%% UNDOCUMENTED AND DEPRECATED OPTION %%%%%%%%%%%%
3143     if freq > 0
3144         coherfreq = mean(freq); % use phase-sort frequency
3145     end
3146     if alpha>0
3147         fprintf('Computing and plotting %g coherence significance level...\n',alpha);
3148 
3149         % [amps,cohers,cohsig,ampsig,allcohers] = ...
3150         %   crosscoher(urdata,data2,length(times),srate,coherfreq,cycles,alpha);
3151 
3152         fprintf('Inter-Trial Coherence significance level: %g\n',cohsig);
3153         fprintf('Amplitude significance levels: [%g %g]\n',ampsig(1),ampsig(2));
3154     else
3155         % [amps,cohers,cohsig,ampsig,allcohers] = ...
3156         %    crosscoher(urdata,data2,length(times),srate,coherfreq,cycles,0);
3157     end
3158     if ~exist('allcohers')
3159         fprintf('erpimage(): allcohers not returned....\n')
3160         return
3161     end
3162     allamps = allcohers; % output variable
3163     % fprintf('Size allcohers = (%d, %d)\n',size(allcohers,1),size(allcohers,2));
3164     % fprintf('#1 Size of allcohers = [%d %d]\n',size(allcohers,1),size(allcohers,2));
3165     base = find(times<=0);
3166     if length(base)<2
3167         base = 1:floor(length(times)/4); % default first quarter-epoch
3168     end
3169 
3170     amps = 20*(log10(amps) - log10(mean(amps))); % convert to dB
3171     %amps = 20*log10(amps); % convert to dB
3172     ampsig = 20*(log10(ampsig) - log10(mean(amps))); % convert to dB
3173     %ampsig = 20*log10(ampsig); % convert to dB
3174 
3175     if isnan(baseamp)
3176         [amps,baseamp] = rmbase(amps,length(times),base); % remove baseline
3177     else
3178         amps = amps - baseamp;
3179     end
3180     % fprintf('#2 Size of allcohers = [%d %d]\n',size(allcohers,1),size(allcohers,2));
3181 
3182     if avewidth > 1 || decfactor > 1
3183         if Nosort == YES
3184             fprintf(...
3185                 'Smoothing the amplitude epochs using a window width of %g epochs '...
3186                 ,avewidth);
3187         else
3188             fprintf(...
3189                 'Smoothing the sorted amplitude epochs with a %g-epoch moving window.'...
3190                 ,avewidth);
3191         end
3192         fprintf('\n');
3193         fprintf('  and a decimation factor of %g\n',decfactor);
3194         % fprintf('4 Size of allcohers = [%d %d]\n',size(allcohers,1),size(allcohers,2));
3195 
3196         [allcohers,outtrials] = movav(allcohers,1:ntrials,avewidth,decfactor,[],[],wt_wind);
3197         % fprintf('5 Size of allcohers = [%d %d]\n',size(allcohers,1),size(allcohers,2));
3198         [outsort,outtrials] = movav(sortvar,1:ntrials,avewidth,decfactor,[],[],wt_wind);
3199 %        fprintf('Output data will be %d frames by %d smoothed trials.\n',...
3200 %            frames,length(outtrials));
3201     else
3202         outtrials = 1:ntrials;
3203         outsort = sortvar;
3204     end
3205     %
3206     %%%%%%%%%%%%%%%%%%%%%%%%% Find color axis limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3207     %
3208     if ~isempty(Caxis)
3209         mindat = Caxis(1);
3210         maxdat = Caxis(2);
3211         fprintf('Using the specified caxis range of [%g,%g].\n',...
3212             mindat,maxdat);
3213     else
3214         mindat = -1;
3215         maxdat = 1
3216         fprintf(...
3217             'The caxis range will be the sym. abs. data range [%g,%g].\n',...
3218             mindat,maxdat);
3219     end
3220     %
3221     %%%%%%%%%%%%%%%%%%%%% Image coherences at coherfreq %%%%%%%%%%%%%%%%%%%%%%%%%%
3222     %
3223     if strcmpi(noshow, 'no')
3224         fprintf('Plotting coherences at freq %g Hz instead of potentials.\n',coherfreq);
3225         if TIMEX
3226             imagesc(times,outtrials,allcohers',[mindat,maxdat]);% plot time on x-axis
3227             set(gca,'Ydir','normal');
3228             axis([timelimits(1) timelimits(2) ...
3229                 min(outtrials) max(outtrials)]);
3230         else
3231             imagesc(outtrials,times,allcohers,[mindat,maxdat]); % plot trials on x-axis
3232             axis([min(outtrials) max(outtrials)...
3233                 timelimits(1) timelimits(2)]);
3234         end
3235         drawnow
3236         hold on
3237     end;
3238 
3239     
3240 end 
3241 
3242 %Change limits on ERPimage y-axis if requested
3243 if ~isempty(sortvar_limits)
3244     if exist('phargs','var'),
3245         fprintf('********* Warning *********\n');
3246         fprintf('Specifying sorting variable limits has no effect when sorting by phase.\n');
3247     elseif exist('valargs','var'),
3248         fprintf('********* Warning *********\n');
3249         fprintf('Specifying sorting variable limits has no effect when sorting by mean EEG voltage.\n');
3250     elseif exist('ampargs','var'),
3251         fprintf('********* Warning *********\n');
3252         fprintf('Specifying sorting variable limits has no effect when sorting by frequency amp.\n'); 
3253     else
3254         v=axis;
3255         img_mn=find_crspnd_pt(sortvar_limits(1),outsort,outtrials);
3256         if isempty(img_mn),
3257             img_mn=1;
3258             sortvar_limits(1)=outsort(1);
3259         end
3260         img_mx=find_crspnd_pt(sortvar_limits(2),outsort,outtrials);
3261         if isempty(img_mx),
3262             img_mx=length(outsort);
3263             sortvar_limits(2)=outsort(img_mx);
3264         end
3265         axis([v(1:2) img_mn img_mx]);
3266         id1=find(sortvar>=sortvar_limits(1));
3267         id2=find(sortvar<=sortvar_limits(2));
3268         id=intersect(id1,id2);
3269         fprintf('%d epochs fall within sortvar limits.\n',length(id));
3270         urdata=urdata(:,id);
3271         if ~isempty(tsurdata),
3272             tsurdata=tsurdata(:,id);
3273         end
3274     end
3275 end
3276 %%%%%%%%%%%%%%%%%%%%%%%%%%% End plot image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3277 
3278 v=axis;
3279 fprintf('Output data will be %d frames by %d smoothed trials.\n',...
3280     frames,v(4)-v(3)+1);
3281 fprintf('Outtrials: %3.2f to %4.2f\n',v(3),v(4));
3282 
3283 %
3284 %%%%%%%%%%%%%%%%%%%%% Compute y-axis tick values and labels (if requested) %%%%%%%%%%%%%%%%%%%%%%%%%%
3285 %
3286 
3287 if ~isempty(img_ylab) && ~strcmpi(img_ylab,'Trials')
3288     %make ERPimage y-tick labels in units of sorting variable
3289     if isempty(sortvar_limits),
3290         mn=min(outsort);
3291         mx=max(outsort);
3292     else
3293         mn=sortvar_limits(1);
3294         mx=sortvar_limits(2);
3295     end
3296     ord=orderofmag(mx-mn);
3297     rng_rnd=round([mn mx]/ord)*ord;
3298     if isempty(img_ytick_lab)
3299         img_ytick_lab=[rng_rnd(1):ord:rng_rnd(2)];
3300         in_range=find((img_ytick_lab>=mn) & (img_ytick_lab<=mx));
3301         img_ytick_lab=img_ytick_lab(in_range);
3302     else
3303         img_ytick_lab=unique(img_ytick_lab); %make sure it is sorted
3304         in_range=find((img_ytick_lab>=mn) & (img_ytick_lab<=mx));
3305         if length(img_ytick_lab)~=length(in_range),
3306             fprintf('\n***Warning***\n');
3307             fprintf('''img_trialax_ticks'' exceed smoothed sorting variable values. Max/min values are %f/%f.\n\n',mn,mx);
3308             img_ytick_lab=img_ytick_lab(in_range);
3309         end
3310     end
3311     n_tick=length(img_ytick_lab);
3312     img_ytick=zeros(1,n_tick);
3313     for tickloop=1:n_tick,
3314         img_ytick(tickloop)=find_crspnd_pt(img_ytick_lab(tickloop),outsort,outtrials);
3315     end
3316 elseif ~isempty(img_ylab), %make ERPimage y-tick labels in units of Trials
3317     if isempty(img_ytick_lab)
3318         v=axis; %note: sorting variable limits have already been used to determine range of ERPimage y-axis
3319         mn=v(3);
3320         mx=v(4);
3321         ord=orderofmag(mx-mn);
3322         rng_rnd=round([mn mx]/ord)*ord;
3323         img_ytick=[rng_rnd(1):ord:rng_rnd(2)];
3324         in_range=find((img_ytick>=mn) & (img_ytick<=mx));
3325         img_ytick=img_ytick(in_range);
3326     else
3327         img_ytick=img_ytick_lab;
3328     end
3329 end
3330 
3331 
3332 %% %%%%%%%%%%%%%%%%%%%%%%%%% plot vert lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3333 if ~isempty(verttimes)
3334     if size(verttimes,1) ~= 1 & size(verttimes,2) == 1 & size(verttimes,1) ~= ntrials
3335         verttimes = verttimes';
3336     end
3337     if size(verttimes,1) ~= 1 & size(verttimes,1) ~= ntrials
3338         fprintf('\nerpimage(): vert arg matrix must have 1 or %d rows\n',ntrials);
3339         return
3340     end;
3341     if strcmpi(noshow, 'no')
3342         if size(verttimes,1) == 1
3343             fprintf('Plotting %d lines at times: ',size(verttimes,2));
3344         else
3345             fprintf('Plotting %d traces starting at times: ',size(verttimes,2));
3346         end
3347         for vt = verttimes % for each column
3348             fprintf('%g ',vt(1));
3349             if isnan(aligntime) % if NOT re-aligned data
3350                 if TIMEX          % overplot vt on image
3351                     if length(vt)==1
3352                         mydotstyle = DOTSTYLE;
3353                         if exist('auxcolors') & ...
3354                                 length(verttimes) == length(auxcolors)
3355                             mydotstyle = auxcolors{find(verttimes == vt)};
3356                         end
3357                         figure(curfig);plot([vt vt],[0 max(outtrials)],mydotstyle,'Linewidth',VERTWIDTH);
3358                     elseif length(vt)==ntrials
3359                         [outvt,ix] = movav(vt,1:ntrials,avewidth,decfactor,[],[],wt_wind);
3360                         figure(curfig);plot(outvt,outtrials,DOTSTYLE,'Linewidth',VERTWIDTH);
3361                     end
3362                 else
3363                     if length(vt)==1
3364                         figure(curfig);plot([0 max(outtrials)],[vt vt],DOTSTYLE,'Linewidth',VERTWIDTH);
3365                     elseif length(vt)==ntrials
3366                         [outvt,ix] = movav(vt,1:ntrials,avewidth,decfactor,[],[],wt_wind);
3367                         figure(curfig);plot(outtrials,outvt,DOTSTYLE,'Linewidth',VERTWIDTH);
3368                     end
3369                 end
3370             else % re-aligned data
3371                 if TIMEX          % overplot vt on image
3372                     if length(vt)==ntrials
3373                         [outvt,ix] = movav(vt,1:ntrials,avewidth,decfactor,[],[],wt_wind);
3374                         figure(curfig);plot(aligntime+outvt-outsort,outtrials,DOTSTYLE,'LineWidth',VERTWIDTH);
3375                     elseif length(vt)==1
3376                         figure(curfig);plot(aligntime+vt-outsort,outtrials,DOTSTYLE,'LineWidth',VERTWIDTH);
3377                     end
3378                 else
3379                     if length(vt)==ntrials
3380                         [outvt,ix] = movav(vt,1:ntrials,avewidth,decfactor,[],[],wt_wind);
3381                         figure(curfig);plot(outtrials,aligntime+outvt-outsort,DOTSTYLE,'LineWidth',VERTWIDTH);
3382                     elseif length(vt)==1
3383                         figure(curfig);plot(outtrials,aligntime+vt-outsort,DOTSTYLE,'LineWidth',VERTWIDTH);
3384                     end
3385                 end
3386             end
3387         end
3388         %end
3389         fprintf('\n');
3390     end;
3391 end
3392 
3393 %
3394 %% %%%%%%%%%%%%%%%%%%%%%%%%% plot horizontal ('horz') lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3395 %
3396 if ~isempty(horzepochs)
3397     if size(horzepochs,1) > 1 & size(horzepochs,1) > 1
3398         fprintf('\nerpimage(): horz arg must be a vector\n');
3399         return
3400     end;
3401     if strcmpi(noshow, 'no')
3402         if ~isempty(img_ylab) && ~strcmpi(img_ylab,'Trials'),
3403             %trial axis in units of sorting variable
3404             mx=max(outsort);
3405             mn=min(outsort);
3406             fprintf('Plotting %d lines at epochs corresponding to sorting variable values: ',length(horzepochs));
3407             for he = horzepochs % for each horizontal line
3408                 fprintf('%g ',he);
3409                 %find trial number corresponding to this value of sorting
3410                 %variable:
3411                 if (he>mn) && (he<mx)
3412                     he_ep=find_crspnd_pt(he,outsort,outtrials);
3413                     if TIMEX          % overplot he_ep on image
3414                         figure(curfig);plot([timelimits(1) timelimits(2)],[he_ep he_ep],LINESTYLE,'Linewidth',HORZWIDTH);
3415                     else
3416                         figure(curfig);plot([he_ep he_ep], [timelimits(1) timelimits(2)],LINESTYLE,'Linewidth',HORZWIDTH);
3417                     end
3418                 end
3419             end
3420             
3421         else %trial axis in units of trials
3422             fprintf('Plotting %d lines at epochs: ',length(horzepochs));
3423             for he = horzepochs % for each horizontal line
3424                 fprintf('%g ',he);
3425                 if TIMEX          % overplot he on image
3426                     figure(curfig);plot([timelimits(1) timelimits(2)],[he he],LINESTYLE,'Linewidth',HORZWIDTH);
3427                 else
3428                     figure(curfig);plot([he he], [timelimits(1) timelimits(2)],LINESTYLE,'Linewidth',HORZWIDTH);
3429                 end
3430             end
3431         end
3432         fprintf('\n');
3433     end;
3434 end
3435 if strcmpi(noshow, 'no')
3436     set(gca,'FontSize',TICKFONT)
3437     hold on;
3438 end;
3439 %
3440 %% %%%%%%%%% plot vertical line at 0 or align time %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3441 %
3442 if strcmpi(noshow, 'no')
3443     if ~isnan(aligntime) % if trials time-aligned
3444         if times(1) <= aligntime & times(frames) >= aligntime
3445             figure(curfig);plot([aligntime aligntime],[min(outtrials) max(outtrials)],...
3446                 'k','Linewidth',ZEROWIDTH); % plot vertical line at time 0
3447             % plot vertical line at aligntime
3448         end
3449     else % trials not time-aligned
3450         if times(1) <= 0 & times(frames) >= 0
3451             figure(curfig);plot([0 0],[min(outtrials) max(outtrials)],...
3452                 'k','Linewidth',ZEROWIDTH); % plot smoothed sortwvar
3453         end
3454     end
3455 end;
3456 
3457 if strcmpi(noshow, 'no') & ( min(outsort) < timelimits(1) ...
3458         |max(outsort) > timelimits(2))
3459     ur_outsort = outsort; % store the pre-adjusted values
3460     fprintf('Not all sortvar values within time vector limits: \n')
3461     fprintf('        outliers will be shown at nearest limit.\n');
3462     i = find(outsort< timelimits(1));
3463     outsort(i) = timelimits(1);
3464     i = find(outsort> timelimits(2));
3465     outsort(i) = timelimits(2);
3466 end
3467 
3468 if strcmpi(noshow, 'no')
3469     if TIMEX
3470         if Nosort == YES
3471             figure(curfig);
3472             l=ylabel(img_ylab);
3473             if ~isempty(img_ylab) && ~strcmpi(img_ylab,'Trials')
3474                 set(gca,'ytick',img_ytick,'yticklabel',img_ytick_lab);
3475             end
3476         else
3477             if exist('phargs','var')
3478                 figure(curfig);l=ylabel('Phase-sorted Trials');
3479             elseif exist('ampargs','var')
3480                 figure(curfig);l=ylabel('Amplitude-sorted Trials');
3481             elseif exist('valargs','var')
3482                 figure(curfig);l=ylabel('Voltage-sorted Trials');
3483             else
3484                 l=ylabel(img_ylab);
3485                 if ~isempty(img_ylab)
3486                       set(gca,'ytick',img_ytick);
3487                 end
3488                 if ~strcmpi(img_ylab,'Trials')
3489                     set(gca,'yticklabel',img_ytick_lab);
3490                 end
3491             end
3492         end
3493     else % if switch x<->y axes
3494         if Nosort == YES & NoTimeflag==NO
3495             figure(curfig);
3496             l=xlabel(img_ylab);
3497             if ~isempty(img_ylab) && ~strcmpi(img_ylab,'Trials')
3498                 set(gca,'xtick',img_ytick,'xticklabel',img_ytick_lab);
3499             end
3500         else
3501             if exist('phargs')
3502                 figure(curfig);l=ylabel('Phase-sorted Trials');
3503             elseif NoTimeflag == NO
3504                 figure(curfig);l=xlabel('Sorted Trials');
3505             else
3506                 l=xlabel(img_ylab);
3507                 if ~isempty(img_ylab) && ~strcmpi(img_ylab,'Trials')
3508                     set(gca,'xtick',img_ytick,'xticklabel',img_ytick_lab);
3509                 end
3510             end
3511         end
3512     end
3513     set(l,'FontSize',LABELFONT);
3514 
3515     if ~strcmpi(plotmode, 'topo')
3516         t=title(titl);
3517         set(t,'FontSize',LABELFONT);
3518     else
3519         NAME_OFFSETX = 0.1;
3520         NAME_OFFSETY = 0.2;
3521         xx = xlim; xmin = xx(1); xdiff = xx(2)-xx(1); xpos = double(xmin+NAME_OFFSETX*xdiff);
3522         yy = ylim; ymax = yy(2); ydiff = yy(2)-yy(1); ypos = double(ymax-NAME_OFFSETY*ydiff);
3523         t=text(xpos, ypos,titl); 
3524         axis off;
3525     end;
3526     
3527     set(gca,'Box','off');
3528     set(gca,'Fontsize',TICKFONT);
3529     set(gca,'color',BACKCOLOR);
3530     if Erpflag == NO & NoTimeflag == NO
3531         if exist('NoTimesPassed')~=1
3532             figure(curfig);l=xlabel('Time (ms)');
3533         else
3534             figure(curfig);l=xlabel('Frames');
3535         end
3536         set(l,'Fontsize',LABELFONT);
3537     end
3538 end;
3539 
3540 %
3541 %% %%%%%%%%%%%%%%%%%% Overplot sortvar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3542 %
3543 if strcmpi(noshow, 'no')
3544 
3545     if Noshow == YES
3546         fprintf('Not overplotting sorted sortvar on data.\n');
3547 
3548     elseif isnan(aligntime) % plot sortvar on un-aligned data
3549 
3550         if Nosort == NO;
3551             fprintf('Overplotting sorted sortvar on data.\n');
3552         end
3553         hold on;
3554         if TIMEX      % overplot sortvar
3555             figure(curfig);plot(outsort,outtrials,'k','LineWidth',SORTWIDTH);
3556         else
3557             figure(curfig);plot(outtrials,outsort,'k','LineWidth',SORTWIDTH);
3558         end
3559         drawnow
3560     else % plot re-aligned zeros on sortvar-aligned data
3561         if Nosort == NO;
3562             fprintf('Overplotting sorted sortvar on data.\n');
3563         end
3564         hold on;
3565         if TIMEX      % overplot re-aligned 0 time on image
3566             figure(curfig);plot([aligntime aligntime],[min(outtrials) max(outtrials)],...
3567                 'k','LineWidth',SORTWIDTH);
3568         else
3569             figure(curfig);plot([[min(outtrials) max(outtrials)],aligntime aligntime],...
3570                 'k','LineWidth',SORTWIDTH);
3571         end
3572         fprintf('Overplotting realigned times-zero on data.\n');
3573         hold on;
3574 
3575         if TIMEX      % overplot realigned sortvar on image
3576             figure(curfig);plot(0+aligntime-outsort,outtrials,'k','LineWidth',ZEROWIDTH);
3577         else
3578             figure(curfig);plot(0+outtrials,aligntime-outsort,'k','LineWidth',ZEROWIDTH);
3579         end
3580         drawnow
3581     end
3582 end;
3583 
3584 %
3585 %% %%%%%%%%%%%%%%%%%% Overplot auxvar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3586 %
3587 if strcmpi(noshow, 'no')
3588     if ~isempty(auxvar)
3589         fprintf('Overplotting auxvar(s) on data.\n');
3590         hold on;
3591         auxtrials = outtrials(:)' ; % make row vector
3592 
3593         if exist('auxcolors')~=1 % If no auxcolors specified
3594             auxcolors = cell(1,size(auxvar,1));
3595             for c=1:size(auxvar,1)
3596                 auxcolors(c) = {'k'};       % plot auxvars as black trace(s)
3597             end
3598         end
3599         if length(auxcolors) < size(auxvar,1)
3600             nauxColors = length(auxcolors);
3601             for k=nauxColors+1:size(auxvar,1)
3602                 auxcolors = { auxcolors{:} auxcolors{1+rem(k-1,nauxColors)}};
3603             end
3604         end
3605         for c=1:size(auxvar,1)
3606             auxcolor = auxcolors{c};
3607             if ~isempty(auxcolor)
3608               if isnan(aligntime) % plot auxvar on un-aligned data
3609                   if TIMEX      % overplot auxvar
3610                     figure(curfig);plot(auxvar(c,:)',auxtrials',auxcolor,'LineWidth',SORTWIDTH);
3611                   else
3612                     figure(curfig);plot(auxtrials',auxvar(c,:)',auxcolor,'LineWidth',SORTWIDTH);
3613                   end
3614                   drawnow
3615               else % plot re-aligned zeros on sortvar-aligned data
3616                   if TIMEX      % overplot realigned 0-time on image
3617                     figure(curfig);plot(auxvar(c,:)',auxtrials',auxcolor,'LineWidth',ZEROWIDTH);
3618                   else
3619                     figure(curfig);plot(0+auxtrials',aligntime-auxvar(c,:)',auxcolor,'LineWidth',ZEROWIDTH);
3620                   end
3621                   drawnow
3622               end % aligntime
3623            end % if auxcolor
3624         end % c
3625     end % auxvar
3626     if exist('outpercent')
3627         for index = 1:length(outpercent)
3628             if isnan(aligntime) % plot auxvar on un-aligned data
3629                 figure(curfig); plot(outpercent{index},outtrials,'k','LineWidth',SORTWIDTH);
3630             else
3631                 figure(curfig); plot(aligntime-outpercent{index},outtrials,'k','LineWidth',SORTWIDTH);
3632             end;
3633         end;
3634     end;
3635 end;
3636 %
3637 %% %%%%%%%%%%%%%%%%%%%%%% Plot colorbar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3638 %
3639 if strcmpi(noshow, 'no')
3640     if Colorbar == YES
3641         pos=get(ax1,'Position');
3642         axcb=axes('Position',...
3643             [pos(1)+pos(3)+0.02 pos(2) ...
3644             0.03 pos(4)]);
3645         figure(curfig);cbar(axcb,0,[mindat,maxdat]); % plot colorbar to right of image
3646         title(cbar_title);
3647         set(axcb,'fontsize',TICKFONT,'xtick',[]);
3648         % drawnow
3649         axes(ax1); % reset current axes to the erpimage
3650     end
3651 end;
3652 
3653 %
3654 %% %%%%%%%%%%%%%%%%%%%%% Compute ERP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3655 %
3656 erp = [];
3657 
3658 if Erpflag == YES
3659     if exist('erpalpha')
3660         if erp_ptiles>1,
3661            fprintf(['\nOnly plotting one ERP (i.e., not plotting ERPs from %d percentile splits) ' ...
3662                'because ''erpalpha'' option was chosen.  You can''t plot both.\n\n'],erp_ptiles)
3663         end
3664         [erp erpsig] = nan_mean(fastif(length(tsurdata) > 0, tsurdata',urdata'), ...
3665             erpalpha);
3666         fprintf('   Mean ERP (p<%g) significance threshold: +/-%g\n', ...
3667             erpalpha,mean(erpsig));
3668     else
3669         %potentially make ERPs of 50%, 33%, or 25% split of trials
3670         n_trials=size(urdata,2);
3671         trials_step=round(n_trials/erp_ptiles);
3672         erp=zeros(erp_ptiles,size(urdata,1));
3673         for ploop=1:erp_ptiles,
3674             ptile_trials=[1:trials_step]+(ploop-1)*trials_step;
3675             if max(ptile_trials)>n_trials,
3676                ptile_trials=ptile_trials(1):n_trials; 
3677             end
3678             erp(ploop,:)=nan_mean(fastif(length(tsurdata) > 0, ...
3679                 tsurdata(:,ptile_trials)',urdata(:,ptile_trials)'));
3680         end
3681        % else
3682             %orig line
3683             %[erp] = nan_mean(fastif(length(tsurdata) > 0, tsurdata', urdata'));
3684         %end
3685     end % compute average ERP, ignoring nan's
3686 end;
3687 
3688 if Erpflag == YES & strcmpi(noshow, 'no')
3689     axes(ax1); % reset current axes to the erpimage
3690     xtick = get(ax1,'Xtick');               % remember x-axis tick locations
3691     xticklabel = get(ax1,'Xticklabel');     % remember x-axis tick locations
3692     set(ax1, 'xticklabel', []);
3693     widthxticklabel = size(xticklabel,2);
3694     xticklabel = cellstr(xticklabel);
3695     for tmpindex = 1:length(xticklabel)
3696         if length(xticklabel{tmpindex}) < widthxticklabel
3697             spaces = char(ones(1,ceil((widthxticklabel-length(xticklabel{tmpindex}))/2) )*32);
3698             xticklabel{tmpindex} = [spaces xticklabel{tmpindex}];
3699         end;
3700     end;
3701     xticklabel = strvcat(xticklabel);
3702     if Erpstdflag == YES
3703         stdev = nan_std(urdata');
3704     end;
3705     %
3706     %%%%%% Plot ERP time series below image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3707     %
3708     if isnan(maxerp)
3709         fac = 10;
3710         maxerp = 0;
3711         while maxerp == 0
3712             maxerp = round(fac*YEXPAND*max(erp))/fac; % minimal decimal places
3713             fac = 10*fac;
3714         end
3715         if Erpstdflag == YES
3716             fac = fac/10;
3717             maxerp = max(maxerp, round(fac*YEXPAND*max(erp+stdev))/fac);
3718         end;
3719         if ~isempty(erpsig)
3720             erpsig = [erpsig;-1*erpsig];
3721             maxerp = max(maxerp, round(fac*YEXPAND*max(erpsig))/fac);
3722         end
3723         maxerp=max(maxerp);
3724     end
3725     if isnan(minerp)
3726         fac = 1;
3727         minerp = 0;
3728         while minerp == 0
3729             minerp = round(fac*YEXPAND*min(erp))/fac; % minimal decimal places
3730             fac = 10*fac;
3731         end
3732         if Erpstdflag == YES
3733             fac = fac/10;
3734             minerp = min(minerp, round(fac*YEXPAND*min(erp-stdev))/fac);
3735         end;
3736         if ~isempty(erpsig)
3737             minerp = min(minerp, round(fac*YEXPAND*min(erpsig))/fac);
3738         end
3739         minerp=min(minerp);
3740     end
3741     limit = [timelimits(1:2) minerp maxerp];
3742     
3743     if ~isnan(coherfreq)
3744         set(ax1,'Xticklabel',[]);     % remove tick labels from bottom of image
3745         ax2=axes('Position',...
3746             [gcapos(1) gcapos(2)+2/3*image_loy*gcapos(4) ...
3747             gcapos(3) (image_loy/3-YGAP)*gcapos(4)]);
3748     else
3749         ax2=axes('Position',...
3750             [gcapos(1) gcapos(2) ...
3751             gcapos(3) image_loy*gcapos(4)]);
3752     end
3753     fprintf('Plotting the ERP trace below the ERP image\n');
3754     if Erpstdflag == YES
3755         if Showwin
3756             tmph = plot1trace(ax2,times,erp,limit,[],stdev,times(winlocs),erp_grid,erp_vltg_ticks); % plot ERP +/-stdev
3757         else
3758             tmph = plot1trace(ax2,times,erp,limit, [], stdev,[],erp_grid,erp_vltg_ticks); % plot ERP +/-stdev
3759         end
3760     elseif ~isempty('erpsig')
3761         if Showwin
3762             tmph = plot1trace(ax2,times,erp,limit,erpsig,[],times(winlocs),erp_grid,erp_vltg_ticks); % plot ERP and 0+/-alpha threshold
3763         else
3764             tmph = plot1trace(ax2,times,erp,limit,erpsig,[],[],erp_grid,erp_vltg_ticks); % plot ERP and 0+/-alpha threshold
3765         end
3766     else % plot ERP alone - no significance or std dev plotted
3767         if Showwin
3768             tmph = plot1trace(ax2,times,erp,limit,[],[],times(winlocs),erp_grid,erp_vltg_ticks); % plot ERP alone
3769         else
3770             tmph = plot1trace(ax2,times,erp,limit,[],[],[],erp_grid,erp_vltg_ticks); % plot ERP alone
3771         end
3772     end;
3773 
3774     if ~isnan(aligntime)
3775         line([aligntime aligntime],[limit(3:4)*1.1],'Color','k','LineWidth',ZEROWIDTH); % x=median sort value
3776         line([0 0],[limit(3:4)*1.1],'Color','k','LineWidth',ZEROWIDTH); % x=median sort value
3777         % remove y axis
3778         if length(tmph) > 1
3779             delete(tmph(end));
3780         end;
3781     end
3782 
3783     set(ax2,'Xtick',xtick);        % use same Xticks as erpimage above
3784     if ~isnan(coherfreq)
3785         set(ax2,'Xticklabel',[]);    % remove tick labels from ERP x-axis
3786     else % bottom axis
3787         set(ax2,'Xticklabel',xticklabel); % add ticklabels to ERP x-axis
3788     end
3789 
3790     if isnan(coherfreq)            % if no amp and coher plots below . . .
3791         if TIMEX & NoTimeflag == NO
3792             if exist('NoTimesPassed')~=1
3793                 figure(curfig);l=xlabel('Time (ms)');
3794             else
3795                 figure(curfig);l=xlabel('Frames');
3796             end
3797             set(l,'FontSize',LABELFONT);
3798             % $$$       else
3799             % $$$         if exist('NoTimesPassed')~=1
3800             % $$$           figure(curfig);l=ylabel('Time (ms)');
3801             % $$$         else
3802             % $$$           figure(curfig);l=ylabel('Frames');
3803             % $$$         end
3804             % $$$         set(l,'FontSize',LABELFONT);
3805         end
3806     end
3807 
3808     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot vert lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3809 
3810     if ~isempty(verttimes)
3811         if size(verttimes,1) == ntrials
3812             vts=sort(verttimes);
3813             vts = vts(ceil(ntrials/2),:); % plot median verttimes values if a matrix
3814         else
3815             vts = verttimes(:)';  % make verttimes a row vector
3816         end
3817         for vt = vts
3818             if isnan(aligntime)
3819                 if TIMEX      % overplot vt on ERP
3820                     mydotstyle = DOTSTYLE;
3821                     if exist('auxcolors') & ...
3822                             length(verttimes) == length(verttimesColors)
3823                         mydotstyle = verttimesColors{find(verttimes == vt)};
3824                     end
3825                     figure(curfig);plot([vt vt],[limit(3:4)],mydotstyle,'Linewidth',VERTWIDTH);
3826                 else
3827                     figure(curfig);plot([min(outtrials) max(outtrials)],[limit(3:4)],DOTSTYLE,...
3828                         'Linewidth',VERTWIDTH);
3829                 end
3830             else
3831                 if TIMEX      % overplot realigned vt on ERP
3832                     figure(curfig);plot(repmat(median(aligntime+vt-outsort),1,2),[limit(3),limit(4)],...
3833                         DOTSTYLE,'LineWidth',VERTWIDTH);
3834                 else
3835                     figure(curfig);plot([limit(3),limit(4)],repmat(median(aligntime+vt-outsort),1,2),...
3836                         DOTSTYLE,'LineWidth',VERTWIDTH);
3837                 end
3838             end
3839         end
3840     end
3841 
3842     limit = double(limit);
3843     ydelta = double(1/10*(limit(2)-limit(1)));
3844     ytextoffset = double(limit(1)-1.1*ydelta);
3845     ynumoffset  = double(limit(1)-0.3*ydelta); % double for Matlab 7
3846 
3847     %Far left axis max and min labels not needed now that there are tick
3848     %marks
3849     %t=text(ynumoffset,0.7*limit(3), num2str(limit(3)));
3850     %set(t,'HorizontalAlignment','right','FontSize',TICKFONT)
3851     %t=text(ynumoffset,0.7*limit(4), num2str(limit(4)));
3852     %set(t,'HorizontalAlignment','right','FontSize',TICKFONT)
3853 
3854     ynum = 0.7*(limit(3)+limit(4))/2;
3855     t=text(ytextoffset,ynum,yerplabel,'Rotation',90);
3856     set(t,'HorizontalAlignment','center','FontSize',LABELFONT)
3857 
3858     if ~exist('YDIR')
3859         error('Default YDIR not read from ''icadefs.m''');
3860     end
3861     if YDIR == 1
3862         set(ax2,'ydir','normal')
3863     else
3864         set(ax2,'ydir','reverse')
3865     end
3866 
3867     set(ax2,'Fontsize',TICKFONT);
3868     set(ax2,'Box','off','color',BACKCOLOR);
3869     drawnow
3870 end
3871 
3872 %
3873 %% %%%%%%%%%%%%%%%%%%% Plot amp, coher time series %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3874 %
3875 if ~isnan(coherfreq)
3876     if freq > 0
3877         coherfreq = mean(freq); % use phase-sort frequency
3878     end
3879     %
3880     %%%%%% Plot amp axis below ERP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3881     %
3882     if ~Allampsflag %%%% don't repeat computation if already done for 'allamps'
3883 
3884         fprintf('Computing and plotting amplitude at %g Hz.\n',coherfreq);
3885 
3886         if ~isnan(signifs) | Cohsigflag==NO % don't compute or plot signif. levels
3887             [amps,cohers] = phasecoher(urdata,size(times,2),srate,coherfreq,cycles);
3888             if ~isnan(signifs)
3889                 ampsig = signifs([1 2]);
3890                 fprintf('Using supplied amplitude significance levels: [%g,%g]\n',...
3891                     ampsig(1),ampsig(2));
3892                 cohsig = signifs(3);
3893                 fprintf('Using supplied coherence significance level: %g\n',cohsig);
3894             end
3895         else % compute amps, cohers with significance
3896             fprintf(...
3897                 'Computing and plotting %g coherence significance level at %g Hz...\n',...
3898                 alpha,                          coherfreq);
3899             [amps,cohers,cohsig,ampsig] = ...
3900                 phasecoher(urdata,size(times,2),srate,coherfreq,cycles,alpha);
3901             fprintf('Coherence significance level: %g\n',cohsig);
3902             ampsig = 20*log10(ampsig); % convert to dB
3903         end
3904         amps = 20*log10(amps); % convert to dB
3905 
3906         if isnan(baseamp) % if baseamp not specified in 'limits'
3907             base = find(times<=DEFAULT_BASELINE_END); % use default baseline end point (ms)
3908             if length(base)<2
3909                 base = 1:floor(length(times)/4); % default first quarter-epoch
3910                 fprintf('Using %g to %g ms as amplitude baseline.\n',...
3911                     times(1),times(base(end)));
3912             end
3913             [amps,baseamp] = rmbase(amps,length(times),base); % remove dB baseline
3914             fprintf('Removed baseline amplitude of %d dB for plotting.\n',baseamp);
3915 
3916         else % if 'basedB' specified in 'limits' (in dB)
3917             fprintf('Removing specified baseline amplitude of %d dB for plotting.\n',...
3918                 baseamp);
3919             amps = amps-baseamp; % remove specified dB baseline
3920         end
3921 
3922         fprintf('Data amplitude levels: [%g %g] dB\n',min(amps),max(amps));
3923 
3924         if alpha>0 % if computed significance levels
3925             ampsig = ampsig - baseamp;
3926             fprintf('Data amplitude significance levels: [%g %g] dB\n',ampsig(1),ampsig(2));
3927         end
3928 
3929     end % ~Allampsflag
3930 
3931     if strcmpi(noshow, 'no')
3932         axis('off') % rm ERP axes axis and labels
3933         
3934         v=axis;
3935         minampERP=v(3);
3936         maxampERP=v(4);
3937         
3938         if ~exist('ynumoffset')
3939             limit = [timelimits(1:2) -max(abs([minampERP maxampERP])) max(abs([minampERP maxampERP]))];
3940             limit = double(limit);
3941             ydelta = double(1/10*(limit(2)-limit(1)));
3942             ytextoffset = double(limit(1)-1.1*ydelta);
3943             ynumoffset  = double(limit(1)-0.3*ydelta); % double for Matlab 7
3944         end
3945 
3946         t=text(ynumoffset,maxampERP,num2str(maxampERP,3));
3947         set(t,'HorizontalAlignment','right','FontSize',TICKFONT);
3948 
3949         t=text(ynumoffset,minampERP, num2str(minampERP,3));
3950         set(t,'HorizontalAlignment','right','FontSize',TICKFONT);    
3951         
3952         ax3=axes('Position',...
3953             [gcapos(1) gcapos(2)+1/3*image_loy*gcapos(4) ...
3954             gcapos(3) (image_loy/3-YGAP)*gcapos(4)]);
3955 
3956         if isnan(maxamp) % if not specified
3957             fac = 1;
3958             maxamp = 0;
3959             while maxamp == 0
3960                 maxamp = floor(YEXPAND*fac*max(amps))/fac; % minimal decimal place
3961                 fac = 10*fac;
3962             end
3963             maxamp = maxamp + 10/fac;
3964             if Cohsigflag
3965                 if ampsig(2)>maxamp
3966                     if ampsig(2)>0
3967                         maxamp = 1.01*(ampsig(2));
3968                     else
3969                         maxamp = 0.99*(ampsig(2));
3970                     end
3971                 end
3972             end
3973         end
3974         if isnan(maxamp), maxamp = 0; end   % In case the above iteration went on
3975                                             % until fac = Inf and maxamp = NaN again.
3976 
3977         if isnan(minamp) % if not specified
3978             fac = 1;
3979             minamp = 0;
3980             while minamp == 0
3981                 minamp = floor(YEXPAND*fac*max(-amps))/fac; % minimal decimal place
3982                 fac = 10*fac;
3983             end
3984             minamp = minamp + 10/fac;
3985             minamp = -minamp;
3986             if Cohsigflag
3987                 if ampsig(1)< minamp
3988                     if ampsig(1)<0
3989                         minamp = 1.01*(ampsig(1));
3990                     else
3991                         minamp = 0.99*(ampsig(1));
3992                     end
3993                 end
3994             end
3995         end
3996         if isnan(minamp), minamp = 0; end   % In case the above iteration went on
3997                                             % until fac = Inf and minamp = NaN again.
3998 
3999         fprintf('Plotting the ERSP amplitude trace below the ERP\n');
4000         fprintf('Min, max plotting amplitudes: [%g, %g] dB\n',minamp,maxamp);
4001         fprintf('     relative to baseamp: %g dB\n',baseamp);
4002         if Cohsigflag
4003             ampsiglims = [repmat(ampsig(1)-mean(ampsig),1,length(times))];
4004             ampsiglims = [ampsiglims;-1*ampsiglims];
4005             plot1trace(ax3,times,amps,[timelimits minamp(1) maxamp(1)],ampsiglims,[],[],0); % plot AMP
4006         else
4007             plot1trace(ax3,times,amps,[timelimits minamp(1) maxamp(1)],[],[],[],0); % plot AMP
4008         end
4009 
4010         if ~isnan(aligntime)
4011             line([aligntime aligntime],[minamp(1) maxamp(1)]*1.1,'Color','k');
4012             % x=median sort value
4013         end
4014         if exist('xtick') % Added -JH
4015             set(ax3,'Xtick',xtick);
4016         end
4017         set(ax3,'Xticklabel',[]);   % remove tick labels from bottom of image
4018         set(ax3,'Yticklabel',[]);   % remove tick labels from left of image
4019         set(ax3,'YColor',BACKCOLOR);
4020         axis('off');
4021         set(ax3,'Box','off','color',BACKCOLOR);
4022 
4023         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot vert marks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4024 
4025         if ~isempty(verttimes)
4026             if size(verttimes,1) == ntrials
4027                 vts=sort(verttimes);
4028                 vts = vts(ceil(ntrials/2),:); % plot median values if a matrix
4029             else
4030                 vts=verttimes(:)';
4031             end
4032             for vt = vts
4033                 if isnan(aligntime)
4034                     if TIMEX      % overplot vt on amp
4035                         mydotstyle = DOTSTYLE;
4036                         if exist('auxcolors') & ...
4037                                 length(verttimes) == length(verttimesColors)
4038                             mydotstyle = verttimesColors{find(verttimes == vt)};
4039                         end
4040 
4041                         figure(curfig);plot([vt vt],[minamp(1) maxamp(1)],mydotstyle,...
4042                             'Linewidth',VERTWIDTH);
4043                     else
4044                         figure(curfig);plot([min(outtrials) max(outtrials)],[minamp(1) maxamp(1)], ...
4045                             DOTSTYLE,...
4046                             'Linewidth',VERTWIDTH);
4047                     end
4048                 else
4049                     if TIMEX      % overplot realigned vt on amp
4050                         figure(curfig);plot(repmat(median(aligntime+vt-outsort),1,2), ...
4051                             [minamp(1),maxamp(1)],DOTSTYLE,...
4052                             'LineWidth',VERTWIDTH);
4053                     else
4054                         figure(curfig);plot([minamp,maxamp],repmat(median(aligntime+vt-outsort),1,2), ...
4055                             DOTSTYLE,...
4056                             'LineWidth',VERTWIDTH);
4057                     end
4058                 end
4059             end
4060         end
4061 
4062         if 0 % Cohsigflag % plot amplitude significance levels
4063             hold on
4064             figure(curfig);plot([timelimits(1) timelimits(2)],[ampsig(1) ampsig(1)] - mean(ampsig),'r',...
4065                 'linewidth',SIGNIFWIDTH);
4066             figure(curfig);plot([timelimits(1) timelimits(2)],[ampsig(2) ampsig(2)] - mean(ampsig),'r',...
4067                 'linewidth',SIGNIFWIDTH);
4068         end
4069 
4070         if ~exist('ynumoffset')
4071             limit = [timelimits(1:2) -max(abs([minamp maxamp])) max(abs([minamp maxamp]))];
4072             limit = double(limit);
4073             ydelta = double(1/10*(limit(2)-limit(1)));
4074             ytextoffset = double(limit(1)-1.1*ydelta);
4075             ynumoffset  = double(limit(1)-0.3*ydelta); % double for Matlab 7
4076         end
4077 
4078         t=text(ynumoffset,maxamp, num2str(maxamp,3));
4079         set(t,'HorizontalAlignment','right','FontSize',TICKFONT);
4080 
4081         t=text(ynumoffset,0, num2str(0));
4082         set(t,'HorizontalAlignment','right','FontSize',TICKFONT);
4083 
4084         t=text(ynumoffset,minamp, num2str(minamp,3));
4085         set(t,'HorizontalAlignment','right','FontSize',TICKFONT);
4086 
4087         t=text(ytextoffset,(maxamp+minamp)/2,'ERSP','Rotation',90);
4088         set(t,'HorizontalAlignment','center','FontSize',LABELFONT);
4089 
4090         axtmp = axis;
4091         dbtxt= text(1/13*(axtmp(2)-axtmp(1))+axtmp(1), ...
4092             11/13*(axtmp(4)-axtmp(3))+axtmp(3), ...
4093             [num2str(baseamp,4) ' dB']);
4094         set(dbtxt,'fontsize',TICKFONT);
4095         drawnow;
4096         set(ax3, 'xlim', timelimits);
4097         set(ax3, 'ylim', [minamp(1) maxamp(1)]);
4098 
4099         %
4100         %%%%%% Make coher axis below amp %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4101         %
4102         ax4=axes('Position',...
4103             [gcapos(1) gcapos(2) ...
4104             gcapos(3) (image_loy/3-YGAP)*gcapos(4)]);
4105         if isnan(maxcoh)
4106             fac = 1;
4107             maxcoh = 0;
4108             while maxcoh == 0
4109                 maxcoh = floor(YEXPAND*fac*max(cohers))/fac; % minimal decimal place
4110                 fac = 10*fac;
4111             end
4112             maxcoh = maxcoh + 10/fac;
4113             if maxcoh>1
4114                 maxcoh=1; % absolute limit
4115             end
4116         end
4117         if isnan(mincoh)
4118             mincoh = 0;
4119         end
4120         fprintf('Plotting the ITC trace below the ERSP\n');
4121         if Cohsigflag % plot coherence significance level
4122             cohsiglims = [repmat(cohsig,1,length(times));zeros(1,length(times))];
4123             coh_handle = plot1trace(ax4,times,cohers,[timelimits mincoh maxcoh],cohsiglims,[],[],0);
4124             % plot COHER, fill sorting window
4125         else
4126             coh_handle = plot1trace(ax4,times,cohers,[timelimits mincoh maxcoh],[],[],[],0); % plot COHER
4127         end
4128         if ~isnan(aligntime)
4129             line([aligntime aligntime],[[mincoh maxcoh]*1.1],'Color','k');
4130             % x=median sort value
4131         end
4132         % set(ax4,'Xticklabel',[]);    % remove tick labels from bottom of
4133         % image
4134         if exist('xtick')
4135             set(ax4,'Xtick',xtick);
4136             set(ax4,'Xticklabel',xticklabel);
4137         end
4138         set(ax4,'Ytick',[]);
4139         set(ax4,'Yticklabel',[]);      % remove tick labels from left of image
4140         set(ax4,'YColor',BACKCOLOR);
4141 
4142         if ~isempty(verttimes)
4143           if size(verttimes,1) == ntrials
4144                 vts=sort(verttimes);
4145                 vts = vts(ceil(ntrials/2),:); % plot median values if a matrix
4146           else
4147                 vts = verttimes(:)';  % make verttimes a row vector
4148           end
4149           for vt = vts
4150               if isnan(aligntime)
4151                if TIMEX      % overplot vt on coher
4152                    mydotstyle = DOTSTYLE;
4153                    if exist('auxcolors') & ...
4154                       length(verttimes) == length(verttimesColors)
4155                       mydotstyle = verttimesColors{find(verttimes == vt)};
4156                    end
4157 
4158                    figure(curfig);plot([vt vt],[mincoh maxcoh],mydotstyle,'Linewidth',VERTWIDTH);
4159                else
4160                    figure(curfig);plot([min(outtrials) max(outtrials)],...
4161                                                    [mincoh maxcoh],DOTSTYLE,'Linewidth',VERTWIDTH);
4162                end
4163              else
4164                if TIMEX      % overplot realigned vt on coher
4165                 figure(curfig);plot(repmat(median(aligntime+vt-outsort),1,2),...
4166                                                    [mincoh,maxcoh],DOTSTYLE,'LineWidth',VERTWIDTH);
4167                else
4168                 figure(curfig);plot([mincoh,maxcoh],repmat(median(aligntime+vt-outsort),1,2),...
4169                                                    DOTSTYLE,'LineWidth',VERTWIDTH);
4170                end
4171              end
4172           end
4173         end
4174 
4175         t=text(ynumoffset,0, num2str(0));
4176         set(t,'HorizontalAlignment','right','FontSize',TICKFONT);
4177 
4178         t=text(ynumoffset,maxcoh, num2str(maxcoh));
4179         set(t,'HorizontalAlignment','right','FontSize',TICKFONT);
4180 
4181         t=text(ytextoffset,maxcoh/2,'ITC','Rotation',90);
4182         set(t,'HorizontalAlignment','center','FontSize',LABELFONT);
4183         drawnow
4184 
4185         %if Cohsigflag % plot coherence significance level
4186         %hold on
4187         %plot([timelimits(1) timelimits(2)],[cohsig cohsig],'r',...
4188         %'linewidth',SIGNIFWIDTH);
4189         %end
4190 
4191         set(ax4,'Box','off','color',BACKCOLOR);
4192         set(ax4,'Fontsize',TICKFONT);
4193         if NoTimeflag==NO
4194             if exist('NoTimesPassed')~=1
4195                 figure(curfig);l=xlabel('Time (ms)');
4196             else
4197                 figure(curfig);l=xlabel('Frames');
4198             end
4199             set(l,'Fontsize',LABELFONT);
4200         end
4201         axtmp = axis;
4202         hztxt=text(10/13*(axtmp(2)-axtmp(1))+axtmp(1), ...
4203             8/13*(axtmp(4)-axtmp(3))+axtmp(3), ...
4204             [num2str(coherfreq,4) ' Hz']);
4205         set(hztxt,'fontsize',TICKFONT);
4206     end;% noshow
4207 else
4208     amps   = [];    % null outputs unless coherfreq specified
4209     cohers = [];
4210 end
4211 axhndls = [ax1 axcb ax2 ax3 ax4];
4212 if exist('ur_outsort')
4213     outsort = ur_outsort; % restore outsort clipped values, if any
4214 end
4215 if nargout<1
4216     data = []; % don't spew out data if no args out and no ;
4217 end
4218 
4219 %
4220 %% %%%%%%%%%%%%% Plot a topoplot() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4221 %
4222 if (~isempty(topomap)) & strcmpi(noshow, 'no')
4223     h(12)=axes('Position',...
4224         [gcapos(1)+0.10*gcapos(3) gcapos(2)+0.92*gcapos(4),...
4225         0.20*gcapos(3) 0.14*gcapos(4)]);
4226     % h(12) = subplot('Position',[.10 .86 .20 .14]);
4227     fprintf('Plotting a topo map in upper left.\n');
4228     eloc_info.plotrad = [];
4229     if length(topomap) == 1
4230         try
4231             topoplot(topomap,eloc_file,'electrodes','off', ...
4232                 'style', 'blank', 'emarkersize1chan', 10, 'chaninfo', eloc_info);
4233         catch
4234             fprintf('topoplot() plotting failed.\n');
4235         end;
4236     else
4237         try
4238             topoplot(topomap,eloc_file,'electrodes','off', 'chaninfo', eloc_info);
4239         catch
4240             fprintf('topoplot() plotting failed.\n');
4241         end;
4242     end;
4243     axis('square')
4244     axhndls = [axhndls h(12)];
4245 end
4246 
4247 %
4248 %% %%%%%%%%%%%%% Plot a spectrum %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4249 %
4250 SPECFONT = 10;
4251 if (~isempty(lospecHz)) & strcmpi(noshow, 'no')
4252     h(13)=axes('Position',...
4253         [gcapos(1)+0.82*gcapos(3) ...
4254         gcapos(2)+0.96*gcapos(4),...
4255         0.15*gcapos(3)*(0.8/gcapos(3))^0.5 ...
4256         0.10*gcapos(4)*(0.8/gcapos(4))^0.5]);
4257 
4258     % h(13) = subplot('Position',[.75 .88 .15 .10]);
4259     fprintf('Plotting the data spectrum in upper right.\n');
4260     winlength = frames;
4261     if winlength > 512
4262         for k=2:5
4263             if rem(winlength,k) == 0
4264                 break
4265             end
4266         end
4267         winlength = winlength/k;
4268     end
4269 
4270     % [Pxx, Pxxc, F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP,P)
4271     if exist('psd') == 2
4272         [Pxx,F] = psd(reshape(urdata,1,size(urdata,1)*size(urdata,2)),...
4273             max(1024,pow2(ceil(log2(frames)))),srate,frames,0,0.05);
4274         % [Pxx,F] = psd(reshape(urdata,1,size(urdata,1)*size(urdata,2)),512,srate,winlength,0,0.05);
4275     else
4276         [Pxx,F] = spec(reshape(urdata,1,size(urdata,1)*size(urdata,2)),...
4277             max(1024,pow2(ceil(log2(frames)))),srate,frames,0);
4278         % [Pxx,F] = spec(reshape(urdata,1,size(urdata,1)*size(urdata,2)),512,srate,winlength,0);
4279     end;
4280 
4281     figure(curfig);plot(F,10*log10(Pxx));
4282     goodfs = find(F>= lospecHz & F <= hispecHz);
4283     maxgfs = max(10*log10(Pxx(goodfs)));
4284     mingfs = min(10*log10(Pxx(goodfs)));
4285     axis('square')
4286     axis([lospecHz hispecHz mingfs-1 maxgfs+1]);
4287     set(h(13),'Box','off','color',BACKCOLOR);
4288     set(h(13),'Fontsize',SPECFONT);
4289     figure(curfig);l=ylabel('dB');
4290     set(l,'Fontsize',SPECFONT);
4291     if ~isnan(coherfreq)
4292         hold on; figure(curfig);plot([coherfreq,coherfreq],[mingfs maxgfs],'r');
4293     end
4294     axhndls = [axhndls h(13)];
4295 end
4296 
4297 %
4298 %% %%%%%%%%%%%%% save plotting limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4299 %
4300 limits = [min(times) max(times) minerp maxerp minamp maxamp mincoh maxcoh];
4301 limits = [limits baseamp coherfreq];  % add coherfreq to output limits array
4302 
4303 %
4304 %% %%%%%%%%%%%%% turn on axcopy() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4305 %
4306 if strcmpi(noshow, 'no')
4307     axcopy(gcf);
4308     % 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';
4309     % axcopy(gcf,eegstr);
4310 end;
4311 
4312 
4313 % returning outsort
4314 if exist('outpercent')
4315     outsort = { outsort outpercent };
4316 end;
4317 
4318 fprintf('Done.\n\n');
4319 
4320 %
4321 %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%  End erpimage() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4322 %
4323 if strcmpi(noshow, 'no')
4324     axes('position',gcapos);
4325     axis off
4326 end;
4327 warning on;
4328 
4329 return
4330 %
4331 %% %%%%%%%%%%%%%%%%% function plot1trace() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4332 %
4333 function [plot_handle] = plot1trace(ax,times,trace,axlimits,signif,stdev,winlocs,erp_grid,erp_vltg_ticks)
4334 %function [plot_handle] = plot1trace(ax,times,trace,axlimits,signif,stdev,winlocs,erp_grid,erp_vltg_ticks)
4335 %                           If signif is [], plot trace +/- stdev
4336 %                           Else if signif, plot trace and signif(1,:)&signif(2,:) fill.
4337 %                           Else, plot trace alone.
4338 %                           If winlocs not [], plot grey back image(s) in sort window
4339 %                                       winlocs(1,1)-> winlocs(1,end) (ms)
4340 %                                        ...
4341 %                                       winlocs(end,1)-> winlocs(end,end) (ms)
4342 
4343 FILLCOLOR    = [.66 .76 1];
4344 WINFILLCOLOR    = [.88 .92 1];
4345 ERPDATAWIDTH = 2;
4346 ERPZEROWIDTH = 2;
4347 axes(ax);
4348 
4349 if nargin<9,
4350     erp_vltg_ticks=[]; %if erp_vltg_ticks is not empty, those will be the ticks used on the y-axis (voltage axis for ERPs)
4351 end
4352 
4353 if ~isempty(winlocs)
4354     for k=1:size(winlocs,1)
4355         winloc = winlocs(k,:);
4356         fillwinx = [winloc winloc(end:-1:1)];
4357         hannwin = makehanning(length(winloc));
4358         hannwin = hannwin./max(hannwin); % make max = 1
4359         hannwin = hannwin(:)'; % make row vector
4360         if ~isempty(axlimits) & sum(isnan(axlimits))==0
4361             % fillwiny = [repmat(axlimits(3),1,length(winloc)) repmat(axlimits(4),1,length(winloc))];
4362             fillwiny = [hannwin*axlimits(3) hannwin*axlimits(4)];
4363         else
4364             % fillwiny = [repmat(min(trace)*1.1,1,length(winloc)) repmat(max(trace)*1.1,1,length(winloc))];
4365             fillwiny = [hannwin*2*min(trace) hannwin*2*max(trace)];
4366         end
4367         fillwh = fill(fillwinx,fillwiny, WINFILLCOLOR); hold on    % plot 0+alpha
4368         set(fillwh,'edgecolor',WINFILLCOLOR-[.00 .00 0]); % make edges NOT highlighted
4369     end
4370 end
4371 if ~isempty(signif);% (2,times) array giving upper and lower signif limits
4372     filltimes = [times times(end:-1:1)];
4373     if size(signif,1) ~=2 | size(signif,2) ~= length(times)
4374         fprintf('plot1trace(): signif array must be size (2,frames)\n')
4375         return
4376     end
4377     fillsignif = [signif(1,:) signif(2,end:-1:1)];
4378     fillh = fill(filltimes,fillsignif, FILLCOLOR); hold on    % plot 0+alpha
4379     set(fillh,'edgecolor',FILLCOLOR-[.02 .02 0]); % make edges slightly highlighted
4380     % [plot_handle] = plot(times,signif, 'r','LineWidth',1); hold on    % plot 0+alpha
4381     % [plot_handle] = plot(times,-1*signif, 'r','LineWidth',1); hold on % plot 0-alpha
4382 end
4383 if ~isempty(stdev)
4384     [st1] = plot(times,trace+stdev, 'r--','LineWidth',1); hold on % plot trace+stdev
4385     [st2] = plot(times,trace-stdev, 'r--','LineWidth',1); hold on % plot trace-stdev
4386 end
4387 %linestyles={'r','m','c','b'};
4388 % 'LineStyle',linestyles{traceloop}); hold on
4389 linecolor={[0 0 1],[.25 0 .75],[.75 0 .25],[1 0 0]};
4390 plot_handle=zeros(1,size(trace,1));
4391 for traceloop=1:size(trace,1),
4392     [plot_handle(traceloop)] = plot(times,trace(traceloop,:),'LineWidth',ERPDATAWIDTH, ...
4393         'color',linecolor{traceloop}); hold on
4394 end
4395 
4396 %Assume that multiple traces are equally sized divisions of data
4397 switch size(trace,1),
4398     case 2
4399         legend('Lower 50%','Higher 50%');
4400     case 3
4401         legend('Lowest 33%','Middle 33%','Highest 33%');
4402     case 4
4403         legend('Lowest 25%','2nd Lowest 25%','3rd Lowest 25%','Highest 25%');
4404 end
4405     
4406 if ~isempty(axlimits) && sum(isnan(axlimits))==0
4407     if axlimits(2)>axlimits(1) && axlimits(4)>axlimits(3)
4408         axis([axlimits(1:2) 1.1*axlimits(3:4)])
4409     end
4410     l1=line([axlimits(1:2)],[0 0],    'Color','k',...
4411         'linewidth',ERPZEROWIDTH); % y=zero-line
4412     timebar=0; 
4413     l2=line([1 1]*timebar,axlimits(3:4)*1.1,'Color','k',...
4414         'linewidth',ERPZEROWIDTH); % x=zero-line
4415 
4416     %y-ticks
4417     if isempty(erp_vltg_ticks),
4418         shrunk_ylimits=axlimits(3:4)*.8;
4419         ystep=(shrunk_ylimits(2)-shrunk_ylimits(1))/4;
4420         if ystep>1,
4421             ystep=round(ystep);
4422         else
4423             ord=orderofmag(ystep);
4424             ystep=round(ystep/ord)*ord;
4425         end
4426         if (sign(shrunk_ylimits(2))*sign(shrunk_ylimits(1)))==1, %y shrunk_ylimits don't include 0
4427             erp_yticks=shrunk_ylimits(1):ystep:shrunk_ylimits(2);
4428         else %y limits include 0
4429             % erp_yticks=0:ystep:shrunk_ylimits(2); %ensures 0 is a tick point
4430             % erp_yticks=[erp_yticks [0:-ystep:shrunk_ylimits(1)]];
4431             erp_yticks=0:ystep:axlimits(4); %ensures 0 is a tick point
4432             erp_yticks=[erp_yticks [0:-ystep:axlimits(3)]];
4433         end
4434         if erp_grid,
4435             set(gca,'ytick',unique(erp_yticks),'ygrid','on');
4436         else
4437             set(gca,'ytick',unique(erp_yticks));
4438         end
4439     else    
4440         set(gca,'ytick',erp_vltg_ticks); 
4441     end
4442 end
4443 
4444 %make ERP traces on top
4445 kids=get(gca,'children')';
4446 for hndl_loop=plot_handle,
4447    id=find(kids==hndl_loop);
4448    kids(id)=[];
4449    kids=[hndl_loop kids];
4450 end
4451 set(gca,'children',kids');
4452 plot_handle=[plot_handle l1 l2];
4453 
4454 %
4455 %% %%%%%%%%%%%%%%%%% function phasedet() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4456 %
4457 % phasedet() - function used in erpimage.m
4458 %              Constructs a complex filter at frequency freq
4459 %
4460 function [ang,amp,win] = phasedet(data,frames,srate,nwin,freq)
4461 % Typical values:
4462 %   frames = 768;
4463 %   srate = 256; % Hz
4464 %   nwin = [200:300];
4465 %   freq = 10; % Hz
4466 
4467 data = reshape(data,[frames prod(size(data))/frames]);
4468 % number of cycles depends on window size
4469 % number of cycles automatically reduced if smaller window
4470 % Note: as the number of cycles changes, the frequency shifts
4471 %       a little -- this should be fixed
4472 
4473 win = exp(2i*pi*freq(:)*[1:length(nwin)]/srate);
4474 win = win .* repmat(makehanning(length(nwin))',length(freq),1);
4475 
4476 %tmp =gcf; figure; plot(real(win)); figure(tmp);
4477 %fprintf('ANY NAN ************************* %d\n', any(any(isnan( data(nwin,:)))));
4478 
4479 tmpdata = data(nwin,:) - repmat(mean(data(nwin,:), 1), [size(data,1) 1]);
4480 resp = win * tmpdata;
4481 ang = angle(resp);
4482 amp = abs(resp);
4483 
4484 %
4485 %% %%%%%%%%%%%% function prctle() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4486 %
4487 function prctl = prctle(data,pc); % return percentile of a distribution
4488 [prows pcols] = size(pc);
4489 if prows ~= 1 & pcols ~= 1
4490     error('pc must be a scalar or a vector.');
4491 end
4492 if any(pc > 100) | any(pc < 0)
4493     error('pc must be between 0 and 100');
4494 end
4495 [i,j] = size(data);
4496 sortdata = sort(data);
4497 if i==1 | j==1 % if data is a vector
4498     i = max(i,j); j = 1;
4499     if i == 1,
4500         fprintf('  prctle() note: input data is a single scalar!\n')
4501         y = data*ones(length(pc),1); % if data is scalar, return it
4502         return;
4503     end
4504     sortdata = sortdata(:);
4505 end
4506 pt = [0 100*((1:i)-0.5)./i 100];
4507 sortdata = [min(data); sortdata; max(data)];
4508 prctl = interp1(pt,sortdata,pc);
4509 %
4510 %% %%%%%%%%%%%%%%%%%%%%% function nan_mean() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4511 %
4512 % nan_mean() - Take the column means of a matrix, ignoring NaN values.
4513 %              Return significance bounds if alpha (0 < alpha< <1) is given.
4514 %
4515 function [out, outalpha]  = nan_mean(in,alpha)
4516 NPERM = 500;
4517 
4518 intrials = size(in,1);
4519 inframes = size(in,2);
4520 nans = find(isnan(in));
4521 in(nans) = 0;
4522 sums = sum(in);
4523 nonnans = ones(size(in));
4524 nonnans(nans) = 0;
4525 nonnans = sum(nonnans);
4526 nononnans = find(nonnans==0);
4527 nonnans(nononnans) = 1;
4528 out = sum(in)./nonnans;
4529 outalpha = [];
4530 if nargin>1
4531     if NPERM < round(3/alpha)
4532         NPERM = round(3/alpha);
4533     end
4534     fprintf('Performing a permuration test using %d permutations to determine ERP significance thresholds... ',NPERM);
4535     permerps = zeros(NPERM,inframes);
4536     for n=1:NPERM
4537         signs = sign(randn(1,intrials)'-0.5);
4538         permerps(n,:) = sum(repmat(signs,1,inframes).*in)./nonnans;
4539         if ~rem(n,50)
4540             fprintf('%d ',n);
4541         end
4542     end
4543     fprintf('\n');
4544     permerps = sort(abs(permerps));
4545     alpha_bnd = floor(2*alpha*NPERM); % two-sided probability threshold
4546     outalpha = permerps(end-alpha_bnd,:);
4547 end
4548 out(nononnans) = NaN;
4549 %
4550 %% %%%%%%%%%%%%%%%%%%%%% function nan_std() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4551 %
4552 function out = nan_std(in)
4553 
4554 nans = find(isnan(in));
4555 in(nans) = 0;
4556 
4557 nonnans = ones(size(in));
4558 nonnans(nans) = 0;
4559 nonnans = sum(nonnans);
4560 nononnans = find(nonnans==0);
4561 nonnans(nononnans) = 1;
4562 
4563 out = sqrt((sum(in.^2)-sum(in).^2./nonnans)./(nonnans-1));
4564 out(nononnans) = NaN;
4565 
4566 % symmetric hanning function
4567 function w = makehanning(n)
4568 if ~rem(n,2)
4569     w = 0.5*(1 - cos(2*pi*(1:n/2)'/(n+1)));
4570     w = [w; w(end:-1:1)];
4571 else
4572     w = 0.5*(1 - cos(2*pi*(1:(n+1)/2)'/(n+1)));
4573     w = [w; w(end-1:-1:1)];
4574 end
4575 %
4576 %% %%%%%%%%%%%%%%%%%%%%% function compute_percentile() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4577 %
4578 function outpercent = compute_percentile(sortvar, percent, outtrials, winsize);
4579 ntrials = length(sortvar);
4580 outtrials=round(outtrials);
4581 sortvar = [ sortvar sortvar sortvar ];
4582 winvals = [round(-winsize/2):round(winsize/2)];
4583 outpercent = zeros(size(outtrials));
4584 for index = 1:length(outtrials)
4585     sortvarval = sortvar(outtrials(index)+ntrials+winvals);
4586     sortvarval = sort(sortvarval);
4587     outpercent(index) = sortvarval(round((length(winvals)-1)*percent)+1);
4588 end;
4589 
4590 %
4591 %% %%%%%%%%%%%%%%%%%%%%% function orderofmag() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4592 %
4593 function ord=orderofmag(val)
4594 %function ord=orderofmag(val)
4595 %
4596 % Returns the order of magnitude of the value of 'val' in multiples of 10
4597 % (e.g., 10^-1, 10^0, 10^1, 10^2, etc ...)
4598 % used for computing erpimage trial axis tick labels as an alternative for
4599 % plotting sorting variable
4600 
4601 val=abs(val);
4602 if val>=1
4603     ord=1;
4604     val=floor(val/10);
4605     while val>=1,
4606         ord=ord*10;
4607         val=floor(val/10);
4608     end
4609     return;
4610 else
4611     ord=1/10;
4612     val=val*10;
4613     while val<1,
4614         ord=ord/10;
4615         val=val*10;
4616     end
4617     return;
4618 end
4619 
4620 %
4621 %% %%%%%%%%%%%%%%%%%%%%% function find_crspnd_pt() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4622 %
4623 function y_pt=find_crspnd_pt(targ,vals,outtrials)
4624 %function id=find_crspnd_pt(targ,vals,outtrials)
4625 %
4626 % Inputs:
4627 %   targ      - desired value of sorting variable
4628 %   vals      - a vector of observed sorting variables (possibly smoothed)
4629 %   outtrials - a vector of y-axis values corresponding to trials in the
4630 %               ERPimage (this will just be 1:n_trials if there's no
4631 %               smoothing)
4632 %
4633 % Output:
4634 %   y_pt  - y-axis value (in units of trials) corresponding to "targ".
4635 %          If "targ" matches more than one y-axis pt, the median point is
4636 %          returned.  If "targ" falls between two points, y_pt is linearly
4637 %          interpolated.
4638 %
4639 % Note: targ and vals should be in the same units (e.g., milliseconds)
4640 
4641     %find closest point above
4642     abv=find(vals>=targ);
4643     if isempty(abv),
4644         %point lies outside of vals range, can't interpolate
4645         y_pt=[];
4646         return
4647     end
4648     abv=abv(1);
4649     
4650     %find closest point below
4651     blw=find(vals<=targ);
4652     if isempty(blw),
4653         %point lies outside of vals range, can't interpolate
4654         y_pt=[];
4655         return
4656     end
4657     blw=blw(end);
4658     
4659     if (vals(abv)==vals(blw)),
4660         %exact match
4661         ids=find(vals==targ);
4662         y_pt=median(outtrials(ids));
4663     else
4664         %interpolate point
4665         
4666         %lst squares inear regression
4667         B=regress([outtrials(abv) outtrials(blw)]',[ones(2,1) [vals(abv) vals(blw)]']);
4668         
4669         %predict outtrial point from target value
4670         y_pt=[1 targ]*B;
4671         
4672     end
4673 
4674