1D Non-derivative Peak Finder

Up-sample and filter noisy data to find peaks without using derivatives.

PEAKFIND general 1D peak finding algorithm
Tristan Ursell, 2013.

peakfind(x_data,y_data)
peakfind(x_data,y_data,upsam)
peakfind(x_data,y_data,upsam,gsize,gstd)
peakfind(x_data,y_data,upsam,htcut,'cuttype')
peakfind(x_data,y_data,upsam,gsize,gstd,htcut,'cuttype')

[xpeaks]=peakfind()
[xout,yout,peakspos]=peakfind()

This function finds peaks without taking first or second derivatives, rather it uses local slope features in a given data set. The function has four basic modes.

Mode 1: peakfind(x_data,y_data) simply finds all peaks in the data
given by 'xdata' and 'ydata'.

Mode 2: peakfind(x_data,y_data,upsam) finds peaks after up-sampling
the data by the integer factor 'upsam' -- this allows for higher
resolution peak finding. The interpolation uses a cubic spline that
does not introduce fictitious peaks.

Mode 3: peakfind(x_data,y_data,upsam,gsize,gstd) up-samples and then
convolves the data with a Gaussian point spread vector of length
gsize (>=3) and standard deviation gstd (>0).

Mode 4: peakfind(x_data,y_data,upsam,htcut,'cuttype') up-samples the
data (upsam=1 analyzes the data unmodified). The string 'cuttype' can
either be 'abs' (absolute) or 'rel' (relative), which specifies a peak
height cutoff which is either:

'abs' - finds peaks that lie an absolute amount 'htcut' above the
lowest value in the data set.

for (htcut > 0) peaks are found if

peakheights > min(yout) + htcut

'rel' - finds peaks that are an amount 'htcut' above the lowest
value in the data set, relative to the full change in y-input
values.

for (0 < htcut < 1) peaks are found if

(peakheights-min(yout))/(max(yout)-min(yout)) > htcut

Upsampling and convolution allows one to find significant peaks in noisy data with sub-pixel resolution. The algorithm also finds peaks in data where the peak is surrounded by zero first derivatives, i.e. the peak is actually a large plateau.

The function outputs the x-position of the peaks in 'xpeaks' or the processed input data in 'xout' and 'yout' with 'peakspos' as the indices of the peaks, i.e. xpeaks = xout(peakspos).

If you want the algorithm to find the position of minima, simply input '-y_data'. Peaks within half the convolution box size of the boundary will be ignored (to avoid this, pad the data before processing).

Example 1:

x_data = -50:50;
y_data =(sin(x_data)+0.000001)./(x_data+0.000001)+1+0.025*(2*rand(1,length(x_data))-1);

[xout,yout,peakspos]=peakfind(x_data,y_data);
plot(x_data,y_data,'r','linewidth',2)
hold on
plot(xout,yout,'b','linewidth',2)
plot(xout(peakspos),yout(peakspos),'g.','Markersize',30)
xlabel('x')
ylabel('y')
title(['Found ' num2str(length(peakspos)) ' peaks.'])
box on

Example 2:

x_data = -50:50;
y_data =(sin(x_data)+0.000001)./(x_data+0.000001)+1+0.025*(2*rand(1,length(x_data))-1);

[xout,yout,peakspos]=peakfind(x_data,y_data,4,6,2,0.2,'rel');

plot(x_data,y_data,'r','linewidth',2)
hold on
plot(xout,yout,'b','linewidth',2)
plot(xout(peakspos),yout(peakspos),'g.','Markersize',30)
xlabel('x')
ylabel('y')
title(['Found ' num2str(length(peakspos)) ' peaks.'])
box on