Group Abstract Group Abstract

Message Boards Message Boards

Wavelet-Based Peak Detection

GROUPS:
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 example

Here 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!

emoticon
POSTED BY: Christopher Wolfram
Answer
1 year ago
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?
POSTED BY: Vitaliy Kaurov
Answer
1 year ago
How can I post it?  I can't just write out the data beacsue it is 5000 data points.
POSTED BY: Christopher Wolfram
Answer
1 year 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 ;-)
POSTED BY: Vitaliy Kaurov
Answer
1 year ago
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]];
POSTED BY: Marcelo De Cicco
Answer
1 year ago
How fast is it when using a large dataset?
POSTED BY: Christopher Wolfram
Answer
1 year ago
I really do not know.

Up to 100 points is really fast.
POSTED BY: Marcelo De Cicco
Answer
1 year ago
I'm going to try 20,000.  Where do you pu tin the list?
POSTED BY: Christopher Wolfram
Answer
1 year 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.
POSTED BY: Marcelo De Cicco
Answer
1 year ago
I got this error
Part::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.
POSTED BY: Christopher Wolfram
Answer
1 year ago
How big is your data (mega, gigas,etcc)?

See the following link: http://mathematica.stackexchange.com/questions/14645/prevent-part-from-trying-to-extract-parts-of-symbolic-expressions
POSTED BY: Marcelo De Cicco
Answer
1 year 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.
POSTED BY: Christopher Wolfram
Answer
1 year 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.
POSTED BY: Hector Zenil
Answer
1 year ago
Let me guess, are you analyzing sound waves to get the current tone (height)?
POSTED BY: Vladimir Grankovsky
Answer
1 year 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 emoticon .
POSTED BY: Christopher Wolfram
Answer
1 year ago
I was doing some pitch recognition with Wavelet transforms, but it was a long time ago and I didn't know Mathematica well. It was very slow. But anyway, I was using this for peak detection (very empirical):
by the way, a useful idea is sample dropping - taking for example every 5th or 15th sample from the initial data (but mind the nyquist frequency)

 cwd=ContinuousWaveletTransform[signal2, GaborWavelet[15], {7, 30},
  SampleRate -> sr];
 freq = (#1[[1]] -> sr/#1[[2]] &) /@ cwd["Scales"];
 
 Maxima[scale_, treshold_: 0.7] :=
  Module[{i, max = -\[Infinity], maxpoz = 0},
   If[Max[scale] <= treshold, Return[{}];];
   For[i = Length[scale], i >= 0, i--,
    If[max < scale[[i]],
    max = scale[[i]];
    maxpoz = i;
    ];
   If[max - scale[[i]] >= treshold,
    Return[{maxpoz}];
    ];
   ];
  {}
  ]

frequencies =
  With[{list =
       2 freq[[Maxima[Abs /@ cwd[All][[50 ;;, 2, #]], 0.1] + 49]][[
         All, 2]]},
     If[Length[list] >= 1, list[[-1]], Null]
     ] & /@ Range[Length[cwd[All][[1, 2]]]];
POSTED BY: Vladimir Grankovsky
Answer
1 year ago
When using the Haar wavelet and only one detail level, isn't this equivalent to finding the crossings of the (upsampled) first differences, which should be efficient too?
POSTED BY: Rui Rojo
Answer
6 months ago