function [tc,yc] = wave_compress(t,y)
% Compress time trace by removing points that are on a straight line
% Replicates part of WAVE_PUT (TDI function) behaviour
assert(iscolumn(t),        'WAVECOMPRESS:t_size',         't must be a row vector')
assert(size(y,1)==numel(t),'WAVECOMPRESS:invalid_y_size', 'y must have as many columns as time points');
assert(ismatrix(y),        'WAVECOMPRESS:y_matrix',       'y must be a vector or matrix')
assert(all(sort(t)==t) && (numel(unique(t))==numel(t)),...
                           'WAVECOMPRESS:t_sorted','t must be sorted and unique')
if numel(t)<3 % trivial case
  yc = y;
  return;
end

yc = zeros(size(y)); tc=zeros(size(t));

tol = 1e-8; % tolerance

ii=1; % current point
jj=2; % test point
kk=1; % index into output trace

% first point
yc(1,:) = y(1,:);
tc(1)   = t(1);

while (jj+1)<=numel(t)
  % if this point jj is on a straight line w.r.t ii and jj+1
  % linear fit parameters
  cross = (t(jj+1)-t(ii))*(y(jj,:)-y(ii,:)) - (t(jj)-t(ii))*(y(jj+1,:)-y(ii,:));
  norm1 = sqrt((t(jj+1)-t(ii)).^2 + (y(jj+1,:)-y(ii,:)).^2);
  norm2 = sqrt((t(jj  )-t(ii)).^2 + (y(jj  ,:)-y(ii,:)).^2);
  if any(abs(cross)./norm1./norm2>tol)
    % store this point, use as starting point for next search
    kk = kk+1;
    yc(kk,:) = y(jj,:);
    tc(kk)   = t(jj,:);
    ii = jj;
    jj = ii+1;
  else
    jj = jj+1; % next point to check for linearity
  end
end
% store last point
% last point, store and exit
kk=kk+1;
yc(kk,:)= y(jj,:);
tc(kk)  = t(jj);

% crop
yc = yc(1:kk,:);
tc = tc(1:kk,:);

return
