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()

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

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

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.