# Wavelet-Based Peak Detection

GROUPS:
 Christopher Wolfram 6 Votes In Mathematica the closest thing (currently) to a peak detect function is a function called MaxDetect, which is unfortunately very slow with large datasets and could be better at finding peaks.  So (for a project which I will post here soon) I decided to write my own peak detection function.I found a nice article from National Instruments here: http://www.ni.com/white-paper/5432/en/I have subsequently implemented this in Mathematica.So the basic idea behind this algorithm is whenever the detail wavelet passes zero, there is a peak.  Let me explain this:(I am assuming that you understand the basic concept of wavelets)When you get the wavelet of a 1-dimensional object, you will get the approximation coefficients, and the detail coefficients.  As their names suggest the approximation coefficients are the basic shape of the thing (used for noise reduction) while the detail coefficients are used for the little bumps and valleys that texture the original wave.For exampleHere is some noisy data:Approximate:Detail:So to access the detail vs approx wavelets from some wavelet data like this:data = Table[Sin[x] + RandomReal[{0, .1}], {x, 0, 2 Pi, .01}];dwd = DiscreteWaveletTransform[data, HaarWavelet[], 3]So I am running a wavelet transform of data with a haar wavelet for 3 layers.I can see the various coefficients like this:dwd["ListPlot"]The coefficients that end in 1 are detail wavelets while those which end in 0 are approx wavelets.(There is a {0} and a {0,0} wavelet but it doesn't show them).The deeper you go into these coefficients the less wavelets are used to construct the wave and therefor they are more basic.In this case we want the highest quality data so we can just look at level 1. (So the coefficient named "{1}").So time to get down to the real code.Here is some example data, in this case a fourier transform.This is a fairly easy case but you can see that the number of datapoints is giant and there is some noise toward the bottom.So the first thing we need to do is get rid of all the little peaks we don't care about.Here "data" is the data with peaks in it and "min" is the minimum size of the peaks that will be detected.dwd=DiscreteWaveletTransform[If[# < min, 0, #] & /@ data, HaarWavelet[],  1];In this case I have set "min" to 1.We can now run the inverse wavelet transform of that for the first detail coefficient:InverseWaveletTransform[dwd, HaarWavelet[], {1}]Which looks like this:Now we detect every time those peaks cross zero.  You may also notice that there is another threshold here which I found to be useful.Normal[CrossingDetect[  If[Abs[#] < ther, 0, #] & /@    InverseWaveletTransform[    DiscreteWaveletTransform[If[# < min, 0, #] & /@ data,      HaarWavelet[], 1], HaarWavelet[], {1}]]]This gets the position of the points in each peak:Flatten[Position[  Normal[CrossingDetect[    If[Abs[#] < ther, 0, #] & /@      InverseWaveletTransform[      DiscreteWaveletTransform[If[# < min, 0, #] & /@ data,        HaarWavelet[], 1], HaarWavelet[], {1}]]], 1]]We then split them into groups for each peak:Split[Flatten[  Position[Normal[    CrossingDetect[     If[Abs[#] < ther, 0, #] & /@       InverseWaveletTransform[       DiscreteWaveletTransform[If[# < min, 0, #] & /@ data,         HaarWavelet[], 1], HaarWavelet[], {1}]]], 1]], Abs[#1 - #2] < cther &]"cther" is the distance 2 points have to be appart before they are counted as sepperate peaks.  In this case I use 50 because I don't want it detecting false peaks right next to the real ones.  But you could set it to split after every run of consecutive 1s.After this it is pretty technical and so I present (drumroll) the final function!!! FindPeaks[data_, ther_: .2, min_: 0, cther_: 50] :=   Function[peaks,     Transpose[{peaks,       data[[peaks]]}]][(#[[Ordering[data[[#]], -1][[1]]]] & /@      Split[Flatten[       Position[        Normal[CrossingDetect[          If[Abs[#] < ther, 0, #] & /@            InverseWaveletTransform[           DiscreteWaveletTransform[If[# < min, 0, #] & /@ data,             HaarWavelet[], 1], HaarWavelet[], {1}]]], 1]],      Abs[#1 - #2] < cther &])]So here is it working on the fourier transform!The blue is the unmodified data and the red is a line from peak to peak that fell within our requirments, but with different settings we could get all of these smaller peaks off to the right.Well that's all I got so if you have any questions about the code or edits that make it better, post it here!
5 years ago
17 Replies
 Very nice approach. I wanted to run the code to play around with it but I don’t think you gave the Fourier transform data. Any chance you can post that?
5 years ago
 How can I post it?  I can't just write out the data beacsue it is 5000 data points.
5 years ago
 I thought you are generating some points with a formula and then taking Fourier transform. But, yes, no way to post 5000 points of pure data ;-)
5 years ago
 Marcelo De Cicco 1 Vote Hi people,Just adding more ideas, I wrote a code that get the maxima/minima in a data range, for example a experimental Absortion light of a material, what it does: read always 3 points and when the middle point is greater or lesser than the two other we can have a maximum, or when it is lesser than the 2 others we can have a minima.Flatten[Take[   Transpose[{x, y}][[    1 + Flatten@      Position[       Partition[Transpose[{x, y}][[All, 2]], 3,         1], {a_, b_, c_} /; b < Min[a, c] || b > Max[a, c]]]], All,    1]];
5 years ago
 How fast is it when using a large dataset?
5 years ago
 I really do not know.Up to 100 points is really fast.
5 years ago
 I'm going to try 20,000.  Where do you pu tin the list?
5 years ago
 Hi,the code is very simple,  in fact is part of another bigger code. So I clean it see below:Flatten[Take[  "your data"[[   1 + Flatten@     Position[      Partition["your data"[[All, 2]], 3, 1], {a_, b_, c_} /;        b > Max[a, c]]]], All, 1]]But I think you should also include some constraints or conditional for the differences between the middle point and two others, otherwise you will get all small fluctuations.For example, you can think a constraint that only considerer differences grater than 1.2 or 2, depends on the scale that you are working with.
5 years ago
 I got this errorPart::partd: Part specification {0.00218645,0.000108345,0.000135432,0.000278259,0.000516786,0.000136859,0.000219704,0.000689318,0.000447265,0.000527103,0.000763647,0.000590564,0.000519036,0.000520092,<<24>>,0.000292454,0.000246734,0.000581422,0.000241731,0.000190411,0.000268275,0.00069401,0.000499151,0.000759404,0.000743745,0.000950466,0.000383596,<<64462>>} is longer than depth of object. >>Also have this looks similar to MinDetect[] and MaxDetect[] but min and max detect are painfully slow.
5 years ago
5 years ago
 Part isn't trying to extract pieces of the expression because I am usion a variable.Also probably kilo of mega but with wavelets you can get the transform for an image which can be giant almost instantlly.
5 years ago
 Christopher,To upload the datapoints have you tried using Compress? If the resultant string is not too long you can post it and people can Decompress it.
5 years ago
 Let me guess, are you analyzing sound waves to get the current tone (height)?
5 years ago
 Almost.  First givin a fourier transform it is actually the position on the x axis which determines the pitch.  And second, I am going to write a community post about some pitch recignition stuff soon .