Finding Local Minima / Maxima in Noisy Data

Question

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:

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 = 
    Flatten@Position[
       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]

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.

Answer 1

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.

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:

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]

Answer 2

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:

ClearAll[findExtrema];
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 = 
       MapThread[
         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:

ClearAll[displayPoints]
displayPoints[pts_, extr_] :=
  ListPlot[pts, 
     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}]]

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).

Answer 3

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]}]

Answer 4

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:

Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/Applications/QuantileRegressionForLocalExtrema.m"]

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

Show[{ListPlot[Join[{data3}, extrema3], 
   PlotStyle -> {{}, {PointSize[Medium], Red}, {PointSize[Medium], 
      Green}}],
  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:

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

Update

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).

Answer 5

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;
tempvaluelist = Sinc[temptimelist] + RandomReal[{-0.02, 0.02}, 200];
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):

Manipulate[
 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}]

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]

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: