% FINDPEAK  ( MatLinks) Find the peaks (and valleys) of a function.
%
%    FINDPEAK(X) finds the peaks in the data vector X.
%
%    FINDPEAK(X,S) ignores S "peaks" from either end of X.  This is necessary
%    when the first and last points of X are not peaks themselves, since
%    ordinarily FINDPEAK(...) considers the endpoints of X to be peaks.  The
%    default value for S=1, which for noiseless oscillations X is sufficient to
%    discard the endpoints as being peaks.  Use S=0 to retain the endpoints as
%    being peaks themselves, or some value S>1 to remove false peaks from
%    oscillations with noisy endpoints.
%
%    FINDPEAK(...) by itself plots X and marks the peaks and valleys with small
%    crosshairs.
%
%    I=FINDPEAK(...) returns the indices I closest to the valleys and peaks.
%    The first column contains the valley indices, the second column the peak
%    indices.  LENGTH(I) will thus give the total number of peaks and valleys.
%
%    [I,J]=FINDPEAK(...) returns the indices I closest to each peak minimum, and
%    the indices J closest to each peak maximum.  LENGTH(I)+LENGTH(J) will thus
%    equal the total number of peaks and valleys.
%
%    If no peaks are detected then I=J=0 is returned.
%
%    See also SETTLING, STEADY, ZCROSS, PHDIFF.
%
%    Type HELP MATLINKS for a full listing of all MatLinks ToolChest functions.
%
function [minima, maxima] = findpeak(x, strip)
%===============================================================================
%  Copyright  1998,2000 Julian Andrew de Marchi, Ph.D. (julian@matlinks.net)
%  Use & distribution covered by GNU General Public License (www.gnu.org)
%===============================================================================

%------------------
% parse the inputs
%------------------
if (nargin == 0),  error('No data vector X supplied.');
elseif (nargin == 1),  strip = 1;
end;


%----------------------------
% find the peak_inx and maxima
%----------------------------
peak_inx = 1;  inx = 2;  n = 1;
while (n < length(x)),

  % find the next positive inflection
  m = n;
  while (n < length(x) & x(n) >= x(n+1)),
    n = n + 1;
  end;
  if (m ~= n & x(m) ~= x(n)),
    peak_inx(inx) = n;  n = n + 1;  inx = inx + 1;
  end;

  % find the next negative inflection
  m = n;
  while (n < length(x) & x(n) <= x(n+1)),
    n = n + 1;
  end;
  if (m ~= n & x(m) ~= x(n)),
    peak_inx(inx) = n;  n = n + 1;  inx = inx + 1;
  end;

end;
if (length(peak_inx) == 1),
  warning('No peaks detected.');
  return;
end;


%----------------------
% strip boundary peaks
%----------------------
if (strip),
  strip = min([strip fix(length(peak_inx)/2)]);
  peak_inx = peak_inx(1+strip : length(peak_inx)-strip);
end;


%-------------------
% parse the outputs
%-------------------
if (nargout ~= 1),                    % extract maxima from peak_inx           %

  if (x(peak_inx(1)) < x(peak_inx(2))),
    maxima = peak_inx(2:2:length(peak_inx));
    minima = peak_inx(1:2:length(peak_inx));
  else
    maxima = peak_inx(1:2:length(peak_inx));
    minima = peak_inx(2:2:length(peak_inx));
  end;

end;

if (nargout == 0),                    % plot the peaks                         %

  hold off, plot(x), hold on, plot(1:length(x), x, 'c.');
  plot(minima, x(minima), 'm+', minima, x(minima), 'y.');
  plot(maxima, x(maxima), 'g+', maxima, x(maxima), 'y.');
  title('Peaks'),
  xlabel('i'), ylabel('x(i)'), hold off, zoom on;

end;

if (nargout <= 1),                    % set combination return value           %
  minima = peak_inx;
end;


%===============================================================================
% End-of-File
%===============================================================================