Detecting peaks with MatLab
For those not familiar to digital signal processing, peak detection is as easy to understand as it sounds: this is the process of finding peaks - we also names them local maxima or local minima - in a signal. The MatLab DSP Toolbox makes this super easy with its
findpeaks function. Saying your want to search local maxima in an audio signal, for example 2000 samples of the Laurent Garnier famous track Cripsy Bacon, all you have to do is:
We can specify filtering options to the function so the peaks that do not interest us are discarded. All this is great, but we need something working in Python. In a perfect world it will give exactly the same output, so we have consistent results between our Python code and the MatLab code.
The Scipy try
Contrary to other MatLab functions that have direct equivalents in the Numpy and Scipy scientific and processing packages, it is no easy task to get the same results from the Scipy
find_peaks_cwt function that from the MatLab
findpeaks. Even worse, the wavelet convolution approach of
find_peaks_cwt is not straightforward to work with: it adds complexity that is of no use for well-filtered and noiseless signals.
Into the wild
Wondering how to make our algorithms works as simply with Python that they were in MatLab, I’ve search around the web for other peak detection algorithms available in Python. Stackoverflow get me to
peakdetect, a translation of a MatLab script. As it is clearly more trivial to use that
find_peaks_cwt, it still won’t give you the same results that the MatLab
findpeaks function. The algorithm don’t find all peaks on low sampled signals or on short samples, and don’t have either a support for minimum peak height filter.
Going ahead I’ve checkout the GNU Octave project, a processing-intended language quite similar to MatLab. The Octave-Forge repository hosts a digital signal processing package with a
findpeaks function. Bingo? Well, yes and no. The function have an appealing interface, with a great filtering support. But this is not a Python project: as you’ll find ways to call your Octave distribution from your Python code (see oct2py), it surely won’t be effective at large scale and makes the requirements for your code more complex.
The PyPi Wonderland
As I was going to code a Python adaptation of the Octave-Force
findpeaks, I finally found what I was searching: a Python native equivalent of the MatLab
findpeaks, with minimum distance and height filtering support. I even found two!
Going through space
The second is published on a jupyter notebook and is written by Marcos Duarte. Not easy to find, but the greatest so far for what I need: the parameters for filtering are modeled after the MatLab
findpeaks ones. It comes as a single source file and only depends on Numpy, so it is no big deal to integrate. And more importantly, it will consistently get you the same results than MalLab
To avoid others the same roaming I’ve put on GitHub an overview of these findings. Let me know if you got another open-source alternatives so we update the list.
Edit 17th November
As Lucas Hermann Negri pointed out on HN, the PeakUtils
indexes function was actually inspired by Marcos Duarte implementation of
detect_peaks, explaining the similar results. Moreover he notes that the PeakUtils comes with other convenient utilities, such as
interpolate function enhances the peak resolution by fitting Gaussians or computing centroids. Taking the previous example, here how you get it work:
This time before the peak resolution, the
baseline function will be very handy in presence of drifting signals or to deal with unwanted low-frequency phenomenon: it kind of high-pass filter the signal. The PeakUtils documentation have a good example of its use.