I'm trying to find local minima / maxima in noisy data, consisting of data values taken at certain time intervals. Ideally, the function should take a pair of lists (one containing time values and one containing observed data values) and return the coordinates of the maxima and minima.

An example of the data is found below:

temptimelist = Range[200]/10;
tempvaluelist = Sinc[#] &@temptimelist + RandomReal[{-1, 1}, 200]*0.02;

While the questions here, here and here have a good range of answers, they don't fit with my requirements since (as far as I can see) they determine whether data is a local maximum / minimum by comparing it with the adjacent values. If I apply such an algorithm, my extrema will be identified as follows:

data with peaks identified

What I've done is to write the following code.

NoisyExtremaFinder = Function[{timeList, valueList, aroundRange},
   (*NoisyExtremaFinder[] takes a pair of lists timeList_ and valueList_ and determines extrema in valueList_, returning a list containing the coordinates (as pairs of time_ and value_) of the minima as the first entry and the maxima as the second entry. 

   As the data is assumed to be noisy, we provide the option aroundRange_ to allow the user to determine the sensitivity of the search. Specifically, when aroundRange_=n, the function will compare each value with the preceding and subsequent n values to determine whether it is an extrema*)

   extremaPosition = 
       Map[#, Partition[valueList, 2*aroundRange + 1, 1, {-(1 + aroundRange), 1 + aroundRange}, {}]] - valueList, 0.] &;
    (*extremaPosition[] is a custom function that determines the position of local Maxima or Minima in valueList_, with the sensitivity determined by the value of aroundRange. When aroundRange_=n, the function will compare each value with the preceding and subsequent n values to determine whether it is an extrema. You can either do extremaPosition[Max] or extremaPosition[Min] *)

   extremaPoints = 
    Transpose@{timeList[[#]], valueList[[#]]} &@extremaPosition[#] &;
  (*extremaPoints[] is a custom function that determines the coordinates of local Maxima or Minima in valueList_

   Custom Functions Used: extremaPosition[]*)

   {extremaPoints[Min], extremaPoints[Max]}];

Basically, instead of deciding whether a point is an extrema by comparing it only with adjacent terms, it decides whether the point is an extrema by comparing with the aroundRange preceeding and aroundRange subsequent points.

Then, by adding the following lines of code, we can see the final results:

NoisyExtremaFinder[temptimelist, tempvaluelist, 10]; (*parameter "10" chosen by estimating the distance from peak to peak*)
ListPlot[Transpose[{temptimelist, tempvaluelist}], 
 Epilog -> {PointSize[Medium], Red, Point[extremaPoints[Min]], Green, 
   Point[extremaPoints[Max]]}, ImageSize -> 700]

better peaks

Is there any better way to code to achieve the goal I have in mind of comparing a point with n preceeding and subsequent points?

I was also advised by rm-rf (in chat) to consider smoothing the noisy data first before trying to find the local minima / maxima. However, I was concerned that I could add unwanted artifacts into the data by the smoothing process. With regards to smoothing, is there an algorithm that is commonly used to filter experimental data? Filters that I have looked at include the Savitsky-Golay filter and the Low-Pass filter.

PS: My data sets have around 10,000 data points each.

  • $\begingroup$ Have you tried MovingAverage? I've used this in the past, successfully to smooth noisy data. $\endgroup$ – Jonathan Shock Apr 23 '13 at 8:40
  • 1
    $\begingroup$ @JonathanShock Thanks for the reminder - I'll try that, and it perhaps could work with my current data set. However, I'm quite sure that it wouldn't work out for data with abrupt changes (to quote an example, velocity data during collisions, where the velocity changes rapidly from positive to negative) $\endgroup$ – Vincent Tjeng Apr 23 '13 at 8:54
  • $\begingroup$ You should know that Sinc[] is listable, and you already know that RandomReal[] can take ranges, so: tempvaluelist = Sinc[temptimelist] + RandomReal[{-0.02, 0.02}, 200]; $\endgroup$ – J. M. will be back soon Apr 23 '13 at 11:26
  • 1
    $\begingroup$ Is the noise really uniformly distributed or was that just an arbitrary choice? Would it be reasonable to assume Gaussian distribution? $\endgroup$ – rm -rf Apr 23 '13 at 13:28
  • 1
    $\begingroup$ @rm-rf I'm not exactly sure how the noise is distributed, but I guess a Gaussian distribution would make more sense. Don't want to make any assumptions as yet though - do you think I should upload my actual data for people to try? $\endgroup$ – Vincent Tjeng Apr 25 '13 at 1:29

An alternative way of doing things is using wavelets. Wavelets are quite good at denoising.

(*Define data and noise*)
temptimelist = Range[200]/10;
data = Sinc[temptimelist];
noise = RandomReal[{-0.02, 0.02}, 200];

(* Define Wavelet for denoising *)
dwd = DiscreteWaveletTransform[data + noise, SymletWavelet[7], 6];
(* Use universal threshold *)
dwd = WaveletThreshold[dwd];
(*Get denoişed values*)
est = InverseWaveletTransform[dwd];

ListLinePlot[{data, est, (data + noise)}, PlotRange -> All]

Once you denoised the data you can find extreme values described in other answers. wavelet-transformed data

For example check this from @Theodros Zelleke answer to find the extremal points in a list.

findExtremaPos[list_List] := 
 Module[{signs, extremaPos, minPos, maxPos}, 
 signs = Sign[Differences[list]];
 signs = signs //. {a___, q_, 0, z__} -> {a, q, q, z};
 extremaPos = 1 + Accumulate@(Length /@ Split[signs]);
 If[First@signs == 1, minPos = extremaPos[[2 ;; -2 ;; 2]]; 
 maxPos = extremaPos[[1 ;; -2 ;; 2]], 
 minPos = extremaPos[[1 ;; -2 ;; 2]]; 
 maxPos = extremaPos[[2 ;; -2 ;; 2]]];
 {minPos, maxPos}]

(*Evaluate extreme points*)
{minPos, maxPos} = findExtremaPos[est];

ListPlot[est, Joined -> True, 
 Epilog -> {PointSize[Large], Red, Point[{#, est[[#]]} & /@ minPos], 
 Blue, Point[{#, est[[#]]} & /@ maxPos]}, PlotRange -> All]

The result:

enter image description here

Edit 1

With the version 10 one can use FindPeaks to find peaks easily.

peaks = FindPeaks[est, .5]

{{2, 1.0147576}, {6, 0.94715655}, {78, 0.13402331}, {143, 0.066621946}, {192, 0.04389471}, {197, 0.041374537}, {199, 0.060916925}}

ListPlot[{est, peaks}, Joined -> {True, False}, 
  PlotStyle -> {Automatic, {PointSize[.01], Red}}, PlotRange -> All]

enter image description here


Here is an alternative based on LinearModelFit and polynomial regression. The idea is to fit your data to a polynomial of some degree, to find approximate extrema points of the smoothed data, and then locate nearest extremal points. Here is the code:

findExtrema[points_List, fitOrder_Integer, around_Integer: 5] :=
   Module[{fit, fn, extrema, x, signs, extPoints},
     fit = LinearModelFit[points, x^Range[0, fitOrder], x];
     fn = fit["Function"];
     extrema = NSolve[fn'[x] == 0, x, Reals];
     signs = Sign[fn''[x] /. extrema];
     extPoints = 
         First@SortBy[Nearest[points, #, around], #2] &,
         {{x, fn[x]} /. extrema, signs /.{1 -> Last, -1 -> (-Last[#] &)}}
     Map[Pick[extPoints , signs, #] &, {1, -1}]

What happens here is that we fit the data with a polynomial of degree fitOrder, then find extremal points as zeros of the derivative, and determine the types of extrema by the sign of the second derivative. Having found extrema for the polynomial, we then find the nearest extremal points in the dataset.

Here is a display function then:

displayPoints[pts_, extr_] :=
     Epilog -> {PointSize[Medium], Red, Point[First@extr], Green,Point[Last@extr]},
     ImageSize -> 700

which we can use as

displayPoints[#, findExtrema[#, 10, 10]] &[Transpose[{temptimelist, tempvaluelist}]]

enter image description here

You can play with the parameters fitOrder and around. Basically, fitOrder needs to be equal or greater than the number of "features" in your data (number of intervals where it is roughly monotonous). It is however not advisable to set fitOrder to very high numbers, so this method will work reaonably well for relatively smooth data. If your data oscillates rapidly, it may make sense to pick a different set of basis functions for regression, such as trigonometric functions (Sin, Cos).

  • $\begingroup$ "polynomial data fitting" - so, generalized Savitzky-Golay, then? :) $\endgroup$ – J. M. will be back soon Apr 23 '13 at 16:45
  • $\begingroup$ @J.M. I know this under the name of polynomial regression, but that is probably the same thing. Regression is not my field, so you probably know better. $\endgroup$ – Leonid Shifrin Apr 23 '13 at 16:49
  • $\begingroup$ @J.M. However, what I suggest is not a usual filter or smoothing kernel, so my method is free of at least some artifacts that can be induced by smoothing. $\endgroup$ – Leonid Shifrin Apr 23 '13 at 16:52
  • $\begingroup$ Hmm, true, you didn't take a window, but you fit to the entire data set. I suppose we should mention the warning not to set fitOrder too high, tho. $\endgroup$ – J. M. will be back soon Apr 23 '13 at 16:57
  • $\begingroup$ @J.M. Yes, that's a good suggestion, will add. In fact, the code above is trivial to generalize to a different set of basis functions, such as trig ones. So, if data is highly oscillatory, that could be a better option. $\endgroup$ – Leonid Shifrin Apr 23 '13 at 16:59

The problem you face is very common in signal- and imageprocessing.

I'm trying to find local minima / maxima in noisy data

Since all the noise introduces small local extrema, the very question is which extrema are still noise and which are signal. What you want to do is to smooth out the noise without loosing signal information.

This task is one of the most common things in signal- and imageprocessing and there are at least four hundred and sixty three million ways to do it. This variety comes from the fact, that the more you know about your noise and your signal, the better you can adjust/choose the smoothing procedure.

A well known method is the Gaussian filter and if you really don't know where to start, use this until you gained more experience. What I propose are two simple steps. First you filter your data with a GaussianFilter and then you interpolate the smoothed data and select the points where the derivative vanishes. Just as you would have done it in school

We start with your sample data

ts = Range[200]/10.;
vals = Sinc[#] &@ts + RandomReal[{-1, 1}, 200]*0.02;

The we interpolate your data, where we have filtered the vals. The magic number here is the filter width of the Gaussian filter. I have chosen 5 here. You have to choose in your real data a value wich removes the noise but leaves the local extrema intact

ip = Interpolation[Transpose[{ts, GaussianFilter[vals, 5]}]]

If you look at your data, the smoothed interpolation and the derivative ip'[t], you see that ip' has zero crossings only where your data has extreme points. The final step is to find each pair {ts[[n]],ts[[n+1]]} of your time data where the derivative of the first time has a different Sign than the second. This indicates, that we crossed the x-axis.

extr = First /@ Select[Partition[ts, 2, 1], Sign[ip'[#[[1]]]] != Sign[ip'[#[[2]]]] &]
(* {4.5, 7.7, 10.8, 14.2, 17., 19.7} *)

Finished and ready to be plotted

ListPlot[Transpose[{ts, vals}],
 Epilog -> {Red, PointSize[.02], Point[{#, ip[#]} & /@ extr]}]

enter image description here

  • $\begingroup$ thank you, and I'll read up more on the Gaussian filter. After running it a few times, I've found that your answer finds the time-coordinate of the "actual" extrema very well (by "actual" extrema I mean the extrema in the data without noise), but not so much for the value-coordinate. $\endgroup$ – Vincent Tjeng Apr 26 '13 at 7:08
  • $\begingroup$ @VincentTjeng Please note, that I used ip[#] for the value because it was short. You must not do that! What you want is the real value from your vals list. You can utilize Position the get the position of the found time from the ts list and then you use this position to extract the correct value from the vals list. $\endgroup$ – halirutan Apr 26 '13 at 10:29
  • $\begingroup$ Yep, noted. I was going to mention that to you at first, but I realized that, ip[#] often gives a better estimator of the "actual" extrema than the real value from the vals list (and I thought you intentionally chose ip[#] because of this reason.) $\endgroup$ – Vincent Tjeng Apr 26 '13 at 10:34
  • $\begingroup$ Yes and no. Since your data is so noisy, the exact value from vals will most probably be useless. On the other hand for ip we smoothed quite a bit. Maybe you want something in between. Take the exact value from vals and collect say 3 surrounding values too. Than you calculate the Median which probably gives better results than ip (which uses basically a weighted Mean). Everything needs then to be adjusted for your real data. $\endgroup$ – halirutan Apr 26 '13 at 10:53

I also got interested in this problem and solved it using Quantile Regression. See my blog post Finding local extrema in noisy data using Quantile Regression .

The proposed Quantile Regression algorithm is a version of the polynomial fitting solution proposed by Leonid Shifrin above, and has the following advantages: (i) it requires less parameter tweaking, and (ii) more importantly it is more robust with very noisy and oscillating data.

That robustness is achieved by using two regression fitted curves: one close to the local minima and another close to the local maxima, computed for low and high quantiles respectively. Quantile Regression is uniquely able to do that -- I look for good examples where this feature is of decisive importance.

Here is an example of very noisy and oscillating data that has 10,000 points:

n = 1000;
xs = N@Rescale[Range[n], {1, n}, {0, 60}];
data3 = Flatten[
   Table[Transpose[{xs + 0.1 RandomReal[{-1, 1}, Length[xs]], 
      Map[5 Sinc[#] + Sin[#] + 4 Sin[1/4 #] &, xs] + 
       1.4 RandomVariate[SkewNormalDistribution[0, 1, 0.9], 
         Length[xs]]}], {10}], 1];
ListPlot[data3, PlotRange -> All]

Here is code for finding the regression quantiles functions, local extrema, and making a plot like the ones below:


{qfuncs3, extrema3} = QRFindExtrema[data3, 24, 2, 150, {0.02, 0.98}];

Show[{ListPlot[Join[{data3}, extrema3], 
   PlotStyle -> {{}, {PointSize[Medium], Red}, {PointSize[Medium], 
  Plot[Through[qfuncs3[x]], {x, Min[data3[[All, 1]]], 
    Max[data3[[All, 1]]]}, PlotStyle -> {Orange}, 
   PerformanceGoal -> "Speed"]}, Frame -> True, PlotRange -> All]

This image shows the data, the fitted regression quantiles, and the estimated local extrema:

enter image description here

Here is another example with more timid data (again ~ 10,000 points):

enter image description here


I would like to point out that since we find two regression quantile curves we can use two nearest neighbors finding functions: one with the points below the low regression quantile, and one with the points above the high regression quantile.

I changed the implementation linked above to also take an option should the Nearest functions for finding the extrema be constructed using all data points or just the outliers (the points outside of the found regression quantiles).

  • $\begingroup$ Thanks for posting this it looks powerful. You have only posted the code for generating the data. Please can you post the code for finding the regression quantiles and the extrema? $\endgroup$ – Hugh Sep 27 '15 at 21:01
  • $\begingroup$ @Hugh Thanks! I updated the answer with such code. $\endgroup$ – Anton Antonov Sep 27 '15 at 21:40
  • $\begingroup$ I haven't had a chance to look fully at the code, but this definitely looks interesting. Thanks for thinking about this problem so long after it was posted! $\endgroup$ – Vincent Tjeng Sep 29 '15 at 3:37
  • $\begingroup$ @VincentTjeng This question was linked to in mathematica.stackexchange.com/questions/61026/… . I thought it would be interesting to see the application of Quantile Regression to it... $\endgroup$ – Anton Antonov Sep 29 '15 at 11:27

Maybe you can undersample your data according to a mean as suggested in the comments, and then try and use one of the methods from the links you provide. The following is semi-manual but it would work reasonably well and quick for a dataset like the one you provide. Starting with the data:

temptimelist = Range[200]/10;
originalData = Transpose[{temptimelist, tempvaluelist}];

you can then check what sampling/averaging works for you (i.e. in the attached I average over six points and undersamplle by a factor of six):

 nearestFunction = 
  Mean@Nearest[Thread[originalData[[All, 1]] -> originalData[[All, 2]]], #, a] &;
 (*nearest function that averages over "a" closest points*)
 newdata = {N@#, Sequence @@ (nearestFunction@#)} & /@ 
   Range[Min@originalData[[All, 1]], 
    Max@originalData[[All, 1]], ((Max@originalData[[All, 1]] - Min@originalData[[All, 1]])/
       Length@originalData[[All, 1]] ) b];
 (*nearest function applied to an undersampling of your data. b = 10 means a tenth of data used etc *)
 ListPlot[{originalData, newdata},
   PlotRange -> All, 
  PlotStyle -> {Directive[PointSize[0.005], Blue, PlotJoined -> True],
     Directive[PointSize[0.014], Red]}], {a, 1, 10, 1}, {b, 1, 20, 1}]

enter image description here

and if that works you can use either one of the other posts or interpolate.

Plot[Interpolation[newdata][x], {x, 0.1, 20}, PlotRange -> All]

enter image description here

from which is easy to find extrema with a variety of methods.

The above has the problem that overaveraging will lose you some of the features of the original function, i.e. here the maximum near zero and near twenty are lost:

enter image description here


Your Answer

