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