erpimages

 erpimages() - Plots (a) ERPimages from two different bins/channels/eim files, 
               (b) the difference between those ERPimages, and (c) the ERPs
               corresponding to those ERPimages in a single figure. Click 
               on individual plots to examine them separately and zoom 
               (using axcopy()). Based on the function erpimage.m.             


 Usage:
  >> figure; [outdata,outsort,outtrials] = erpimages(data1,data2,sortvar1,sortvar2,chan1,chan2,varargin)

  
 Required Inputs:
   data1     = [vector or matrix] Single-channel input data from chan1 to image.
               Format (time points,trials)
   data2     = [vector or matrix] Single-channel input data from chan2 to image.
               Format (time points,trials)
   sortvar1  = [vector] Variable to sort data1 epochs on (length(sortvar) = # of trials in data1)
               Example: sortvar may by subject response time in each epoch (in msec).            
   sortvar2  = [vector] Variable to sort data2 epochs on (length(sortvar) = # of trials in data2)
   chan1     = [string] the name of channel corresponding to data1 (e.g., 'RLOc')
   chan2     = [string] the name of channel corresponding to data2 (e.g., 'LLOc')

 Unordered Optional Inputs:
   times         = [vector | []] vector of latencies (in msec) for each epoch time point.
                   Else [] -> 1:# of time points {default: []}
   title         = [string] Figure title {default: none}
   img1_title    = [string] title of the data1 ERPimage {default: 'Img1'}
   img2_title    = [string] title of the data2 ERPimage {default: 'Img2'}
   stdev         = [scalar>0] The standard deviation (in units of epochs) 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: 1/7, no smoothing}
   decfactor     = Factor to decimate|interpolate ntrials by (may be non-integer)
                   Else, if this is large (> sqrt(# of epochs)), output this many epochs.
                   {default: no decimation} 
   cbarlimits    = [min max] Two element vector specifying the minimum and
                   maximum value of the ERPimage color scale. {default:
                   -/+ max(abs(data)))}
   cbartitle     = [string] Title of color scale {default: '\muV'}
   yimglabel     = [string] Label for ERPimage y-axes.  If not 'Trials', yticks are automatically
                   set to units of sorting variable. {default: 'Trials'}
   yimgticks     = [vector] Points on ERPimage y-axes at which to make
                   tick marks. Provide values in units of y-axes label (either 'Trials' or
                   sorting variable). {default: MATLAB default ticks}
   timelimits    = [min max] Two element vector specifying the minimum and
                   maximum value of all time axes (in ms). {default: [min(times)
                   max(times)]
   erplimits     = [min max] Two element vector specifying the minimum and
                   maximum value (in uV) of the voltage axis for data1's and 
                   data2's ERPs. {default: -/+ max(max(abs([erp1 erp2])))}
   erpdiflimits  = [min max] Two element vector specifying the minimum and
                   maximum value (in uV) of the voltage axis for the difference 
                   wave between data1 and data2. {default: -/+ max(max(abs([erp1-erp2])))}
   sortvarlimits = [min max] Two element vector specifying the minimum and
                   maximum value (in units of 'Trials' or the sorting variable, depending on 
                   yimglabel) of the ERPimage y-axis.  ERPs are computed
                   only from epochs that fall within this range.
                   {default: ERPimage min and max}
   yerplabel     = The label for the voltage axis on ERP plots {default:
                   'ERP (\muV)'}
   erpgrid       = ['on' | 'off'] If 'on', horizontal dashed lines will
                   be plot at voltage axis tick marks on ERP plots. {default: 'on'}
   topos         = ['on' | 'off'] If 'on', two cartoon heads will appear
                   above ERPimages to illustrate locations of chan1 and
                   chan2.  'chanlocs1' and 'chanlocs2' (see below) must be
                   specified to use this option. {default: 'on'}               
   chanlocs1     = EEGLAB chanlocs struct that specifies electrode
                   coordinates for chan1. {default: none}
   chanlocs2     = EEGLAB chanlocs struct that specifies electrode
                   coordinates for chan2. {default: chanlocs2=chanlocs1}
   renorm        = ['on'|'off'| formula] Normalize sorting variable to epoch
                   time range. 'on'= autoscale. Formula must be a linear 
                   transformation in the format 'a*x+b' (e.g., '3*x+2'). 
                   {default: 'off'}
   rmerp         = ['on'|'off'] Subtract the average ERP from each trial before
                    processing {default: 'off'}
   showsortvar   = ['on'|'off'] Plot sortvar on top of ERPimages {default: 'on'}
   marktimes     = [vector] Plot vertical dashed lines at specified latencies (in ms).
                   erpimage.m's "vert" option. {default: ignored}
   marktrials    = [vector] Plot horizontal lines at specified ERPimage y-axis
                   values.  Depending on value of 'yimglabel' argument (above),
                   vector should be in units of either 'Trials' or the
                   sorting variable.  Analogous to erpimage.m's "horz" option.
                   {default: ignored}
   srate         = [positive scalar] sampling rate. Necessary for 'filt'
                   option (below). {default: none}
   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 high pass filter.  If high_boundary=srate/2, then 
                   the filter is a low 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 (see below) 
                   as well since the mean prestimulus baseline may no longer be 0 after 
                   the filter is applied. {default: no filtering}
   baseline      = [low_boundary high_boundary] a time window (in msec) whose mean 
                   amplitude in each epoch will be removed from each epoch (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}


 Optional outputs:
    outdata  = (times,epochsout,img) data matrix (after smoothing).
               img=1 corresponds to data1. img=2 corresponds to data2. img=3 is 
               data1-data2 (i.e., the difference ERPimage).
     outsort = (1,epochsout) smoothed sorting variable values (corresponding to each row of the ERPimage).
   outtrials = (1,epochsout) smoothed trial numbers

 Example: 
 >>load eeglab_data_epochs_ica.set -MAT % this file should be included with EEGLAB 
 >>epochs_with_rt=setdiff(1:80,[4 27 46 71 76]);
 >>data1=squeeze(EEG.data(25,:,epochs_with_rt));
 >>data2=squeeze(EEG.data(29,:,epochs_with_rt));
 >>rt=[]; for a=epochs_with_rt, rt=[rt EEG.epoch(a).eventlatency{end}]; end;
 >>figure; erpimages(data1,data2,rt,rt,'PO7','PO8','chanlocs1',EEG.chanlocs,'times',EEG.times,'stdev',1);

 Author: 
 David Groppe
 Kutaslab, 11/2009 

 But largely based on erpimage.m function written by Scott Makeig, Tzyy-Ping Jung 
 & Arnaud Delorme, (CNL/Salk Institute)

 Notes: 
 -Epochs with identical sorting variable values are replaced with their
 mean. Thus, if there are such "tied" epochs, there is smoothing in addition
 to the moving average window.  #/% of epochs that tie other epochs are reported
 in the MATLAB command line window.

0001 % erpimages() - Plots (a) ERPimages from two different bins/channels/eim files,
0002 %               (b) the difference between those ERPimages, and (c) the ERPs
0003 %               corresponding to those ERPimages in a single figure. Click
0004 %               on individual plots to examine them separately and zoom
0005 %               (using axcopy()). Based on the function erpimage.m.
0006 %
0007 %
0008 % Usage:
0009 %  >> figure; [outdata,outsort,outtrials] = erpimages(data1,data2,sortvar1,sortvar2,chan1,chan2,varargin)
0010 %
0011 %
0012 % Required Inputs:
0013 %   data1     = [vector or matrix] Single-channel input data from chan1 to image.
0014 %               Format (time points,trials)
0015 %   data2     = [vector or matrix] Single-channel input data from chan2 to image.
0016 %               Format (time points,trials)
0017 %   sortvar1  = [vector] Variable to sort data1 epochs on (length(sortvar) = # of trials in data1)
0018 %               Example: sortvar may by subject response time in each epoch (in msec).
0019 %   sortvar2  = [vector] Variable to sort data2 epochs on (length(sortvar) = # of trials in data2)
0020 %   chan1     = [string] the name of channel corresponding to data1 (e.g., 'RLOc')
0021 %   chan2     = [string] the name of channel corresponding to data2 (e.g., 'LLOc')
0022 %
0023 % Unordered Optional Inputs:
0024 %   times         = [vector | []] vector of latencies (in msec) for each epoch time point.
0025 %                   Else [] -> 1:# of time points {default: []}
0026 %   title         = [string] Figure title {default: none}
0027 %   img1_title    = [string] title of the data1 ERPimage {default: 'Img1'}
0028 %   img2_title    = [string] title of the data2 ERPimage {default: 'Img2'}
0029 %   stdev         = [scalar>0] The standard deviation (in units of epochs) of the Gaussian
0030 %                   window used to smooth (vertically) with a moving-average.  Gaussian
0031 %                   window extends three standard deviations below and three standard
0032 %                   deviations above window center (trials beyond window are not incorporated
0033 %                   into average). {default: 1/7, no smoothing}
0034 %   decfactor     = Factor to decimate|interpolate ntrials by (may be non-integer)
0035 %                   Else, if this is large (> sqrt(# of epochs)), output this many epochs.
0036 %                   {default: no decimation}
0037 %   cbarlimits    = [min max] Two element vector specifying the minimum and
0038 %                   maximum value of the ERPimage color scale. {default:
0039 %                   -/+ max(abs(data)))}
0040 %   cbartitle     = [string] Title of color scale {default: '\muV'}
0041 %   yimglabel     = [string] Label for ERPimage y-axes.  If not 'Trials', yticks are automatically
0042 %                   set to units of sorting variable. {default: 'Trials'}
0043 %   yimgticks     = [vector] Points on ERPimage y-axes at which to make
0044 %                   tick marks. Provide values in units of y-axes label (either 'Trials' or
0045 %                   sorting variable). {default: MATLAB default ticks}
0046 %   timelimits    = [min max] Two element vector specifying the minimum and
0047 %                   maximum value of all time axes (in ms). {default: [min(times)
0048 %                   max(times)]
0049 %   erplimits     = [min max] Two element vector specifying the minimum and
0050 %                   maximum value (in uV) of the voltage axis for data1's and
0051 %                   data2's ERPs. {default: -/+ max(max(abs([erp1 erp2])))}
0052 %   erpdiflimits  = [min max] Two element vector specifying the minimum and
0053 %                   maximum value (in uV) of the voltage axis for the difference
0054 %                   wave between data1 and data2. {default: -/+ max(max(abs([erp1-erp2])))}
0055 %   sortvarlimits = [min max] Two element vector specifying the minimum and
0056 %                   maximum value (in units of 'Trials' or the sorting variable, depending on
0057 %                   yimglabel) of the ERPimage y-axis.  ERPs are computed
0058 %                   only from epochs that fall within this range.
0059 %                   {default: ERPimage min and max}
0060 %   yerplabel     = The label for the voltage axis on ERP plots {default:
0061 %                   'ERP (\muV)'}
0062 %   erpgrid       = ['on' | 'off'] If 'on', horizontal dashed lines will
0063 %                   be plot at voltage axis tick marks on ERP plots. {default: 'on'}
0064 %   topos         = ['on' | 'off'] If 'on', two cartoon heads will appear
0065 %                   above ERPimages to illustrate locations of chan1 and
0066 %                   chan2.  'chanlocs1' and 'chanlocs2' (see below) must be
0067 %                   specified to use this option. {default: 'on'}
0068 %   chanlocs1     = EEGLAB chanlocs struct that specifies electrode
0069 %                   coordinates for chan1. {default: none}
0070 %   chanlocs2     = EEGLAB chanlocs struct that specifies electrode
0071 %                   coordinates for chan2. {default: chanlocs2=chanlocs1}
0072 %   renorm        = ['on'|'off'| formula] Normalize sorting variable to epoch
0073 %                   time range. 'on'= autoscale. Formula must be a linear
0074 %                   transformation in the format 'a*x+b' (e.g., '3*x+2').
0075 %                   {default: 'off'}
0076 %   rmerp         = ['on'|'off'] Subtract the average ERP from each trial before
0077 %                    processing {default: 'off'}
0078 %   showsortvar   = ['on'|'off'] Plot sortvar on top of ERPimages {default: 'on'}
0079 %   marktimes     = [vector] Plot vertical dashed lines at specified latencies (in ms).
0080 %                   erpimage.m's "vert" option. {default: ignored}
0081 %   marktrials    = [vector] Plot horizontal lines at specified ERPimage y-axis
0082 %                   values.  Depending on value of 'yimglabel' argument (above),
0083 %                   vector should be in units of either 'Trials' or the
0084 %                   sorting variable.  Analogous to erpimage.m's "horz" option.
0085 %                   {default: ignored}
0086 %   srate         = [positive scalar] sampling rate. Necessary for 'filt'
0087 %                   option (below). {default: none}
0088 %   filt          = [low_boundary high_boundary] a two element vector indicating
0089 %                   the frequency cut-offs for a 3rd order Butterworth filter that
0090 %                   will be applied to each trial of data.  If low_boundary=0, the
0091 %                   filter is a high pass filter.  If high_boundary=srate/2, then
0092 %                   the filter is a low pass filter.  If both boundaries are between
0093 %                   0 and srate/2, then the filter is a bandpass filter. If both
0094 %                   boundaries are between 0 and -srate/2, then the filter is a bandstop
0095 %                   filter (with boundaries equal to the absolute values of low_boundary
0096 %                   and high_boundary).  Note, using this option requires the 'srate'
0097 %                   option to be specified and the signal processing toolbox function
0098 %                   butter.m.  You should probably use the 'baseline' option (see below)
0099 %                   as well since the mean prestimulus baseline may no longer be 0 after
0100 %                   the filter is applied. {default: no filtering}
0101 %   baseline      = [low_boundary high_boundary] a time window (in msec) whose mean
0102 %                   amplitude in each epoch will be removed from each epoch (e.g.,
0103 %                   [-100 0]) after filtering. Useful in conjunction with 'filt' option
0104 %                   to re-basline trials after they have been filtered. Not necessary
0105 %                   if data have already been baselined and ERPimage processing does not
0106 %                   affect baseline amplitude {default: no further baselining of data}
0107 %
0108 %
0109 % Optional outputs:
0110 %    outdata  = (times,epochsout,img) data matrix (after smoothing).
0111 %               img=1 corresponds to data1. img=2 corresponds to data2. img=3 is
0112 %               data1-data2 (i.e., the difference ERPimage).
0113 %     outsort = (1,epochsout) smoothed sorting variable values (corresponding to each row of the ERPimage).
0114 %   outtrials = (1,epochsout) smoothed trial numbers
0115 %
0116 % Example:
0117 % >>load eeglab_data_epochs_ica.set -MAT % this file should be included with EEGLAB
0118 % >>epochs_with_rt=setdiff(1:80,[4 27 46 71 76]);
0119 % >>data1=squeeze(EEG.data(25,:,epochs_with_rt));
0120 % >>data2=squeeze(EEG.data(29,:,epochs_with_rt));
0121 % >>rt=[]; for a=epochs_with_rt, rt=[rt EEG.epoch(a).eventlatency{end}]; end;
0122 % >>figure; erpimages(data1,data2,rt,rt,'PO7','PO8','chanlocs1',EEG.chanlocs,'times',EEG.times,'stdev',1);
0123 %
0124 % Author:
0125 % David Groppe
0126 % Kutaslab, 11/2009
0127 %
0128 % But largely based on erpimage.m function written by Scott Makeig, Tzyy-Ping Jung
0129 % & Arnaud Delorme, (CNL/Salk Institute)
0130 %
0131 % Notes:
0132 % -Epochs with identical sorting variable values are replaced with their
0133 % mean. Thus, if there are such "tied" epochs, there is smoothing in addition
0134 % to the moving average window.  #/% of epochs that tie other epochs are reported
0135 % in the MATLAB command line window.
0136 
0137 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%%
0138 %
0139 %12-7-09 Function did not take into account unused trials (NaN values for
0140 %sorting variable) when comparing data1 and data2.  Thus it would not let
0141 %you compare two bins at the same channel since it thought you were trying
0142 %to compare the exact same two ERPimages.  This has been fixed now (thanks
0143 %Priya!).
0144 
0145 %%%%%%%%%%%%%%%% POSSIBLE FUTURE DEVELOPMENT %%%%%%%%%%%%%%%%
0146 %-lots of erpimage.m functions are not currently implemented but might be useful.
0147 % For example:
0148 %  Timewarping
0149 %  Coherence
0150 %  Spectra
0151 %  Align time
0152 %  Auxvar
0153 %  Ampsort
0154 %  ERPalpha
0155 %  avg_type: (Gaussian always assumed)
0156 %  nosort: (not necessary, just don't specify sortvar)
0157 %  cbar (probably not useful)
0158 %  erpstd (probably not useful)
0159 %  TIMEX (time is always assumed to be x-axis)
0160 %  plotmode
0161 %  NoTimeFlag
0162 %  coherfreq
0163 % Note also that function "prctle" was removed since it was no longer being
0164 % used
0165 %
0166 %-Would it make more sense to more sortvar renorming until after is smoothed with the moving
0167 %avg?  Currently it is renormed before smoothing because that's how erpimage.m does it.
0168 %
0169 %-If topos are not plot (i.e., 'topos' is 'off'), figure background is grey.  Otherwise,
0170 %figure background is EEGLAB blue.  It would be nice to have it be consistent.
0171 %
0172 %-Implement verblevel to control amount of run-time reporting done by
0173 %function
0174 
0175 %% Function Outline
0176 % I. Error Check Arguments
0177 % II. Pre-Process Sortvar
0178 % III. Construct Gaussian window to weight trials
0179 % IV. Compute ERPimage y-axis tick values and labels (if requested)
0180 % V. First EEG data loop: process ERPimage data to get max & min
0181 % (e.g., frequency domain filtering, averaging across epochs)
0182 % VI. Find color bar axes limits
0183 % VII. Image the aligned/sorted/smoothed data
0184 % VIII. Plot Topos/electrode locations
0185 % IX. Plot ERPs
0186 %
0187 % Embedded functions tagged onto end of code
0188 % 1. plot1trace (for plotting ERPs)
0189 % 2. nan_mean (for computing ERPs)
0190 % 3. compute_percentile (for percentile option which is currently not implemented)
0191 % 4. orderofmag
0192 % 5. function find_crspnd_pt()
0193 % 6. watchit
0194 
0195 
0196 function [outdata,outsort,outtrials] = erpimages(data1,data2,sortvar1,sortvar2,chan1,chan2,varargin)
0197 
0198 
0199 %% Parse Input Arguments
0200 
0201 p=inputParser;
0202 %Note: the order of the required arguments needs to match precisely their
0203 %order in the function definition (which is also the order used by p.parse
0204 %below)
0205 p.addRequired('data1',@isnumeric);
0206 p.addRequired('data2',@isnumeric);
0207 p.addRequired('sortvar1',@isnumeric);
0208 p.addRequired('sortvar2',@isnumeric);
0209 p.addRequired('chan1',@ischar);
0210 p.addRequired('chan2',@ischar);
0211 p.addParamValue('stdev',1/7,@(x) isnumeric(x) && (x>=0));
0212 p.addParamValue('times',[],@isnumeric);
0213 p.addParamValue('decfactor',0,@isnumeric);
0214 p.addParamValue('title',[],@ischar);
0215 p.addParamValue('img1_title','Img1',@(x) ischar(x) || isempty(x));
0216 p.addParamValue('img2_title','Img2',@(x) ischar(x) || isempty(x));
0217 p.addParamValue('cbarlimits',[],@(x) (isnumeric(x) && (length(x)==2)) || isempty(x));
0218 p.addParamValue('cbartitle','\muV',@ischar);
0219 p.addParamValue('yimglabel','Trials',@ischar);
0220 p.addParamValue('yimgticks',[],@isnumeric);
0221 p.addParamValue('timelimits',[],@(x) (isnumeric(x) && (length(x)==2)) || isempty(x));
0222 p.addParamValue('erplimits',[],@(x) (isnumeric(x) && (length(x)==2)) || isempty(x));
0223 p.addParamValue('erpdiflimits',[],@(x) (isnumeric(x) && (length(x)==2)) || isempty(x));
0224 p.addParamValue('sortvarlimits',[],@(x) (isnumeric(x) && (length(x)==2)) || isempty(x));
0225 p.addParamValue('yerplabel','ERP (\muV)',@ischar);
0226 p.addParamValue('erpgrid','on',@(x) strcmpi(x,'on') || strcmpi(x,'off'));
0227 p.addParamValue('topos','on',@(x) strcmpi(x,'on') || strcmpi(x,'off'));
0228 p.addParamValue('chanlocs1',[],@isstruct); %?? maybe make a filename a legal possibility too?
0229 p.addParamValue('chanlocs2',[],@isstruct); %?? maybe make a filename a legal possibility too?
0230 p.addParamValue('renorm','off',@ischar); %check for legality of value later
0231 p.addParamValue('rmerp','off',@(x) strcmpi(x,'on') || strcmpi(x,'off'));
0232 p.addParamValue('showsortvar','on',@(x) strcmpi(x,'on') || strcmpi(x,'off'));
0233 p.addParamValue('srate',[],@(x) isnumeric(x) && (length(x)==1) && (x>0));
0234 p.addParamValue('marktimes',[],@isnumeric);
0235 p.addParamValue('marktrials',[],@isnumeric);
0236 p.addParamValue('baseline',[],@(x) (isnumeric(x) && (length(x)==2)) || isempty(x));
0237 p.addParamValue('filt',[],@(x) (isnumeric(x) && (length(x)==2)) || isempty(x));
0238 
0239 p.parse(data1,data2,sortvar1,sortvar2,chan1,chan2,varargin{:});
0240 
0241 %
0242 %%%%%%%%% I. Error Check Arguments %%%%%%%%%%%%
0243 %
0244 
0245 %Are data1 and data2 compatible?
0246 if size(data1,1)~=size(data2,1),
0247     error('data1 and data2 need to have the same number of time points.');
0248 else
0249     ntrials1=size(data1,2);
0250     ntpts=size(data1,1);
0251     ntrials2=size(data2,2);
0252     mn_ntrials=min([ntrials1 ntrials2]);
0253     not_nan_trials1=find(~isnan(sortvar1));
0254     not_nan_trials2=find(~isnan(sortvar2));
0255     if length(not_nan_trials1)==length(not_nan_trials2),
0256        if ~unique(data1(:,not_nan_trials1)-data2(:,not_nan_trials2))
0257            error('The two datasets you''re trying to contrast are exactly the same.');
0258        end
0259     end
0260 end
0261 
0262 
0263 %Are data and sortvar compatible
0264 if length(sortvar1)~=ntrials1,
0265     error('sortvar1 and data1 have to have the same number of trials.');
0266 end
0267 if length(sortvar2)~=ntrials2,
0268     error('sortvar2 and data2 have to have the same number of trials.');
0269 end
0270 
0271 
0272 %Time units & limits
0273 if isempty(p.Results.times)
0274     times = 1:ntpts;
0275     xlab='Time Points';
0276 else
0277     times=p.Results.times;
0278     if length(times)~=ntpts,
0279         error('The number time points in ''times'', %d, does not equal the number of time points in the data. %d.',length(times),ntpts)
0280     end
0281     xlab='Time (msec)';
0282 end
0283 if isempty(p.Results.timelimits)
0284     timelimits = [min(times) max(times)];
0285 else
0286     timelimits = p.Results.timelimits;
0287     if (timelimits(2)<timelimits(1)),
0288         error('First element of ''timelimits'' needs to be less than or equal to the second element.');
0289     end
0290     
0291     mn_time=min(times);
0292     if timelimits(1)<mn_time,
0293         watchit(sprintf('''timelimits'' extends before each trial''s earliest time point %d.',mn_time));
0294         fprintf('Setting earliest time limit to %d msec.\n',mn_time);
0295         timelimits(1)=mn_time;
0296         
0297         %check upper boundary as well
0298         if timelimits(2)<mn_time,
0299             timelimits(2)=min(setdiff(times,mn_time)); %second earliest time
0300             fprintf('Setting latest time limit to %d msec.\n',timelimits(2));
0301         end
0302     end
0303     
0304     mx_time=max(times);
0305     if timelimits(2)>mx_time,
0306         watchit(sprintf('''timelimits'' extends after each trial''s latest time point %d.',mx_time));
0307         fprintf('Setting latest time limit to %d msec.\n',mx_time);
0308         timelimits(2)=mx_time;
0309         
0310         %check lower boundary as well
0311         if timelimits(1)>mx_time,
0312             timelimits(1)=max(setdiff(times,mx_time)); %second latest time
0313             fprintf('Setting earliest time limit to %d msec.\n',timelimits(1));
0314         end
0315     end
0316 end
0317 if ~isempty(p.Results.erplimits),
0318     if p.Results.erplimits(1)>=p.Results.erplimits(2),
0319         error('First element of ''erplimits'' needs to be less than the second element.');
0320     end
0321 end
0322 if ~isempty(p.Results.erpdiflimits),
0323     if p.Results.erpdiflimits(1)>=p.Results.erpdiflimits(2),
0324         error('First element of ''erpdiflimits'' needs to be less than the second element.');
0325     end
0326 end
0327 
0328 %Decimation factor
0329 if p.Results.decfactor > mn_ntrials
0330     fprintf('Setting variable decfactor to max %d.\n',mn_ntrials)
0331     decfactor = mn_ntrials;
0332 else
0333     decfactor = p.Results.decfactor;
0334 end
0335 
0336 %Chanlocs
0337 if strcmpi(p.Results.topos,'on') && isempty(p.Results.chanlocs1),
0338     error('If ''topos'' is set to ''on'', you need to provide ''chanlocs'' (i.e., channel locations) as well.');
0339 end
0340 if isempty(p.Results.chanlocs2),
0341     chanlocs2=p.Results.chanlocs1;
0342 else
0343     chanlocs2=p.Results.chanlocs2;
0344 end
0345 topomap1=0; %the 1st channel to plot in a topomap
0346 for a=1:length(p.Results.chanlocs1),
0347     if strcmpi(strtok(chan1),strtok(p.Results.chanlocs1(a).labels))
0348         topomap1=a;
0349         break;
0350     end
0351 end
0352 topomap2=0; %the 2nd channel to plot in a topomap
0353 for a=1:length(chanlocs2),
0354     if strcmpi(strtok(chan2),strtok(chanlocs2(a).labels))
0355         topomap2=a;
0356         break;
0357     end
0358 end
0359 
0360 
0361 %Plotting Options
0362 img_ytick_lab=p.Results.yimgticks;
0363 verttimes=p.Results.marktimes;
0364 horzepochs=p.Results.marktrials;
0365 aligntime=NaN; % fix later ?? -DG, as is, aligning to sortvar is not an option
0366 percentiles=[]; % add later if useful ?? -DG
0367 
0368 %Other defaults from erpimage.m
0369 YES = 1;  % logical variables
0370 NO  = 0;
0371 
0372 curfig = gcf;
0373 BACKCOLOR = [0.8 0.8 0.8]; % grey background
0374 try, icadefs; catch, end;
0375 % read BACKCOLOR for plot from defs file (edit this)
0376 % read DEFAULT_SRATE for coher,phase,allamps, etc.
0377 % read YDIR for plotting ERP
0378 % Fix plotting text and line style parameters
0379 SORTWIDTH = 2.5;    % Linewidth of plotted sortvar
0380 ZEROWIDTH = 3.0;    % Linewidth of vertical 0 line
0381 VERTWIDTH = 2.5;    % Linewidth of optional vertical lines
0382 HORZWIDTH = 2.1;    % Linewidth of optional vertical lines
0383 SIGNIFWIDTH = 1.9;  % Linewidth of red significance lines for amp, coher
0384 DOTSTYLE   = 'k--'; % line style to use for vetical dotted/dashed lines
0385 LINESTYLE = '-';    % solid line
0386 LABELFONT = 10;     % font sizes for axis labels, tick labels
0387 TITLEFONT = 12; %subplot titles
0388 FIG_TITLEFONT = 16; %overall figure title
0389 TICKFONT  = 10; 
0390 
0391 PLOT_HEIGHT = 0.2;  % fraction of y dim taken up by each time series axes
0392 YGAP = 0.03;        % fraction gap between time axes
0393 YEXPAND = 1.3;      % expansion factor for y-axis about erp, amp data limits
0394 
0395 
0396 %
0397 %% %%%%%%%%%%%%%%%%% II. Pre-Process Sortvars %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0398 %
0399 
0400 
0401 
0402 %%%%%%%%%%%%%%% if sortvar(i)=NaN, remove ith trial %%%%%
0403 %for sortvar1
0404 if any(isnan(sortvar1))
0405     nantrials = find(isnan(sortvar1));
0406     fprintf('Removing data1 %d trials with NaN sortvar1 values.\n', length(nantrials));
0407     data1(:,nantrials) = [];
0408     sortvar1(nantrials) = [];
0409     ntrials1 = length(sortvar1);
0410     if ~ntrials1,
0411         error('No trials with legal values for sortvar1');
0412     end
0413 end
0414 
0415 
0416 %for sortvar2
0417 if any(isnan(sortvar2))
0418     nantrials = find(isnan(sortvar2));
0419     fprintf('Removing data2 %d trials with NaN sortvar2 values.\n', length(nantrials));
0420     data2(:,nantrials) = [];
0421     sortvar2(nantrials) = [];
0422     ntrials2 = length(sortvar2);
0423     if ~ntrials2,
0424         error('No trials with legal values for sortvar2');
0425     end
0426 end
0427 
0428 %%% find sorting variable values that appear in both bins
0429 cmn_vals=intersect(sortvar1,sortvar2); %intersect also sorts the values returned
0430 if isempty(cmn_vals),
0431    error('The two sets of data have no common values of the sorting variable.  You cannot compare them.');
0432 end
0433 cmn_data1=zeros(ntpts,mn_ntrials)*NaN;
0434 cmn_data2=zeros(ntpts,mn_ntrials)*NaN;
0435 sortedvar=zeros(mn_ntrials)*NaN;
0436 n_ties=0;
0437 dist_ties=zeros(1,length(cmn_vals));
0438 
0439 cmn_ct=0;
0440 dist_ct=0;
0441 for rept_loop=cmn_vals,
0442     ids1=find(sortvar1==rept_loop);
0443     ids2=find(sortvar2==rept_loop);
0444     n_ids1=length(ids1);
0445     n_ids2=length(ids2);
0446     mn_ids=min([n_ids1 n_ids2]);
0447     mx_ids=max([n_ids1 n_ids2]);
0448     dist_ct=dist_ct+1;
0449     if mx_ids>1,
0450         avg1=mean(data1(:,ids1),2); %avg all trials with the same value of sorting variable
0451         cmn_data1(:,cmn_ct+1:mn_ids+cmn_ct)=repmat(avg1,1,mn_ids); %add the minimum number trials to the erpimage (so that both bins have an eqal number of trials)
0452         avg2=mean(data2(:,ids2),2); %same thing for bin2
0453         cmn_data2(:,cmn_ct+1:mn_ids+cmn_ct)=repmat(avg2,1,mn_ids);
0454         sortedvar(cmn_ct+1:mn_ids+cmn_ct)=rept_loop;
0455         n_ties=n_ties+mn_ids;
0456         cmn_ct=cmn_ct+mn_ids;
0457         dist_ties(dist_ct)=mn_ids;
0458     else
0459         cmn_ct=cmn_ct+1;
0460         cmn_data1(:,cmn_ct)=data1(:,ids1);
0461         cmn_data2(:,cmn_ct)=data2(:,ids2);
0462         sortedvar(cmn_ct)=rept_loop;
0463     end
0464 end
0465 
0466 % remove NaN trials
0467 cmn_data1=cmn_data1(:,1:cmn_ct);
0468 cmn_data2=cmn_data2(:,1:cmn_ct);
0469 sortedvar=sortedvar(1:cmn_ct);
0470 fprintf('%.2f%c of the used trials (i.e., %d out of %d) have the same sortvar value as at least one other trial.\n', ...
0471     (n_ties/cmn_ct)*100,37,n_ties,cmn_ct);
0472 fprintf('Distribution of number ties per unique value of sortvar:\n');
0473 fprintf('Min: %d, 25th ptile: %d, Median: %d, 75th ptile: %d, Max: %d\n',min(dist_ties),prctile(dist_ties,25), ...
0474     median(dist_ties),prctile(dist_ties,75),max(dist_ties));
0475 fprintf('Trials with tied sorting values will be replaced by their mean.\n\n');
0476 sortidx=1:cmn_ct;
0477 
0478 
0479 
0480 %%%%%%%%%%%%%%%%%%% Renormalize sortvar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0481 switch lower(p.Results.renorm)
0482     case 'on',
0483         fprintf('\nSorting variable renormalized.\n\n');
0484         sortedvar = (sortedvar-min(sortedvar)) / (max(sortedvar) - min(sortedvar)) * ...
0485             0.5 * (max(times) - min(times)) + min(times) + 0.4*(max(times) - min(times));
0486     case 'off',;
0487     otherwise,
0488         locx = findstr('x', lower(p.Results.renorm));
0489         if length(locx) ~= 1, error('erpimage2chan: unrecognized renormalizing formula'); end;
0490         try
0491             eval( [ 'sortedvar =' p.Results.renorm(1:locx-1) 'sortedvar' p.Results.renorm(locx+1:end) ';'] );
0492         catch
0493             error('Invalid sorting variable renormalization formula.');
0494         end
0495 end;
0496 
0497 
0498 %
0499 %% III. Construct Gaussian window to weight trials
0500 %
0501 stdev=p.Results.stdev;
0502 if stdev==0,
0503    stdev=1/7; %default value
0504 end
0505 wt_wind=exp(-0.5*([-3*stdev:3*stdev]/stdev).^2)';
0506 wt_wind=wt_wind/sum(wt_wind); %normalize to unit sum
0507 avewidth=length(wt_wind);
0508 
0509 if avewidth==0,
0510     wt_wind=1;
0511     avewidth=1; %should be a window with one time point (smallest possible)
0512     watchit('Moving average smoothing window is smaller than a time point.  No smoothing will be done.');
0513 elseif avewidth > cmn_ct,
0514     watchit('Moving average window is too big for this number of trials (i.e., ''stdev'' is too big).\n');
0515     %try a smaller stdev
0516     stdev=floor((cmn_ct-1)/6);
0517     fprintf('\nTrying a Gaussian window with stdev of %d trials.\n',stdev);
0518     wt_wind=exp(-0.5*([-3*stdev:3*stdev]/stdev).^2)';
0519     wt_wind=wt_wind/sum(wt_wind); %normalize to unit sum
0520     avewidth=length(wt_wind);
0521     
0522     %recheck size
0523     if avewidth==0,
0524         wt_wind=1; %a window with one time point (smallest possible), no smoothing
0525         avewidth=1;
0526         fprintf('Attempt failed. Moving average smoothing window is smaller than a time point.  No smoothing will be done.\n');
0527     elseif avewidth > cmn_ct,
0528         fprintf('Attempt failed. Try picking a new stdev by hand.  Currently, no smoothing will be done.\n');
0529         wt_wind=1;
0530         avewidth=1;
0531     end
0532 end
0533 
0534 %smooth sorting variable if requested
0535 if avewidth > 1 || decfactor > 1
0536     fprintf('\nSmoothing the sorting variable using a window width of %g epochs ',avewidth);
0537     fprintf('\n  and a decimation factor of %g\n',decfactor);
0538     [outsort,outtrials] = movav(sortedvar,1:cmn_ct,avewidth,decfactor,[],[],wt_wind);
0539    
0540     for index=1:length(percentiles)
0541         outpercent{index} = compute_percentile( sortedvar, percentiles(index), outtrials, avewidth);
0542     end;
0543 else % don't smooth
0544     outtrials = 1:cmn_ct;
0545     outsort = sortedvar;
0546 end
0547 
0548 %
0549 %%%%%%%%%%%%%%%%%%%%% IV. Compute ERPimage y-axis tick values and labels (if requested) %%%%%%%%%%%%%%%%%%%%%%%%%%
0550 %
0551 if ~isempty(p.Results.yimglabel) && ~strcmpi(p.Results.yimglabel,'Trials')
0552     mn=min(outsort);
0553     mx=max(outsort);
0554     ord=orderofmag(mx-mn);
0555     rng_rnd=round([mn mx]/ord)*ord;
0556     if isempty(img_ytick_lab)
0557         img_ytick_lab=[rng_rnd(1):ord:rng_rnd(2)];
0558         in_range=find((img_ytick_lab>=mn) & (img_ytick_lab<=mx));
0559         img_ytick_lab=img_ytick_lab(in_range);
0560     else
0561         img_ytick_lab=unique(img_ytick_lab); %make sure it is sorted
0562         in_range=find((img_ytick_lab>=mn) & (img_ytick_lab<=mx));
0563         if length(img_ytick_lab)~=length(in_range),
0564             fprintf('\n***Warning***\n');
0565             fprintf('''img_trialax_ticks'' exceed smoothed sorting variable values. Max/min values are %f/%f.\n\n',mn,mx);
0566             img_ytick_lab=img_ytick_lab(in_range);
0567         end
0568     end
0569     n_tick=length(img_ytick_lab);
0570     img_ytick=zeros(1,n_tick);
0571     for tickloop=1:n_tick,
0572         img_ytick(tickloop)=find_crspnd_pt(img_ytick_lab(tickloop),outsort,outtrials);
0573     end
0574 end
0575 
0576 %
0577 %% V. First EEG data loop: process ERPimage data to get max & min
0578 %
0579 
0580 %pre-allocate memory
0581 erps=zeros(3,ntpts);
0582 mindat=zeros(1,3);
0583 maxdat=mindat;
0584 
0585 for eimloop=1:3,
0586     switch eimloop
0587         case 1
0588             data=cmn_data1;
0589         case 2
0590             data=cmn_data2;
0591         case 3
0592             data=cmn_data1-cmn_data2;
0593     end
0594     
0595     %
0596     %% %%%%%%%%%%%%%%  Check for NaN's in EEG %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0597     %
0598     nans = find(isnan(data));
0599     if ~isempty(nans)
0600         if eimloop==1,
0601             error('Some values of ''data1'' are NaN (i.e., not a number). There''s been some error in data processing (e.g., division by 0).');
0602         elseif eimloop==2,
0603             error('Some values of ''data2'' are NaN (i.e., not a number). There''s been some error in data processing (e.g., division by 0).');
0604         end
0605     end
0606     
0607     
0608     %
0609     %% Filter data with Butterworth filter if requested %%%%%%%%%%%%
0610     %
0611     if ~isempty(p.Results.filt),
0612         %error check
0613         if isempty(p.Results.srate),
0614             watchit('To use the ''filt'' option you need to specify the sampling rate with the ''srate'' parameter.');
0615             fprintf('Not filtering data.\n\n');
0616         else
0617             if length(p.Results.filt)~=2,
0618                 error('''filt'' parameter argument should be a two element vector.');
0619             elseif max(p.Results.filt)>(p.Results.srate/2),
0620                 error('''filt'' parameters need to be less than or equal to sampling rate/2 (i.e., %f).',p.Results.srate/2);
0621             elseif (p.Results.filt(2)==(p.Results.srate/2)) && (p.Results.filt(1)==0),
0622                 error('If second element of ''filt'' parameter is srate/2, then the first element must be greater than 0.');
0623             elseif abs(p.Results.filt(2))<=abs(p.Results.filt(1)),
0624                 error('Second element of ''filt'' parameters must be greater than first in absolute value.');
0625             elseif (p.Results.filt(1)<0) || (p.Results.filt(2)<0),
0626                 if (p.Results.filt(1)>=0) || (p.Results.filt(2)>=0),
0627                     error('BOTH parameters of ''filt'' need to be greater than or equal to zero OR need to be negative.');
0628                 end
0629                 if min(p.Results.filt)<=(-p.Results.srate/2),
0630                     error('''filt'' parameters need to be greater than sampling rate/2 (i.e., -%f) when creating a stop band.',p.Results.srate/2);
0631                 end
0632             end
0633             
0634             msg=sprintf('\nFiltering data with 3rd order Butterworth filter: ');
0635             if (p.Results.filt(1)==0),
0636                 %lowpass filter the data
0637                 [B A]=butter(3,p.Results.filt(2)*2/p.Results.srate,'low');
0638                 msg=[msg sprintf('lowpass at %.0f Hz\n',p.Results.filt(2))];
0639             elseif (p.Results.filt(2)==(p.Results.srate/2)),
0640                 %higpass filter the data
0641                 [B A]=butter(3,p.Results.filt(1)*2/p.Results.srate,'high');
0642                 msg=[msg sprintf('highpass at %.0f Hz\n',p.Results.filt(1))];
0643             elseif (p.Results.filt(1)<0)
0644                 %bandstop filter the data
0645                 [B A]=butter(3,-p.Results.filt*2/p.Results.srate,'stop');
0646                 msg=[msg sprintf('stopband from %.0f to %.0f Hz\n',-p.Results.filt(1),-p.Results.filt(2))];
0647             else
0648                 %bandpass filter the data
0649                 [B A]=butter(3,p.Results.filt*2/p.Results.srate);
0650                 msg=[msg sprintf('bandpass from %.0f to %.0f Hz\n',p.Results.filt(1),p.Results.filt(2))];
0651             end
0652             for trial=1:cmn_ct,
0653                 data(:,trial)=filtfilt(B,A,data(:,trial));
0654             end
0655             if isempty(p.Results.baseline)
0656                 msg=[msg sprintf('Note, you might want to re-baseline the data using the ''baseline'' option.\n\n')];
0657             end
0658             if eimloop==1,
0659                disp(msg); 
0660             end
0661         end
0662     end
0663     
0664     %% Mean Baseline Each Trial (if requested) %%
0665     if ~isempty(p.Results.baseline),
0666         %check argument values for errors
0667         if p.Results.baseline(2)<p.Results.baseline(1),
0668             error('First element of ''baseline'' argument needs to be less than or equal to second argument.');
0669         elseif p.Results.baseline(2)<times(1),
0670             error('Second element of ''baseline'' argument needs to be greater than or equal to epoch start time %.1f.',times(1));
0671         elseif p.Results.baseline(1)>times(end),
0672             error('First element of ''baseline'' argument needs to be less than or equal to epoch end time %.1f.',times(end));
0673         end
0674         
0675         %convert msec into time points
0676         if p.Results.baseline(1)<times(1),
0677             strt_pt=1;
0678         else
0679             strt_pt=find_crspnd_pt(p.Results.baseline(1),times,1:ntpts);
0680             strt_pt=ceil(strt_pt);
0681         end
0682         if p.Results.baseline(end)>times(end),
0683             end_pt=length(times);
0684         else
0685             end_pt=find_crspnd_pt(p.Results.baseline(2),times,1:ntpts);
0686             end_pt=floor(end_pt);
0687         end
0688         if eimloop==1,
0689             fprintf('\nRemoving pre-stimulus mean baseline from %.1f to %.1f msec.\n\n',times(strt_pt),times(end_pt));
0690         end
0691         bsln_mn=mean(data(strt_pt:end_pt,:),1);
0692         data=data-repmat(bsln_mn,ntpts,1);
0693     end
0694 
0695     if isempty(p.Results.sortvarlimits), %if sortvar limits are specified, ERPs are made later
0696         erps(eimloop,:)=mean(data'); %make ERP
0697     end
0698     
0699     %
0700     %% %%%%%%%%%%%%%%%%%%%% Remove the ERP (if requested) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0701     %
0702     if strcmpi(p.Results.rmerp, 'on')
0703         data = data - repmat(nan_mean(data')', [1 size(data,2)]);
0704     end;
0705     
0706     
0707     %
0708     %% %%%%%%%%%%%%%%%% Smooth data using Gaussian weighted moving average %%%%%%%%%%%%%%%%%%%%%%%%%%%
0709     %
0710     if avewidth > 1 || decfactor > 1
0711         if eimloop==1,
0712             fprintf('Smoothing the sorted epochs with a %g-epoch moving window.',...
0713                 avewidth);
0714             fprintf('\n');
0715             fprintf('  and a decimation factor of %g\n',decfactor);
0716         end
0717         
0718         [data,outtrials] = movav(data,1:cmn_ct,avewidth,decfactor,[],[],wt_wind);
0719         % Note: movav() here sorts using square window
0720         [outsort,outtrials] = movav(sortedvar,1:cmn_ct,avewidth,decfactor,[],[],wt_wind);
0721         
0722         for index=1:length(percentiles)
0723             outpercent{index} = compute_percentile( sortedvar, percentiles(index), outtrials, avewidth);
0724         end;
0725         if eimloop==1,
0726             fprintf('Output data will be %d ntpts by %d smoothed trials.\n',...
0727                 ntpts,length(outtrials));
0728             fprintf('Outtrials: %3.2f to %4.2f\n',min(outtrials),max(outtrials));
0729         end
0730     else % don't smooth
0731         outtrials = 1:cmn_ct;
0732         outsort = sortedvar;
0733     end
0734     outdata(:,:,eimloop)=data;
0735     mindat(eimloop) = min(min(data));
0736     maxdat(eimloop) = max(max(data));
0737     
0738 end % first EEG data processing loop
0739 
0740 %
0741 %%%%%%%%%%%%%%%%%%%%%%%%% VI. Find color bar axes limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0742 %
0743 if ~isempty(p.Results.cbarlimits)
0744     img_min = p.Results.cbarlimits(1);
0745     img_max = p.Results.cbarlimits(2);
0746     if img_min>img_max,
0747         watchit('First element of ''cbarlimits'' needs to be less than the second element.');
0748         img_min = min(min(data));
0749         img_max = max(max(data));
0750         img_max =  max(abs([img_min img_max])); % make symmetrical about 0
0751         img_min = -img_max;
0752         fprintf('The color scale imits will be the symmetric abs. data range -> [%g,%g].\n', img_min, img_max);
0753     else
0754         fprintf('Using the specified caxis range of [%g,%g].\n', img_min, img_max);
0755     end
0756     img_mxDIF=max(abs([mindat(3) maxdat(3)]));
0757     fprintf('The color scale limits for the difference between erpimages will be the symmetric abs. data range -> [%g,%g].\n', -img_mxDIF, img_mxDIF);
0758 else
0759     img_max=max(abs([mindat(1:2) maxdat(1:2)]));
0760     img_min=-img_max;
0761     fprintf('The cbarlimits for individual erpimages will be the symmetric abs. data range -> [%g,%g].\n', img_min, img_max);
0762     img_mxDIF=max(abs([mindat(3) maxdat(3)]));
0763     fprintf('The cbarlimits for the difference between erpimages will be the symmetric abs. data range -> [%g,%g].\n', -img_mxDIF, img_mxDIF);
0764 end
0765 
0766 %
0767 %% %%%%%%%%%%% VII. Image the aligned/sorted/smoothed data %%%%%%%%%%%%%%%%%%%%%%%%%%
0768 %
0769 
0770 %timelimits were set above
0771 fprintf('Data will be plotted between %g and %g ms.\n',timelimits(1),timelimits(2));
0772 
0773 for eimloop=1:3,
0774     
0775     %%%%%%%% Plot ERPimage %%%%%%%%%%
0776     switch (eimloop),
0777         case 1
0778             img_pos1=[.09 .53 .35 .35];
0779             img_ax=axes('position',img_pos1);
0780             imagesc(times,outtrials,outdata(:,:,eimloop)',[img_min img_max]);
0781         case 2
0782             img_pos2=[.54 .53 .35 .35];
0783             img_ax=axes('position',img_pos2);
0784             imagesc(times,outtrials,outdata(:,:,eimloop)',[img_min img_max]);
0785         case 3
0786             img_pos3=[.09 .09 .35 .35];
0787             img_ax=axes('position',img_pos3);
0788             imagesc(times,outtrials,outdata(:,:,eimloop)',[-img_mxDIF,img_mxDIF]); %difference between erpimages
0789     end
0790     set(gca,'Ydir','normal');
0791     axis([timelimits(1) timelimits(2) ...
0792         min(outtrials) max(outtrials)]);
0793     hold on
0794     drawnow
0795     
0796     %Y-AXIS LABEL & TICKS
0797     if eimloop~=2, %2nd ERPimage can use same ylabel as first
0798         l=ylabel(p.Results.yimglabel);
0799         set(l,'FontSize',LABELFONT);
0800         if ~isempty(p.Results.yimgticks) || (~isempty(p.Results.yimglabel) && ~strcmpi(p.Results.yimglabel,'Trials')),
0801             set(gca,'ytick',img_ytick,'yticklabel',img_ytick_lab);
0802         end
0803     elseif ~isempty(p.Results.yimglabel) && ~strcmpi(p.Results.yimglabel,'Trials'),
0804         set(gca,'ytick',img_ytick,'yticklabel',[]);
0805     end
0806     
0807     
0808     %% %%%%%%%%%%%%%%%%%%%%%%%%% plot vert lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0809     if ~isempty(verttimes)
0810         if (size(verttimes,1) ~= 1) && (size(verttimes,2) == 1) && (size(verttimes,1) ~= cmn_ct)
0811             verttimes = verttimes';
0812         end
0813         if (size(verttimes,1) ~= 1) && (size(verttimes,1) ~= cmn_ct)
0814             fprintf('\nerpimage2chan(): vert arg matrix must have 1 or %d rows\n',cmn_ct);
0815             return
0816         end;
0817         
0818         if size(verttimes,1) == 1
0819             fprintf('Plotting %d lines at times: ',size(verttimes,2));
0820         else
0821             fprintf('Plotting %d traces starting at times: ',size(verttimes,2));
0822         end
0823         for vt = verttimes % for each column
0824             fprintf('%g ',vt(1));
0825             if length(vt)==1
0826                 mydotstyle = DOTSTYLE;
0827                 if exist('auxcolors') && (length(verttimes) == length(auxcolors))
0828                     mydotstyle = auxcolors{find(verttimes == vt)};
0829                 end
0830                 figure(curfig);plot([vt vt],[0 max(outtrials)],mydotstyle,'Linewidth',VERTWIDTH);
0831             elseif length(vt)==cmn_ct
0832                 [outvt,ix] = movav(vt,1:cmn_ct,avewidth,decfactor,[],[],wt_wind);
0833                 figure(curfig);plot(outvt,outtrials,DOTSTYLE,'Linewidth',VERTWIDTH);
0834             end
0835         end
0836         fprintf('\n');
0837     end
0838     
0839     %
0840     %% %%%%%%%%%%%%%%%%%%%%%%%%% plot horizontal ('horz') lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0841     %
0842     if ~isempty(horzepochs)
0843         if size(horzepochs,1) > 1 && size(horzepochs,1) > 1
0844             watchit('Argument ''marktrials'' must be a vector.  Ignoring ''marktrials''.');
0845         else
0846             if ~isempty(p.Results.yimglabel) && ~strcmpi(p.Results.yimglabel,'Trials'),
0847                 %trial axis in units of sorting variable
0848                 mx=max(outsort);
0849                 mn=min(outsort);
0850                 fprintf('Plotting %d lines at epochs corresponding to sorting variable values: ',length(horzepochs));
0851                 for he = horzepochs % for each horizontal line
0852                     fprintf('%g ',he);
0853                     %find trial number corresponding to this value of sorting
0854                     %variable:
0855                     if (he>mn) && (he<mx)
0856                         he_ep=find_crspnd_pt(he,outsort,outtrials);
0857                         figure(curfig);plot([timelimits(1) timelimits(2)],[he_ep he_ep],LINESTYLE,'Linewidth',HORZWIDTH);
0858                     end
0859                 end
0860                 
0861             else %trial axis in units of trials
0862                 fprintf('Plotting %d lines at epochs: ',length(horzepochs));
0863                 for he = horzepochs % for each horizontal line
0864                     fprintf('%g ',he);
0865                     % overplot he on image
0866                     figure(curfig);plot([timelimits(1) timelimits(2)],[he he],LINESTYLE,'Linewidth',HORZWIDTH);
0867                 end
0868             end
0869             fprintf('\n');
0870         end
0871     end
0872     set(gca,'FontSize',TICKFONT)
0873     hold on;
0874 
0875     
0876     %
0877     %% %%%%%%%%% plot vertical line at 0 or align time %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0878     %
0879     if times(1) <= 0 & times(ntpts) >= 0
0880         figure(curfig);plot([0 0],[min(outtrials) max(outtrials)],...
0881             'k','Linewidth',ZEROWIDTH); % plot smoothed sortvar
0882     end
0883     
0884     if ( min(outsort) < timelimits(1) ...
0885             |max(outsort) > timelimits(2))
0886         ur_outsort = outsort; % store the pre-adjusted values
0887         fprintf('Not all sortvar values within time vector limits: \n')
0888         fprintf('        outliers will be shown at nearest limit.\n');
0889         i = find(outsort< timelimits(1));
0890         outsort(i) = timelimits(1);
0891         i = find(outsort> timelimits(2));
0892         outsort(i) = timelimits(2);
0893     end
0894     
0895     % ERPIMAGE TITLES
0896     switch (eimloop),
0897         case 1
0898             l=title(p.Results.img1_title);
0899             set(l,'FontSize',TITLEFONT,'fontweight','bold');
0900         case 2
0901             l=title(p.Results.img2_title);
0902             set(l,'FontSize',TITLEFONT,'fontweight','bold');
0903         case 3
0904             l=title([p.Results.img1_title '-' p.Results.img2_title]);
0905             set(l,'FontSize',TITLEFONT,'fontweight','bold');
0906     end
0907     
0908     set(gca,'Box','off');
0909     set(gca,'Fontsize',TICKFONT);
0910     set(gca,'color',BACKCOLOR);
0911     if eimloop==3,
0912         figure(curfig);l=xlabel(xlab);
0913         set(l,'Fontsize',LABELFONT);
0914     else
0915         set(gca,'xticklabel',[]);
0916     end
0917     
0918     
0919     %
0920     %% %%%%%%%%%%%%%%%%%% Overplot sortvar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0921     %
0922     if strcmpi(p.Results.showsortvar,'off')
0923         fprintf('Not overplotting sorted sortvar on data.\n');
0924     elseif isnan(aligntime) % plot sortvar on un-aligned data
0925         hold on;
0926         figure(curfig);plot(outsort,outtrials,'k','LineWidth',SORTWIDTH);
0927         drawnow
0928     end
0929     
0930     %
0931     %% %%%%%%%%%%%%% turn on axcopy() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0932     %
0933     axcopy(gcf);
0934     
0935     % Uncomment the code below if you ever decide to re-activate percentile option
0936     %     if (eimloop==1) && exist('outpercent')
0937     %         outsort = { outsort outpercent };
0938     %     end;
0939 
0940     %Change limits on ERPimage y-axis if requested
0941     if ~isempty(p.Results.sortvarlimits)
0942         sortvarlimits=p.Results.sortvarlimits;
0943         v=axis;
0944         img_mn=find_crspnd_pt(sortvarlimits(1),outsort,outtrials);
0945         if isempty(img_mn),
0946             img_mn=1;
0947             sortvarlimits(1)=outsort(1);
0948         end
0949         img_mx=find_crspnd_pt(sortvarlimits(2),outsort,outtrials);
0950         if isempty(img_mx),
0951             img_mx=length(outsort);
0952             sortvarlimits(2)=outsort(img_mx);
0953         end
0954         axis([v(1:2) img_mn img_mx]);
0955         id1=find(sortvar>=sortvarlimits(1));
0956         id2=find(sortvar<=sortvarlimits(2));
0957         id=intersect(id1,id2);
0958         fprintf('%d epochs fall within sortvar limits.\n',length(id));
0959         erps(eimloop,:)=mean(outdata(:,round(img_mn):round(img_mx),eimloop),2)';
0960     end
0961     
0962     %
0963     %% %%%%%%%%%%%%%%%%%%%%%% Plot colorbar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0964     %
0965     if (eimloop>1) %first ERPimage can use same colorbar as second ERPimage
0966         pos=get(img_ax,'Position');
0967         axcb=axes('Position',...
0968             [pos(1)+.97*pos(3)+0.02 pos(2) ...
0969             0.03 pos(4)]);
0970         if eimloop==3, %difference ERPimage
0971             figure(curfig);cbar(axcb,0,[-img_mxDIF,img_mxDIF]); % plot colorbar to right of image
0972         else
0973             figure(curfig);cbar(axcb,0,[img_min img_max]); % plot colorbar to right of image
0974         end
0975         title(p.Results.cbartitle);
0976         set(axcb,'fontsize',TICKFONT,'xtick',[]);
0977     end
0978     
0979 end %end ERPimage loop
0980 warning on;
0981 
0982 
0983 %
0984 %% %%%%%%%%%%%%% VIII. Plot Topos/electrode locations %%%%%%%%%%
0985 %
0986 if strcmpi(p.Results.topos,'on'),
0987     axes('Position',[img_pos1(1) img_pos1(2)+.99*img_pos1(4) ...
0988         .1 .1]);
0989     fprintf('Plotting topo maps to show electrode location.\n');
0990     eloc_info.plotrad = [];
0991     try
0992         topoplot(topomap1,p.Results.chanlocs1,'electrodes','off', ...
0993             'style', 'blank', 'emarkersize1chan', 10, 'chaninfo', eloc_info);
0994     catch
0995         fprintf('topoplot() plotting failed.\n');
0996     end;
0997     axis('square');
0998     
0999     axes('Position',[img_pos2(1) img_pos2(2)+.99*img_pos2(4) .1 .1]);
1000     try
1001         topoplot(topomap2,chanlocs2,'electrodes','off', ...
1002             'style', 'blank', 'emarkersize1chan',10,'emarkercolor1chan','b', 'chaninfo', eloc_info);
1003     catch
1004         fprintf('topoplot() plotting failed.\n');
1005     end;
1006     axis('square');
1007 end
1008 
1009 
1010 %
1011 %% %%%%%%%%%%%%% IX. Plot ERPs %%%%%%%%%%
1012 %
1013 if ~strcmpi(p.Results.rmerp,'on'),
1014     ERPWIDTH=1;
1015     
1016     %individual channel ERPs
1017     if isempty(p.Results.erplimits)
1018         fac = 10;
1019         maxerp = 0;
1020         while maxerp == 0
1021             maxerp = round(fac*YEXPAND*max(max(erps(1:2,:))))/fac; % minimal decimal places
1022             fac = 10*fac;
1023         end
1024         maxerp=max(maxerp);
1025         
1026         fac = 1;
1027         minerp = 0;
1028         while minerp == 0
1029             minerp = round(fac*YEXPAND*min(min(erps(1:2,:))))/fac; % minimal decimal places
1030             fac = 10*fac;
1031         end
1032         minerp=min(minerp);
1033     else
1034         minerp=p.Results.erplimits(1);
1035         maxerp=p.Results.erplimits(2);
1036     end
1037     limit = [timelimits minerp maxerp];
1038     erpAX=axes('position',[.54 .31 .35 .15]);
1039     fprintf('Plotting the ERPs\n');
1040     plot1trace(erpAX,times,erps(1:2,:),limit,[],[],[],p.Results.erpgrid,{'r','b'});
1041     erp_limits=axis;
1042     axis([timelimits erp_limits(3:4)]);
1043     set(gca,'xticklabel',[],'yaxislocation','right','Fontsize',TICKFONT);
1044     ylabel(p.Results.yerplabel);
1045     l=title(['{\color{red}' p.Results.img1_title '    \color{blue}' p.Results.img2_title '}']);
1046     set(l,'fontsize',TITLEFONT);
1047     axcopy(gcf); %clicking on axis, makes a copy of just that axis in a new figure
1048     
1049     %difference wave
1050     if isempty(p.Results.erpdiflimits)
1051         fac = 10;
1052         maxdif = 0;
1053         while maxdif == 0
1054             maxdif = round(fac*YEXPAND*max(max(erps(3,:))))/fac; % minimal decimal places
1055             fac = 10*fac;
1056         end
1057         maxdif=max(maxdif);
1058         
1059         fac = 1;
1060         mindif = 0;
1061         while mindif == 0
1062             mindif = round(fac*YEXPAND*min(min(erps(3,:))))/fac; % minimal decimal places
1063             fac = 10*fac;
1064         end
1065         mindif=min(mindif);
1066     else
1067         mindif=p.Results.erpdiflimits(1);
1068         maxdif=p.Results.erpdiflimits(2);
1069     end
1070     limit = [timelimits mindif maxdif];
1071    
1072     
1073     purple=[148 0 211]/255;
1074     difAX=axes('position',[.54 .09 .35 .15]);
1075     plot1trace(difAX,times,erps(3,:),limit,[],[],[],p.Results.erpgrid,{purple});
1076     dif_limits=axis;
1077     axis([timelimits dif_limits(3:4)]);
1078     xlabel(xlab);
1079     ylabel(p.Results.yerplabel);
1080     set(gca,'yaxislocation','right','Fontsize',TICKFONT);
1081     l=title(['\color[rgb]{.58 0 .828}' p.Results.img1_title '-' p.Results.img2_title ]);
1082     set(l,'fontsize',TITLEFONT);
1083     axcopy(gcf); %clicking on axis, makes a copy of just that axis in a new figure
1084     if ~exist('YDIR')
1085         error('Default YDIR not read from ''icadefs.m''');
1086     end
1087     if YDIR == 1
1088         set(erpAX,'ydir','normal');
1089         set(difAX,'ydir','normal');
1090     else
1091         set(erpAX,'ydir','reverse');
1092         set(difAX,'ydir','reverse');
1093     end
1094     
1095     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot vert lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1096     if ~isempty(verttimes)
1097         if size(verttimes,1) == cmn_ct
1098             vts=sort(verttimes);
1099             vts = vts(ceil(cmn_ct/2),:); % plot median verttimes values if a matrix
1100         else
1101             vts = verttimes(:)';  % make verttimes a row vector
1102         end
1103         for vt = vts
1104             if isnan(aligntime)
1105                 mydotstyle = DOTSTYLE;
1106                 if exist('auxcolors') & ...
1107                         length(verttimes) == length(verttimesColors)
1108                     mydotstyle = verttimesColors{find(verttimes == vt)};
1109                 end
1110                 axes(erpAX);plot([vt vt],erp_limits(3:4),mydotstyle,'Linewidth',ERPWIDTH);
1111                 axes(difAX);plot([vt vt],dif_limits(3:4),mydotstyle,'Linewidth',ERPWIDTH);
1112             else
1113                 axes(erpAX);plot(repmat(median(aligntime+vt-outsort),1,2),erp_limits(3:4),...
1114                     DOTSTYLE,'LineWidth',ERPWIDTH);
1115                 axes(difAX);plot(repmat(median(aligntime+vt-outsort),1,2),dif_limits(3:4),...
1116                     DOTSTYLE,'LineWidth',ERPWIDTH);
1117             end
1118         end
1119     end
1120 end %ERP plot chunk
1121 
1122 %Title for entire figure (if specified)
1123 if ~isempty(p.Results.title),
1124    fig_title=textsc(p.Results.title,'title');
1125    set(fig_title,'fontsize',FIG_TITLEFONT,'fontweight','bold');
1126 end
1127 
1128 return %end erpimage2chan main body
1129 
1130 
1131 %
1132 %% %%%%%%%%%%%%%%%%% 1. function plot1trace() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1133 %
1134 function [plot_handle] = plot1trace(ax,times,trace,axlimits,signif,stdev,winlocs,erp_grid,erp_colors)
1135 %function [plot_handle] = plot1trace(ax,times,trace,axlimits,signif,stdev,winlocs,erp_grid)
1136 %                           If signif is [], plot trace +/- stdev
1137 %                           Else if signif, plot trace and signif(1,:)&signif(2,:) fill.
1138 %                           Else, plot trace alone.
1139 %                           If winlocs not [], plot grey back image(s) in sort window
1140 %                                       winlocs(1,1)-> winlocs(1,end) (ms)
1141 %                                        ...
1142 %                                       winlocs(end,1)-> winlocs(end,end) (ms)
1143 
1144 FILLCOLOR    = [.66 .76 1];
1145 WINFILLCOLOR    = [.88 .92 1];
1146 ERPDATAWIDTH = 1;
1147 ERPZEROWIDTH = 1;
1148 axes(ax);
1149 
1150 if ~isempty(winlocs)
1151     for k=1:size(winlocs,1)
1152         winloc = winlocs(k,:);
1153         fillwinx = [winloc winloc(end:-1:1)];
1154         hannwin = makehanning(length(winloc));
1155         hannwin = hannwin./max(hannwin); % make max = 1
1156         hannwin = hannwin(:)'; % make row vector
1157         if ~isempty(axlimits) & sum(isnan(axlimits))==0
1158             % fillwiny = [repmat(axlimits(3),1,length(winloc)) repmat(axlimits(4),1,length(winloc))];
1159             fillwiny = [hannwin*axlimits(3) hannwin*axlimits(4)];
1160         else
1161             % fillwiny = [repmat(min(trace)*1.1,1,length(winloc)) repmat(max(trace)*1.1,1,length(winloc))];
1162             fillwiny = [hannwin*2*min(trace) hannwin*2*max(trace)];
1163         end
1164         fillwh = fill(fillwinx,fillwiny, WINFILLCOLOR); hold on    % plot 0+alpha
1165         set(fillwh,'edgecolor',WINFILLCOLOR-[.00 .00 0]); % make edges NOT highlighted
1166     end
1167 end
1168 if ~isempty(signif);% (2,times) array giving upper and lower signif limits
1169     filltimes = [times times(end:-1:1)];
1170     if size(signif,1) ~=2 | size(signif,2) ~= length(times)
1171         fprintf('plot1trace(): signif array must be size (2,ntpts)\n')
1172         return
1173     end
1174     fillsignif = [signif(1,:) signif(2,end:-1:1)];
1175     fillh = fill(filltimes,fillsignif, FILLCOLOR); hold on    % plot 0+alpha
1176     set(fillh,'edgecolor',FILLCOLOR-[.02 .02 0]); % make edges slightly highlighted
1177     % [plot_handle] = plot(times,signif, 'r','LineWidth',1); hold on    % plot 0+alpha
1178     % [plot_handle] = plot(times,-1*signif, 'r','LineWidth',1); hold on % plot 0-alpha
1179 end
1180 if ~isempty(stdev)
1181     [st1] = plot(times,trace+stdev, 'r--','LineWidth',1); hold on % plot trace+stdev
1182     [st2] = plot(times,trace-stdev, 'r--','LineWidth',1); hold on % plot trace-stdev
1183 end
1184 %linestyles={'r','m','c','b'};
1185 % 'LineStyle',linestyles{traceloop}); hold on
1186 plot_handle=zeros(1,size(trace,1));
1187 for traceloop=1:size(trace,1),
1188     [plot_handle(traceloop)] = plot(times,trace(traceloop,:),'LineWidth',ERPDATAWIDTH, ...
1189         'color',erp_colors{traceloop}); hold on
1190 end
1191 
1192 if ~isempty(axlimits) && sum(isnan(axlimits))==0
1193     if axlimits(2)>axlimits(1) && axlimits(4)>axlimits(3)
1194         axis([axlimits(1:2) 1.1*axlimits(3:4)])
1195     end
1196     l1=line([axlimits(1:2)],[0 0],    'Color','k',...
1197         'linewidth',ERPZEROWIDTH); % y=zero-line
1198     timebar=0;
1199     l2=line([1 1]*timebar,axlimits(3:4)*1.1,'Color','k',...
1200         'linewidth',ERPZEROWIDTH); % x=zero-line
1201     
1202     %y-ticks
1203     shrunk_ylimits=axlimits(3:4)*.8;
1204     ystep=(shrunk_ylimits(2)-shrunk_ylimits(1))/4;
1205     if ystep>1,
1206         ystep=round(ystep);
1207     else
1208         ord=orderofmag(ystep);
1209         ystep=round(ystep/ord)*ord;
1210     end
1211 
1212     if (sign(shrunk_ylimits(2))*sign(shrunk_ylimits(1)))==1, %y shrunk_ylimits don't include 0
1213         erp_yticks=shrunk_ylimits(1):ystep:shrunk_ylimits(2);
1214     else %y limits include 0
1215         % erp_yticks=0:ystep:shrunk_ylimits(2); %ensures 0 is a tick point
1216         % erp_yticks=[erp_yticks [0:-ystep:shrunk_ylimits(1)]];
1217         erp_yticks=0:ystep:axlimits(4); %ensures 0 is a tick point
1218         erp_yticks=[erp_yticks [0:-ystep:axlimits(3)]];
1219     end
1220     if strcmpi(erp_grid,'on')
1221         set(gca,'ytick',unique(erp_yticks),'ygrid','on');
1222     else
1223         set(gca,'ytick',unique(erp_yticks));
1224     end
1225 end
1226 
1227 
1228 %make ERP traces on top
1229 kids=get(gca,'children')';
1230 for hndl_loop=plot_handle,
1231    id=find(kids==hndl_loop);
1232    kids(id)=[];
1233    kids=[hndl_loop kids];
1234 end
1235 set(gca,'children',kids');
1236 plot_handle=[plot_handle l1 l2];
1237 
1238 
1239 %
1240 %% %%%%%%%%%%%%%%%%%%%%% 2. function nan_mean() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1241 %
1242 % nan_mean() - Take the column means of a matrix, ignoring NaN values.
1243 %              Return significance bounds if alpha (0 < alpha< <1) is given.
1244 %
1245 function [out, outalpha]  = nan_mean(in,alpha)
1246 NBOOT = 500;
1247 
1248 intrials = size(in,1);
1249 inframes = size(in,2);
1250 nans = find(isnan(in));
1251 in(nans) = 0;
1252 sums = sum(in);
1253 nonnans = ones(size(in));
1254 nonnans(nans) = 0;
1255 nonnans = sum(nonnans);
1256 nononnans = find(nonnans==0);
1257 nonnans(nononnans) = 1;
1258 out = sum(in)./nonnans;
1259 outalpha = [];
1260 if nargin>1
1261     if NBOOT < round(3/alpha)
1262         NBOOT = round(3/alpha);
1263     end
1264     fprintf('Computing %d bootstrap ERP values... ',NBOOT);
1265     booterps = zeros(NBOOT,inframes);
1266     for n=1:NBOOT
1267         signs = sign(randn(1,intrials)'-0.5);
1268         booterps(n,:) = sum(repmat(signs,1,inframes).*in)./nonnans;
1269         if ~rem(n,50)
1270             fprintf('%d ',n);
1271         end
1272     end
1273     fprintf('\n');
1274     booterps = sort(abs(booterps));
1275     alpha = 1+floor(2*alpha*NBOOT); % one-sided probability threshold
1276     outalpha = booterps(end+1-alpha,:);
1277 end
1278 out(nononnans) = NaN;
1279 
1280 %
1281 %% %%%%%%%%%%%%%%%%%%%%% 3. function compute_percentile() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1282 %
1283 function outpercent = compute_percentile(sortvar, percent, outtrials, winsize);
1284 ntrials = length(sortvar);
1285 outtrials=round(outtrials);
1286 sortvar = [ sortvar sortvar sortvar ];
1287 winvals = [round(-winsize/2):round(winsize/2)];
1288 outpercent = zeros(size(outtrials));
1289 for index = 1:length(outtrials)
1290     sortvarval = sortvar(outtrials(index)+ntrials+winvals);
1291     sortvarval = sort(sortvarval);
1292     outpercent(index) = sortvarval(round((length(winvals)-1)*percent)+1);
1293 end;
1294 
1295 %
1296 %% %%%%%%%%%%%%%%%%%%%%% 4. function orderofmag() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1297 %
1298 function ord=orderofmag(val)
1299 %function ord=orderofmag(val)
1300 %
1301 % Returns the order of magnitude of the value of 'val' in multiples of 10
1302 % (e.g., 10^-1, 10^0, 10^1, 10^2, etc ...)
1303 % used for computing ERPimage trial axis tick labels as an alternative for
1304 % plotting sorting variable
1305 
1306 if val>=1
1307     ord=1;
1308     val=floor(val/10);
1309     while val>=1,
1310         ord=ord*10;
1311         val=floor(val/10);
1312     end
1313     return;
1314 else
1315     ord=1/10;
1316     val=val*10;
1317     while val<1,
1318         ord=ord/10;
1319         val=val*10;
1320     end
1321     return;
1322 end
1323 
1324 %
1325 %% %%%%%%%%%%%%%%%%%%%%% 5. function find_crspnd_pt() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1326 %
1327 function y_pt=find_crspnd_pt(targ,vals,outtrials)
1328 %function id=find_crspnd_pt(targ,vals,outtrials)
1329 %
1330 % Inputs:
1331 %   targ      - desired value of sorting variable
1332 %   vals      - a vector of observed sorting variables (possibly smoothed)
1333 %   outtrials - a vector of y-axis values corresponding to trials in the
1334 %               ERPimage (this will just be 1:n_trials if there's no
1335 %               smoothing)
1336 %
1337 % Output:
1338 %   y_pt  - y-axis value (in units of trials) corresponding to "targ".
1339 %          If "targ" matches more than one y-axis pt, the median point is
1340 %          returned.  If "targ" falls between two points, y_pt is linearly
1341 %          interpolated.
1342 %
1343 % Note: targ and vals should be in the same units (e.g., milliseconds)
1344 
1345 %find closest point above
1346 abv=find(vals>=targ);
1347 if isempty(abv),
1348     %point lies outside of vals range, can't interpolate
1349     y_pt=[];
1350     return
1351 end
1352 abv=abv(1);
1353 
1354 %find closest point below
1355 blw=find(vals<=targ);
1356 if isempty(blw),
1357     %point lies outside of vals range, can't interpolate
1358     y_pt=[];
1359     return
1360 end
1361 blw=blw(end);
1362 
1363 if (vals(abv)==vals(blw)),
1364     %exact match
1365     ids=find(vals==targ);
1366     y_pt=median(outtrials(ids));
1367 else
1368     %interpolate point
1369     
1370     %lst squares inear regression
1371     B=regress([outtrials(abv) outtrials(blw)]',[ones(2,1) [vals(abv) vals(blw)]']);
1372     
1373     %predict outtrial point from target value
1374     y_pt=[1 targ]*B;
1375     
1376 end
1377 return
1378 %
1379 %% function 6. watchit
1380 %
1381 function watchit(msg)
1382 %function watchit(msg)
1383 %
1384 % Displays a warning message on the Matlab command line.
1385 %
1386 
1387 disp(' ');
1388 disp('****************** Warning ******************');
1389 disp(msg);
1390 disp(' ');
1391 return;