function hL = dsplot(x, y, numPoints)
  
  %DSPLOT Create down sampled plot.
  %   This function creates a down sampled plot to improve the speed of
  %   exploration (zoom, pan).
  %
  %   DSPLOT(X, Y) plots Y versus X by downsampling if there are large number
  %   of elements. X and Y needs to obey the following:
  %     1. X must be a monotonically increasing vector.
  %     2. If Y is a vector, it must be the same size as X.
  %     3. If Y is a matrix, one of the dimensions must line up with X.
  %
  %   DSPLOT(Y) plots the columns of Y versus their index.
  %
  %   hLine = DSPLOT(X, Y) returns the handles of the line. Note that the
  %   lines may be downsampled, so they may not represent the full data set.
  %
  %   DSPLOT(X, Y, NUMPOINTS) or DSPLOT(Y, [], NUMPOINTS) specifies the
  %   number of points (roughly) to display on the screen. The default is
  %   50000 points (~390 kB doubles). NUMPOINTS can be a number greater than
  %   500.
  %
  %   It is very likely that more points will be displayed than specified by
  %   NUMPOINTS, because it will try to plot any outlier points in the range.
  %   If the signal is stochastic or has a lot of sharp changes, there will
  %   be more points on plotted on the screen.
  %
  %   The figure title (name) will indicate whether the plot shown is
  %   downsampled or is the true representation.
  %
  %   The figure can be saved as a .fig file, which will include the actual
  %   data. The figure can be reloaded and the actual data can be exported to
  %   the base workspace via a menu.
  %
  %   Run the following examples and zoom/pan to see the performance.
  %
  %  Example 1: (with small details)
  %   x  = linspace(0, 2*pi, 1000000);
  %   y1 = sin(x)+.02*cos(200*x)+0.001*sin(2000*x)+0.0001*cos(20000*x);
  %   dsplot(x,y1);title('Down Sampled');
  %   % compare with
  %   figure;plot(x,y1);title('Normal Plot');
  %
  %  Example 2: (with outlier points)
  %   x  = linspace(0, 2*pi, 1000000);
  %   y1 = sin(x) + .01*cos(200*x) + 0.001*sin(2000*x);
  %   y2 = sin(x) + 0.3*cos(3*x)   + 0.001*randn(size(x));
  %   y1([300000, 700000, 700001, 900000]) = [0, 1, -2, 0.5];
  %   y2(300000:500000) = y2(300000:500000) + 1;
  %   y2(500001:600000) = y2(500001:600000) - 1;
  %   y2(800000) = 0;
  %   dsplot(x, [y1;y2]);title('Down Sampled');
  %   % compare with
  %   figure;plot(x, [y1;y2]);title('Normal Plot');
  %
  %  See also PLOT.
  
  %  Version:
  %   v1.0 - first version (Aug 1, 2007)
  %   v1.1 - added CreateFcn for the figure so that when the figure is saved
  %          and re-loaded, the zooming and panning works. Also added a menu
  %          item for saving out the original data back to the base
  %          workspace. (Aug 10, 2007)
  %
  %  Jiro Doke
  %  August 1, 2007
  
  debugMode = false;
  
  %--------------------------------------------------------------------------
  % Error checking
  error(nargchk(1, 3, nargin, 'struct'));
  if nargin < 3
    % Number of points to show on the screen. It's quite possible that more
    % points will be displayed if there are outlier points
    numPoints = 50000;  % ~390 kB for doubles
  end
  if nargin == 1 || isempty(y)
    noXVar = true;
    y = x;
    x = [];
  else
    noXVar = false;
  end
  myErrorCheck;
  %--------------------------------------------------------------------------
  
  if size(x, 2) > 1  % it's a row vector -> transpose
    x = x';
    y = y';
    varTranspose = true;
  else
    varTranspose = false;
  end
  
  % Number of lines
  numSignals = size(y, 2);
  
  % If the number of lines is greater than the number of data points per
  % line, it's possible that the user may have mistaken the matrix
  % orientation.
  if numSignals > size(y, 1)
    s = input(sprintf('Are you sure you want to plot %d lines? (y/n) ', ...
      numSignals), 's');
    if ~strcmpi(s, 'y')
      disp('Canceled. You may want to transpose the matrix.');
      if nargout == 1
        hL = [];
      end
      return;
    end
  end
  
  % Attempt to find outliers. Use a running average technique
  filterWidth = ceil(min([50, length(x)/10])); % max window size of 50
  a  = y - filter(ones(filterWidth,1)/filterWidth, 1, y);
  [iOutliers, jOutliers] = find(abs(a - repmat(mean(a), size(a, 1), 1)) > ...
    repmat(4 * std(a), size(a, 1), 1));
  clear a;
  
  % Always create new figure because it messes around with zoom, pan,
  % datacursors.
  hFig    = figure;
  figName = '';
  
  % Create template plot using NaNs
  hLine   = plot(NaN(2, numSignals), NaN(2, numSignals));
  set(hLine, 'tag', 'dsplot_lines');
  
  % Define CreateFcn for the figure
  set(hFig, 'CreateFcn', @mycreatefcn);
  mycreatefcn();
  
  % Create menu for exporting data
  hMenu = uimenu(hFig, 'Label', 'Data');
  uimenu(hMenu, ...
    'Label'   , 'Export data to workspace.', ...
    'Callback', @myExportFcn);
  
  % Update lines
  updateLines([min(x), max(x)]);
  
  % Deal with output argument
  if nargout == 1
    hL = hLine;
  end
  
  %--------------------------------------------------------------------------
    function myExportFcn(varargin)
      % This callback allows for extracting the actual data from the figure.
      % This means that if you save this figure and load it back later, you
      % can get back the data.
      
      % Determine the variable name
      allVarNames = evalin('base', 'who');
      newVarName = genvarname('dsplotData', allVarNames);
      
      % X
      if ~noXVar
        if varTranspose
          dat.x = x';
        else
          dat.x = x;
        end
      end
      
      % Y
      if varTranspose
        dat.y = y';
      else
        dat.y = y;
      end
      
      assignin('base', newVarName, dat);
      
      msgbox(sprintf('Data saved to the base workspace as ''%s''.', ...
        newVarName), 'Saved', 'modal');
      
    end
  
  %--------------------------------------------------------------------------
    function mycreatefcn(varargin)
      % This callback defines the custom zoom/pan functions. It is defined as
      % the CreateFcn of the figure, so it allows for saving and reloading of
      % the figure.
  
      if nargin > 0
        hFig = varargin{1};
      end
      hLine = findobj(hFig, 'type', 'axes');
      hLine(strmatch('legend', get(hLine, 'tag'))) = [];
      hLine = get(hLine, 'Children');
      
      % Create Zoom, Pan, Datacursor objects
      hZoom = zoom(hFig);
      hPan  = pan(hFig);
      hDc   = datacursormode(hFig);
      set(hZoom, 'ActionPostCallback', @mypostcallback);
      set(hPan , 'ActionPostCallback', @mypostcallback);
      set(hDc  , 'UpdateFcn'         , @myDCupdatefcn);
  
    end
  
  %--------------------------------------------------------------------------
    function mypostcallback(obj, evd) %#ok
      % This callback that gets called when the mouse is released after
      % zooming or panning.
  
      % single or double-click
      switch get(hFig, 'SelectionType')
        case {'normal', 'alt'}
          updateLines(xlim(evd.Axes));
  
        case 'open'
          updateLines([min(x), max(x)]);
  
      end
  
    end
  
  %--------------------------------------------------------------------------
    function updateLines(rng)
      % This helper function is for determining the points to plot on the
      % screen based on which portion is visible in the current limits.
  
      % find indeces inside the range
      id = find(x >= rng(1) & x <= rng(2));
  
      % if there are more points than we want
      if length(id) > numPoints / numSignals
  
        % see how many outlier points are in this range
        blah = iOutliers > id(1) & iOutliers < id(end);
  
        % determine indeces of points to plot. 
        idid = round(linspace(id(1), id(end), round(numPoints/numSignals)))';
  
        x2 = cell(numSignals, 1);
        y2 = x2;
        for iSignals = 1:numSignals
          % add outlier points
          ididid = unique([idid; iOutliers(blah & jOutliers == iSignals)]);
          x2{iSignals} = x(ididid);
          y2{iSignals} = y(ididid, iSignals);
        end
  
        if debugMode
          figName = ['downsampled - ', sprintf('%d, ', cellfun('length', y2))];
        else
          figName = 'downsampled';
        end
  
      else % no need to down sample
        figName = 'true';
  
        x2 = repmat({x(id)}, numSignals, 1);
        y2 = mat2cell(y(id, :), length(id), ones(1, numSignals))';
  
      end
  
      % Update plot
      set(hLine, {'xdata', 'ydata'} , [x2, y2]);
      set(hFig, 'Name', figName);
  
    end
  
  %--------------------------------------------------------------------------
    function txt = myDCupdatefcn(empt, event_obj) %#ok
      % This function displays appropriate data cursor message based on the
      % display type
  
      pos = get(event_obj,'Position');
      switch figName
        case 'true'
          txt = {['X: ',num2str(pos(1))],...
            ['Y: ',num2str(pos(2))]};
        otherwise
          txt = {['X: ',num2str(pos(1))],...
            ['Y: ',num2str(pos(2))], ...
            'Warning: Downsampled', ...
            'May not be accurate'};
      end
    end
  
  %--------------------------------------------------------------------------
    function myErrorCheck
      % Do some error checking on the input arguments.
  
      if ~isa(numPoints, 'double') || numel(numPoints) > 1 || numPoints < 500
        error('Third argument must be a scalar greater than 500');
      end
      if ~isnumeric(x) || ~isnumeric(y)
        error('Arguments must be numeric');
      end
      if length(size(x)) > 2 || length(size(y)) > 2
        error('Only 2-D data accepted');
      end
      
      % If only one input, create index vector X
      if isempty(x)
        if ismember(1, size(y))
          x = reshape(1:numel(y), size(y));
        else
          x = (1:size(y, 1))';
        end
      end
      
      if ~ismember(1, size(x))
        error('First argument has to be a vector');
      end
      if ~isequal(size(x, 1), size(y, 1)) && ~isequal(size(x, 2), size(y, 2))
        error('One of the dimensions of the two arguments must match');
      end
      if any(diff(x) <= 0)
        error('The first argument has to be a monotonically increasing vector');
      end
    end
  
  end