public
Last active

Peak detection in Python

  • Download Gist
peakdet.m
Matlab
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
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
Python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
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()
readme.md
Markdown

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

Justin

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.

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.

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

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.

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.