Peak detection in Python

peakdetect.py

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

def peakdet(v, delta, x = None):
    """
    Converted from MATLAB script at http://billauer.co.il/peakdet.html
    
    Currently returns two lists of tuples, but maybe arrays would be better
    
    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 maxtab, mintab

if __name__=="__main__":
    series = [0,0,0,2,0,0,0,-2,0,0,0,2,0,0,0,-2,0]
    print peakdet(series,1)

Thanks that is pretty useful!

Found an embarrassing error in the peakdetect_zero_crossing function, where it flipped the x and y values in the return. Also added some functions for fitting model functions to the waveform and a fft based interpolation, although at this time I only recommend the peakdetect and peakdetect_zero_crossing function. The reason for this is that I don't yet know if the other functions actually increases the accuracy of the peak detection. Being a metrologist I want to know the uncertainty of my peak detection funtion. For raw peak detection function this can easily be calculated as a function of the time resolution and signal frequency.

Oh the new fitting functions look good. I migth give those a try. I found a Savitzky-Golay filter pretty helpfull for my data. (Analytical chemistry)
Thanks for the update

Correct me if I'm wrong, but isn't the if index+lookahead >= length unnecessary in the peakdetect function, considering that the loop is from [0,length-lookahead) ?

This is very useful! Thanks!

great piece of code

great code, just a small issue correct me if I'm wrong, should
line 140 be mn = y
line 158 be mx = y

thanks

very useful, nice!

Thanks for the code.

What is the license associated with this code?

I have a negative signal that has a peak, but this code lists all the values of x and y

Hi Sixtenbe! I implemented the peakdetect.py into the algorithm of my master thesis and it works like a charm:). Therefore I wanted to ask you what the license for this code is and how you want me to reference you? Thanks!

great tool. i am using it for an acoustic experiment designed to behave like a quantum state. lots of data sets with lots of peaks; lots of manual work done by your script.
Thanks

highly appreciate your work.

Great work, @sixtenbe! Would you mind adding a license to the script? I would like to include it in a BSD-licensed library for metabolomics, https://github.com/metabolite-atlas/metatlas.

This works really well for identifying peaks in a saw-tooth pattern from atomic force microscopy. I would also want to ask for the license.

Hi, @sixtenbe, thank you for coding this. Do you release this under the MIT license?

Killer code. Thanks a bunch for sharing!

Just saw that I didn't actually get any kind of notifications for this code and lots of people are asking for the license.

The original piece of code that this was based on is in public domain see http://billauer.co.il/peakdet.html. The contributions I've made to the code can be considered licensed under http://www.wtfpl.net/, which means that you can do whatever you want with the code.

I might do some more with this code as I want to implement some automatic tests of it to verify that everything is returning sane values for a few analytic waveforms that will be based on an IEC standard for high voltage measurement systems.

Refactored a bit of the code and added some automatic testing of the peakdetection functions so that any new function can be verified to function to some standard of bare minimum.

I also changed the name of peakdetect_parabole to peakdetect_parabola, but added a legacy function so that the old name will still work.
I changed the name of the window variable in zero_crossings to window_len and added the variable window_f to enable specification of window type.

Future improvements that I really should do are to make all the function handle offsets without breaking and possibly add some further tests for triangle waveshapes ac defined in ACV_5 or ACV_6, but for this I would need to decide how I define what the maxima and minima are.

please how to use this with PyOpenCl ?

@algerianmaster I'm not sure what you're asking as to use something with PyOpenCl you would also have to be interested in making calls to the OpenCl API and why would even want from that. To use this with PyOpenCl you would have to rewrite all the appropriate math to use OpenCl calls and once again: Why?

Very useful code, and I have been quite happy with the peakdetect function itself. But after trying the peakdetect_parabola, I am wondering if the zero_crossing function has some prerequisites for the data? In particular, my data has a DC offset. So the line

indices = np.where(np.diff(np.sign(y_axis)))[0]
# check if zero-crossings are valid
diff = np.diff(indices)

results in indices and diff being empty arrays. This fails the test in the next line (nan) and finally comes out with the ValueError:

    if len(indices) < 1:
        raise ValueError("No zero crossings found")

So it seems the offset correction in the middle is ineffective.

EDIT: Just noticed your note that this function is ineffective for large offset. So this is not a bug.

@subhacom A quick fix you can do on your end is to simply do a offset = np.mean(np.max(my_data), np.min(my_data))
and subtract that from your data and you should end up with zero crossing that are close enough to the real ones for all peak detect methods to find their peaks.

As you probably found from the code it currently only takes DC offsets less than half the amplitude as it must have a zero crossing or it assumes it's non-periodic data.

And you can also use the peakdetect function, which only finds the raw values and will work with any offsets. It can take a constant DC offset and even a linearly increasing offset as it's searching for local extremes. The downside is that it's a bit slower than the other functions due to how it crawls over all data points.

Any chance to get this published to pypi? Would cool to make this installable!

published: https://pypi.python.org/pypi/analytic-wfm

Thanks. The code saved me a lot of time. @cancan101 Is the pypi version being updated regularly? Not that I'm going to use it, but I might have to suggest the pypi version to my friends.

@Sakthi-G33k it not like this gist is seeing much in terms of updates, so I don't see any pypi getting regular updates either. There isn't much to update, unless adding additional logic like processing the indata for frequency domain padding for time domain interpolation to reduce edge effects, but personally I would prefer to just sample 2**n samples over x amount of whole periods of the signal.

Revised and published on PyPI and GitHub: https://pypi.org/project/peakdetect/

Hi,
The code is not working under python 3.8.6.
Below I put the diff for changes I made to make it work.

diff --git a/peakdetect.py b/peakdetect.py
index e388408..5d070fa 100644
--- a/peakdetect.py
+++ b/peakdetect.py
@@ -333,7 +333,7 @@ def peakdetect_fft(y_axis, x_axis, pad_len = 20):
     #max_peaks, min_peaks = peakdetect_zero_crossing(y_axis_ifft, x_axis_ifft)
     
     # store one 20th of a period as waveform data
-    data_len = int(np.diff(zero_indices).mean()) / 10
+    data_len = int(np.diff(zero_indices).mean()) // 10
     data_len += 1 - data_len & 1
diff --git a/test.py b/test.py
index 680c2a8..c7c6dfa 10644
--- a/test.py
+++ b/test.py
@@ -31,7 +31,7 @@ def prng():
 
 
 def _write_log(file, header, message):
-    with open(file, "ab") as f:
+    with open(file, "a") as f:
         f.write(header)
         f.write("\n")
         f.writelines(message)
@@ -241,6 +241,8 @@ class _Test_peakdetect_template(unittest.TestCase):
             *self.args, **self.kwargs
             )
         #check if the correct amount of peaks were discovered
+        max_p = list(max_p)
+        min_p = list(min_p)
         self.assertIn(len(max_p), [4,5])
         self.assertIn(len(min_p), [4,5])