Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

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.

share|improve this question
    
Have you tried MovingAverage? I've used this in the past, successfully to smooth noisy data. –  Jonathan Shock Apr 23 '13 at 8:40
    
@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) –  Vincent Tjeng Apr 23 '13 at 8:54
    
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]; –  J. M. Apr 23 '13 at 11:26
    
@J.M. thanks! sometimes when I edit code I leave all the silly bits in :) –  Vincent Tjeng Apr 23 '13 at 12:32
1  
Is the noise really uniformly distributed or was that just an arbitrary choice? Would it be reasonable to assume Gaussian distribution? –  rm -rf Apr 23 '13 at 13:28

4 Answers 4

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

share|improve this answer
    
Thank you for your answer. Just wondering, could you explain to me what the wavelets do? I'm a bit uncomfortable using smoothing techniques that I don't understand. –  Vincent Tjeng Apr 26 '13 at 6:38
    
you can check the help of mathematica there under applications you can see the denoising part. reference.wolfram.com/mathematica/ref/… –  s.s.o May 2 '13 at 11:16
    

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

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

share|improve this answer
    
"polynomial data fitting" - so, generalized Savitzky-Golay, then? :) –  J. M. Apr 23 '13 at 16:45
    
@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. –  Leonid Shifrin Apr 23 '13 at 16:49
    
@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. –  Leonid Shifrin Apr 23 '13 at 16:52
    
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. –  J. M. Apr 23 '13 at 16:57
    
@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. –  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

share|improve this answer
    
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. –  Vincent Tjeng Apr 26 '13 at 7:08
    
@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. –  halirutan Apr 26 '13 at 10:29
    
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.) –  Vincent Tjeng Apr 26 '13 at 10:34
    
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. –  halirutan Apr 26 '13 at 10:53

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

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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