Message Boards Message Boards

Wavelet-Based Peak Detection

Posted 11 years ago
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
17 Replies
Posted 8 years ago

So Chris, for fast pitch recognition use a Sliding DFT (SDFT) to calculate a continuous time-frequency spectrogram. The sliding DFT allows you to recursively generate your next spectrum bin using the previous bin giving you a new spectrum every sample. It also has the added benefit of providing any frequency resolution lower than the nyquist of course, without the size constraints & overhead of the FFT (64, 128,256,512,1024, ... , 2^n) requirements on size (due to the spider web multiplication). This would calculate your spectrum quickly & continuously. This is what I did for a project last year here & here using the STM32F7 discovery board. The main focus of that project was to implement just a RT audio spectrogram, but we were able to extend it to recognize the instrument pitch which is not as straight forward as peak detection, due to the timbre creating harmonics in the spectrum (the harmonic's amplitude sometimes is greater than the fundamental). We ended up using a harmonic product spectrum which basically takes the fundamental & multiplies it with it's respective harmonics. We did this for all the notes & their respective 5 harmonics. Then we just found the max value after the harmonic product spectrum, which for the most part worked nicely. Here's a link with the guitar note recognition which is a little rusty since I overhauled the note display portion with an all night-er coding link. Any way I hope that helps you, any questions feel free to ask, but before I forget, thank you for this write-up!!!! I think i can use this to identify the real peaks in the RT spectrum to better visualize the FM de-modulation for a better demonstration. Tell your dad i'm looking for a job to if yah could. ; ) LOL

POSTED BY: Nathan Imig
Posted 10 years 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
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]]]];
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 .
Let me guess, are you analyzing sound waves to get the current tone (height)?
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
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.
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
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.
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
I really do not know.

Up to 100 points is really fast.
POSTED BY: Marcelo De Cicco
I'm going to try 20,000.  Where do you pu tin the list?
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
How fast is it when using a large dataset?
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
How can I post it?  I can't just write out the data beacsue it is 5000 data points.
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
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract