Ytukyg

I am interested finding the peaks in data without using the function FindPeaks. Any suggestion on how to do it effectively?

data1=RandomInteger[{1, 3997670}, 3936];

peaks = FindPeaks[data1]

ListLinePlot[

 data1, 

 Epilog -> {Red, PointSize[0.01], Point[peaks]}

 ]

You could try something like:

pks = {};



If[data1[[1]] > data1[[2]],

  AppendTo[pks, {1, data1[[1]]}]];



For[i = 2, i <= Length[data1] - 1, i++,

  If[data1[[i - 1]] < data1[[i]] > data1[[i + 1]],

   AppendTo[pks, {i, data1[[i]]}]]];



If[data1[[-1]] > data1[[-2]],

  AppendTo[pks, {Length[data1], data1[[-1]]}]];

It's a simple and rather inelegant way of finding peaks. Basically, if a point is higher than the two points on either side, it's considered to be a peak. In the case of the first and last point, it's a peak if it's higher than the closest point. I'm not sure if that's what you mean by finding all the peaks.

Is there a reason you can't use FindPeaks? It's almost always the best way. I didn't use it in my answer because you asked for an answer without, but I can get the exact same output with the following code:

pks2 = FindPeaks[data1, 0, 0, -Infinity];

If you want to find the local maximum over some range of values, you could do something like:

Max[data1[[100;;200]]]

It's position in the list can be found with:

Position[data1, Max[data1[[100;;200]]]]

Since you have a list, this will return the absolute local maximum. Calling Max on the entire list would of course return the global max.

The documentation describes a lot of how FindPeaks works, especially the section Scope > Parameters. It can all be done very efficiently with convolutions and vectorized operations.

This is basic peak finding:

findPeaks1[data_] := Module[{rightDiff, leftDiff, peaks},

  rightDiff = ListConvolve[{0, 1, -1}, data];

  leftDiff = ListConvolve[{-1, 1, 0}, data];

  Unitize[Sign[leftDiff] + Sign[rightDiff] - 2]

  ]

The reason for the 1 at the end of the name will become clear later. The convolutions in this function perform almost the same job as Differences, in the end, we get a list that is 0 in those positions that correspond to values that are larger than both its neighbors.

findPeaks1[data_, sigma_] := Module[{blurred},

  blurred = GaussianFilter[data, {3, sigma}];

  Unitize[findPeaks1[data] + findPeaks1[blurred]]

  ]

In this step, we add the second argument, which corresponds to smoothing the data using a Gaussian filter with variance sigma. A zero in this list means that the corresponding element is a peak, using the definition in the previous function, in both the smoothed data and the original data.

findPeaks1[data_, sigma_, s_] := Module[{sharpness},

  sharpness = ListConvolve[{1, -2, 1}, data];

  Unitize[findPeaks1[data, sigma] + 1 - Unitize[sharpness, s]]

  ]

In this piece of code, we implemented the third argument, which says that the peaks that don't have a negative second derivative that is larger than s should be discarded. The convolution we use here is the definition for the negative discrete second derivative.

findPeaks1[data_, sigma_, s_, t_] := Unitize[

  findPeaks1[data, sigma, s] + Unitize[UnitStep[Most@Rest@data - t] - 1]

  ]

Finally, in this code block, we implement the fourth argument, which is a threshold. It says that the value of a peak must be larger than t.

For efficiency reasons and because it made the implementation short and simple, these functions all pass lists of integers to each other. They have a 1 suffixed to their names to distinguish them from the interface which we shall now define:

outpformat[data_, sel_] := Module[{pos, values},

  pos = Position[sel, 0] + 1;

  values = Extract[data, pos];

  Transpose[{Flatten[pos], values}]

  ]



findPeaks[data_] := outpformat[data, findPeaks1[data]]

findPeaks[data_, sigma_] := outpformat[data, findPeaks1[data, sigma]]

findPeaks[data_, sigma_, s_] := outpformat[data, findPeaks1[data, sigma, s]]

findPeaks[data_, sigma_, s_, t_] := outpformat[data, findPeaks1[data, sigma, s, t]]

We now have a function findPeaks that does more or less the same thing as FindPeaks. I find that the result differs a bit for the Gaussian filter, they seem to do it a bit differently, but this still captures the idea.

Of course, there are many other signal processing techniques that we could use to find peaks, but there is no universal answer to how the signal should be processed. If you find that this function is not enough, perhaps because the signal is noisy or something else, then you have to look at the signal to decide what filters and such might be appropriate to make the peaks you are looking for more distinguishable.