Instantly share code, notes, and snippets.

Embed
What would you like to do?
Peak detection in Python [Eli Billauer]
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
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()
@wyojustin

This comment has been minimized.

wyojustin commented Jul 9, 2011

Thanks for sharing your script. This saves me a huge headache!

Justin

@dusita

This comment has been minimized.

dusita commented Aug 10, 2011

Thanks a lot for sharing this script. It works well for me too.
This script causes me register to github as I would like to say thank you.

@endolith

This comment has been minimized.

Owner

endolith commented Aug 10, 2011

@sixtenbe

This comment has been minimized.

sixtenbe commented Aug 30, 2011

Great script. Gave me a few good ideas about how to find peaks, but I think you should note that it if your signal starts with a minima you will systematically miss that one. Just remove the need for the: 'If lookoformax' statement and replace with some other logic and you can find the first peak, even if it's a minima. Doing this also means that you almost always will have a false hit on the first or one of the first samples and will need to pop that peak from the results.

Changed this logic, added another test criteria to find the peaks and I also made a peak finder routine which divided the signal into different sections by looking for zero crossings and finding the peak in each section.

@endolith

This comment has been minimized.

Owner

endolith commented Sep 9, 2011

That's great. I'll probably use that version in the future.

For the "repeatable sinusoidal signals with some amount of RMS noise tolerable" version, why would you want the peak of the signal+noise? I would think you would want the peak of the signal itself?

Also see https://gist.github.com/255291#file_parabolic.py

@sixtenbe

This comment has been minimized.

sixtenbe commented Sep 9, 2011

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]
@endolith

This comment has been minimized.

Owner

endolith commented Sep 9, 2011

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.

@sixtenbe

This comment has been minimized.

sixtenbe commented Sep 9, 2011

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.

@v3nz3n

This comment has been minimized.

v3nz3n commented Oct 25, 2013

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

This comment has been minimized.

Zeokat commented Mar 5, 2014

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

@endolith

This comment has been minimized.

Owner

endolith commented Mar 18, 2014

@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.

@Alymantara

This comment has been minimized.

Alymantara commented Aug 27, 2014

Works perfectly. Thanks for sharing!!!

@San08

This comment has been minimized.

San08 commented Sep 2, 2014

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?

@iuribeferrari

This comment has been minimized.

iuribeferrari commented Jan 15, 2015

Thanks for sharing! A great help.

@MartinSamouiller

This comment has been minimized.

MartinSamouiller commented Sep 1, 2015

Thanks you !

@bancek

This comment has been minimized.

bancek commented Oct 11, 2015

Thank you.

@adowaconan

This comment has been minimized.

adowaconan commented Mar 21, 2016

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

@Ruthygg

This comment has been minimized.

Ruthygg commented Mar 22, 2016

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

@golobor

This comment has been minimized.

golobor commented Jul 1, 2016

thanks for the code!

@algerianmaster

This comment has been minimized.

algerianmaster commented Aug 3, 2016

please how to use this with PyOpenCl ?

@mpiannucci

This comment has been minimized.

mpiannucci commented Sep 20, 2016

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

@Rishie13

This comment has been minimized.

Rishie13 commented Dec 3, 2016

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

@atikalidyad

This comment has been minimized.

atikalidyad commented May 23, 2017

can i implement it in GNURadio? Thank you very much

@JorgeHormaza

This comment has been minimized.

JorgeHormaza commented Dec 4, 2017

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

@PZiso

This comment has been minimized.

PZiso commented Mar 29, 2018

Great coding ! Thanks a lot!

@LuisAli22

This comment has been minimized.

LuisAli22 commented Jun 20, 2018

Thank you so much

@15738897318

This comment has been minimized.

15738897318 commented Oct 7, 2018

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

@primehaxor

This comment has been minimized.

primehaxor commented Nov 17, 2018

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment