function P=peaksat(x,y,threshold)
% P=peaksat(x,y,threshold)
% peaksat.m (Peaks Above Threshold) lists every y value that (a) has lower
% y values on both sides and (b) is above the specified threshold. Returns
% a 2 x n matrix P with the x and y values of each peak, where n is the
% number of peaks detected. (Noisy data must be smoothed first).
% Example:
%  x=2:38
%  y=abs(randn(size(x)));
%  P=peaksat(x,y,1);
%  disp('           peak      position      height    ')
%  disp(P)
%  plot(x,y)
%  text(P(:,2),P(:,3),num2str(P(:,1))) % Numbers peaks
%
peak=0;
for k = 2:length(x)-1
    if y(k) > y(k-1)
        if y(k) > y (k+1)
            if y(k)>threshold
                peak = peak + 1;
                P(peak,:)=[peak x(k) y(k)];
            end
        end
    end
end