5
$\begingroup$

I have experimental data which is noisy. My objective is to find the peaks. Here is a short part of the data:

data = Uncompress[FromCharacterCode[
  Flatten[ImageData[Import["http://i.stack.imgur.com/nCHC5.png"],"Byte"]]]];

A simple application of FindPeaks leads to the following:

peaks = FindPeaks[data, 50];
ListLinePlot[data, Epilog -> {Red, Point[peaks]}]

Mathematica graphics

The two peaks each have a double point. I would have thought there should be just one for each.

To explore this more closely I have put together the following DynamicModule to explore the use of FindPeaks

ClearAll[myFindPeaks];
SetAttributes[myFindPeaks, HoldFirst];
myFindPeaks[op_, data_] := 
 DynamicModule[{peakPts, σ = 0, s = 0, n1, n2, nn, n1p, n2p, 
   f1, f2, f0, finc},
    n1 = 1; n2 = nn = Length[data]; n1p = 0; n2p = 1;
  peakPts = FindPeaks[data, σ];
 Column[{
  Row[{"Gaussian Bluring " Slider[
    Dynamic[σ, (σ = #; 
       peakPts = FindPeaks[data, σ, s];) &], {0, 100, 1}], 
  Dynamic[σ],
  " Sharpness ", 
  Slider[Dynamic[
    s, (s = #; peakPts = FindPeaks[data, σ, s];) &], {0, 
    50}],
  Button["Output Peaks", op = peakPts]}],
Row[{"Zoom Interval", 
  IntervalSlider[
   Dynamic[{n1p, 
     n2p}, ({n1p, n2p} = #; n1 = Round[(nn - 1) n1p + 1]; 
      n2 = Round[(nn - 1) n2p + 1];) &], {0, 1}, 
   ImageSize -> {10 72, 20}] }],
   Dynamic[
   ListLinePlot[{data[[n1 ;; n2]], 
   GaussianFilter[data[[n1 ;; n2]], {3, σ}]}, 
  PlotRange -> All, Frame -> True, 
  FrameLabel -> {"Frequency/Hz", "FRF Modulus"},
  PlotLegends -> {"Data", "Smoothed Data"},
  Epilog -> {
    Red, PointSize[0.005], 
    Point[{#[[1]] - (n1 - 1), #[[2]]} & /@ peakPts],
    Black, 
    Table[Text[
      ToString[n], {#[[1]] - (n1 - 1), #[[2]]} &[
       peakPts[[n]]], {0, -2}], {n, Length[peakPts]}]
       },
    ImageSize -> 12 72]
    ]
    }]
   ]


myFindPeaks[pks, data]

Mathematica graphics

When the Gaussian blurring is zero the noise generates many peaks. At about 21 the noise has gone but the peaks at about 200 and 400 each have double peaks. When the Gaussian blurring reaches 58 the two peaks on the peak at about 400 disappear. Should not just one peak have disappeared? How do I find my peaks?

$\endgroup$
  • $\begingroup$ FindPeaks only searches for peaks, smoothing is you task. You can give to FindPeaks smoothed curve using MovingAverage, LowpassFiler, MedianFilter, etc. Thus you will remove small noise and get the accurate peaks. $\endgroup$ – funnyp0ny Nov 12 '14 at 22:10
  • $\begingroup$ MaxDetect[data, 0.01] where 0.01 is a "noise threshold" is also worth a try - the resulting vector is 1 for (possibly extended) maxima, and 0 elsewhere. $\endgroup$ – Niki Estner Nov 13 '14 at 7:05
  • $\begingroup$ @funnypony Thank you but I need to be convinced. Is Gaussian blurring not a form of smoothing? I will reply in more detail against your answer. $\endgroup$ – Hugh Nov 13 '14 at 11:09
  • $\begingroup$ @nikie Thanks. This is a new module to me. I will investigate. $\endgroup$ – Hugh Nov 13 '14 at 11:23
3
$\begingroup$

Going further with my comment:

peaks = FindPeaks[MedianFilter[data, 2]]
(*{{1, 0.253492}, {389/2, 0.489684}, {377, 0.313048}}*)

then you have 3 nice peaks. Instead of MedianFilter you can use LowpassFiler, MovingAverage, etc. - this won't disturb your data but will give you peaks positions.

$\endgroup$
  • $\begingroup$ It may turn out that it is necessary to do my own smoothing. However, I thought that FindPeaks had a smoothing system called Gaussian blurring. This does seem to work to some extent as illustrated by my DynamicModule. Also the examples make reference to Gaussian blurring. I show in my DynamicModule the curve smoothed by Gaussian blurring and the blurred curve only has one peak. I may well do the flitering as you suggest but the problem is that it will change my data. So if there is built in filtering why is it not working? Why do I have to do extra filtering? $\endgroup$ – Hugh Nov 13 '14 at 11:13
4
$\begingroup$

You may choose among many smoothing techniques. A way to do is using wavelets, which is not so different from my other answers [1,2]. You can use those with minor changes in threshold type and value. However, you say there are two peeks and there should be one instead. It's up to you to decide which one is the right peak to find. Below, I choose the one on the right (in your image its number 3). But if you need the left one just change SymletWavelet[5] to SymletWavelet[2].

dwd = DiscreteWaveletTransform[data, SymletWavelet[5], Padding -> "Fixed"];
(*Use SmoothGarrote threshold*)
dwdt = WaveletThreshold[dwd, {"SmoothGarrote", 0.005, 5}];
(*Get filtered values*)
smooth = InverseWaveletTransform[dwdt];
(*Find peeks of smoothed data*)
peaks = FindPeaks[smooth, 50]
{{1, 0.25371993}, {195, 0.49051315}, {378, 0.31279138}}
(*Plot data and filtered values*)
ListLinePlot[{data, smooth}, Epilog -> {Red, Point[peaks]}]

enter image description here

$\endgroup$
  • $\begingroup$ Thanks for this. The general suggestion is that I must do my own filtering (smoothing) and can't use the smoothing build into FindPeaks. I have curves with many peaks some close to each other. Smoothing is necessary but on a length scale that is smaller than my peak seperation. I will look at wavelet smoothing. $\endgroup$ – Hugh Nov 13 '14 at 11:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.