Peak detection in Python [Eli Billauer]

readme.md

I just directly translated to Python from this MATLAB script, which I've also included here. The original is public domain, so my translation is, too.

Also see Frequency estimation methods in Python for interpolating to get sharp intersample peaks

sixtenbe has posted a more powerful version here

and there's a PyPI repo

peakdet.m

function [maxtab, mintab]=peakdet(v, delta, x)
%PEAKDET Detect peaks in a vector
%        [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local
%        maxima and minima ("peaks") in the vector V.
%        MAXTAB and MINTAB consists of two columns. Column 1
%        contains indices in V, and column 2 the found values.
%      
%        With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices
%        in MAXTAB and MINTAB are replaced with the corresponding
%        X-values.
%
%        A point is considered a maximum peak if it has the maximal
%        value, and was preceded (to the left) by a value lower by
%        DELTA.

% Eli Billauer, 3.4.05 (Explicitly not copyrighted).
% This function is released to the public domain; Any use is allowed.

maxtab = [];
mintab = [];

v = v(:); % Just in case this wasn't a proper vector

if nargin < 3
  x = (1:length(v))';
else 
  x = x(:);
  if length(v)~= length(x)
    error('Input vectors v and x must have same length');
  end
end
  
if (length(delta(:)))>1
  error('Input argument DELTA must be a scalar');
end

if delta <= 0
  error('Input argument DELTA must be positive');
end

mn = Inf; mx = -Inf;
mnpos = NaN; mxpos = NaN;

lookformax = 1;

for i=1:length(v)
  this = v(i);
  if this > mx, mx = this; mxpos = x(i); end
  if this < mn, mn = this; mnpos = x(i); end
  
  if lookformax
    if this < mx-delta
      maxtab = [maxtab ; mxpos mx];
      mn = this; mnpos = x(i);
      lookformax = 0;
    end  
  else
    if this > mn+delta
      mintab = [mintab ; mnpos mn];
      mx = this; mxpos = x(i);
      lookformax = 1;
    end
  end
end

peakdetect.py

import sys
from numpy import NaN, Inf, arange, isscalar, asarray, array

def peakdet(v, delta, x = None):
    """
    Converted from MATLAB script at http://billauer.co.il/peakdet.html
    
    Returns two arrays
    
    function [maxtab, mintab]=peakdet(v, delta, x)
    %PEAKDET Detect peaks in a vector
    %        [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local
    %        maxima and minima ("peaks") in the vector V.
    %        MAXTAB and MINTAB consists of two columns. Column 1
    %        contains indices in V, and column 2 the found values.
    %      
    %        With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices
    %        in MAXTAB and MINTAB are replaced with the corresponding
    %        X-values.
    %
    %        A point is considered a maximum peak if it has the maximal
    %        value, and was preceded (to the left) by a value lower by
    %        DELTA.
    
    % Eli Billauer, 3.4.05 (Explicitly not copyrighted).
    % This function is released to the public domain; Any use is allowed.
    
    """
    maxtab = []
    mintab = []
       
    if x is None:
        x = arange(len(v))
    
    v = asarray(v)
    
    if len(v) != len(x):
        sys.exit('Input vectors v and x must have same length')
    
    if not isscalar(delta):
        sys.exit('Input argument delta must be a scalar')
    
    if delta <= 0:
        sys.exit('Input argument delta must be positive')
    
    mn, mx = Inf, -Inf
    mnpos, mxpos = NaN, NaN
    
    lookformax = True
    
    for i in arange(len(v)):
        this = v[i]
        if this > mx:
            mx = this
            mxpos = x[i]
        if this < mn:
            mn = this
            mnpos = x[i]
        
        if lookformax:
            if this < mx-delta:
                maxtab.append((mxpos, mx))
                mn = this
                mnpos = x[i]
                lookformax = False
        else:
            if this > mn+delta:
                mintab.append((mnpos, mn))
                mx = this
                mxpos = x[i]
                lookformax = True

    return array(maxtab), array(mintab)

if __name__=="__main__":
    from matplotlib.pyplot import plot, scatter, show
    series = [0,0,0,2,0,0,0,-2,0,0,0,2,0,0,0,-2,0]
    maxtab, mintab = peakdet(series,.3)
    plot(series)
    scatter(array(maxtab)[:,0], array(maxtab)[:,1], color='blue')
    scatter(array(mintab)[:,0], array(mintab)[:,1], color='red')
    show()

What I mean by "some amount of RMS noise tolerable'" is that the function also works when you have a noisy signal, I'm not adding any noisy if that's what you think. The sensitivity of the function can be set with the window variable. For a very noisy signal with maybe a thousand or a few thousand samples per period a window of about 20 should be enough, but I set it quite high to get a good margin and it doesn't effect the final result anyway, as long as it can find the zero-crossings properly. Should probably change the description to "repeatable signals, where some noise is tolerated" or something to that end.

As for performing interpolation I'm a little split about that, theoretically I think it might be a good practice, but experience hasn't always shown this. A colleague had a labView program for analysing waveforms, where he adapted a sinewave to every peak to find the actual peak smoothing away noise, but when we investigated the peaks we observed that it consistently choose a value lower than the actual peak and offset in time.

So the solution we use is to simply find the highest measured peak without trying to do anything smart. Although at another company they have fewer sample per period than we, but they perform an fft of the signal and then padds the fft-array with zero values and calculates the ifft of it to interpolate new values. This way they can have a longer aperture time for each sample increasing their accuracy for each sample and should theoretically retain resolution in time.

I would a define a realistic triangle as the following and by realistic I mean that you might actually encounter it on an actual power grid in a house or a factory:

triangle = lambda T: [1000*sqrt(2)*(sin(t)+0.05*sin(t*3-pi) + 0.05*sin(t*5)+0.02*sin(t*7-pi)+0.01*sin(t*9)) for t in T]
#or if you want a fixed frequency of the waveshape and not radians:
triangle = lambda T, Hz: [1000*sqrt(2)*(sin(Hz*pi*t)+0.05*sin(Hz*pi*t*3-pi) + 0.05*sin(Hz*pi*t*5)+0.02*sin(Hz*pi*t*7-pi)+0.01*sin(Hz*pi*t*9)) for t in T]

Yeah, that's basically what I was thinking with fitting a sine wave to the noisy sine wave and finding the peak of that instead. And yes, the parabolic interpolation is for sharp narrow peaks in a spectrum, not for finding peaks of sine waves.

Oh I notice that I managed to remove the context for my definition of a triangle. The point of the triangle is that a triangle and a sine wave, with some noise can be a good way of testing any function for fitting or interpolating a peak.

As for fitting sine waves, as I said I don't think it's worthwhile to fit any sine waves to the peak or interpolating it. For a pure sine it would also be good to compare a RMS calculation of the waveform with the Vpp/sqrt(8) where Vpp = difference between positive and negative peak. This is a good test to see if a function can find peaks for a pure sine wave. The triangle is useful when performing an optical inspection of the peak finding function.

For just finding the maximum respecitve the minimum value as done in my peakdetect_zero_crossing function I've verified that on a real pure sine wave it will give results of no worse than 60 ppm, with 1000 samples per period and calculated over 5 periods.

thanks for the script.

to make it work I had to add asarray and pylab imports:

from numpy import NaN, Inf, arange, isscalar, array, asarray
from pylab import *

as well as actually show the matplotlib graph by adding (in the main function):

plot(series)
show()

This piece of code saved Zeokat lots of hours. Thanks.

@v3nz3n sorry, I'm not good at listing imports because I just did from pylab import *, but I'm trying to be better now. Spyder has pyflakes built in, which helps.

Works perfectly. Thanks for sharing!!!

I have a negative signal with a peak pointing downwards. When I use this code maxtab gives me the start of the peak and mintab gives me the maximum position of the peak.
I would like to get the start and end of the peak. Is there a way to find?

Thanks for sharing! A great help.

Thanks you !

Thank you.

Thank you! Great help for search target waves in EEG data.

Thanks a lot for this, do you have this in R?

thanks for the code!

please how to use this with PyOpenCl ?

Here is a go version I converted https://github.com/mpiannucci/peakdetect

Is there a way I can preprocess csv files and show output for accelerometer?

can i implement it in GNURadio? Thank you very much

Mate, I just want to say Thank you so much because you save my life with this script.

Great coding ! Thanks a lot!

Thank you so much

有关峰值点检测的研究促进这一领域的快速发展，我希望有志之士团结合作共同推进，同时我正在做峰值检测的学术课题，希望与有意者交流探讨，联系方式：m15738897318@gail.com,谢谢

Thank you for this snippet, it works perfectly in my implementation.

Thanks man!

Thanks for sharing your code! It works very well for my project.