Find a local maximum in list without using built-in function FindPeaks
up vote
3
down vote
favorite
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]}
]
list-manipulation
add a comment |
up vote
3
down vote
favorite
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]}
]
list-manipulation
add a comment |
up vote
3
down vote
favorite
up vote
3
down vote
favorite
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]}
]
list-manipulation
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]}
]
list-manipulation
list-manipulation
edited Nov 24 at 17:14
kglr
173k8195400
173k8195400
asked Nov 24 at 16:21
Kiril Danilchenko
714316
714316
add a comment |
add a comment |
2 Answers
2
active
oldest
votes
up vote
4
down vote
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.
New contributor
If your data has noise, defining a peak as a point that is larger than its neighbors might not be useful. In the presence of noise, you could a) smooth the data, b) have a min threshold between the peak point and the closest neighbor, and/or c) require a peak point be a peak in some window in which points on the left increase while points on the right decrease (with some heuristic for constant sequences). Note that (c) is different from MassDefect's maximum over some range of values method. Other methods exist.
– Robert Jacobson
Nov 24 at 19:00
@RobertJacobson You're absolutely right. My algorithm is pretty weak. I'm not entirely clear on what kind of peaks the poster is looking for.
– MassDefect
Nov 25 at 22:53
add a comment |
up vote
4
down vote
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.
add a comment |
2 Answers
2
active
oldest
votes
2 Answers
2
active
oldest
votes
active
oldest
votes
active
oldest
votes
up vote
4
down vote
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.
New contributor
If your data has noise, defining a peak as a point that is larger than its neighbors might not be useful. In the presence of noise, you could a) smooth the data, b) have a min threshold between the peak point and the closest neighbor, and/or c) require a peak point be a peak in some window in which points on the left increase while points on the right decrease (with some heuristic for constant sequences). Note that (c) is different from MassDefect's maximum over some range of values method. Other methods exist.
– Robert Jacobson
Nov 24 at 19:00
@RobertJacobson You're absolutely right. My algorithm is pretty weak. I'm not entirely clear on what kind of peaks the poster is looking for.
– MassDefect
Nov 25 at 22:53
add a comment |
up vote
4
down vote
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.
New contributor
If your data has noise, defining a peak as a point that is larger than its neighbors might not be useful. In the presence of noise, you could a) smooth the data, b) have a min threshold between the peak point and the closest neighbor, and/or c) require a peak point be a peak in some window in which points on the left increase while points on the right decrease (with some heuristic for constant sequences). Note that (c) is different from MassDefect's maximum over some range of values method. Other methods exist.
– Robert Jacobson
Nov 24 at 19:00
@RobertJacobson You're absolutely right. My algorithm is pretty weak. I'm not entirely clear on what kind of peaks the poster is looking for.
– MassDefect
Nov 25 at 22:53
add a comment |
up vote
4
down vote
up vote
4
down vote
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.
New contributor
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.
New contributor
New contributor
answered Nov 24 at 16:52
MassDefect
1312
1312
New contributor
New contributor
If your data has noise, defining a peak as a point that is larger than its neighbors might not be useful. In the presence of noise, you could a) smooth the data, b) have a min threshold between the peak point and the closest neighbor, and/or c) require a peak point be a peak in some window in which points on the left increase while points on the right decrease (with some heuristic for constant sequences). Note that (c) is different from MassDefect's maximum over some range of values method. Other methods exist.
– Robert Jacobson
Nov 24 at 19:00
@RobertJacobson You're absolutely right. My algorithm is pretty weak. I'm not entirely clear on what kind of peaks the poster is looking for.
– MassDefect
Nov 25 at 22:53
add a comment |
If your data has noise, defining a peak as a point that is larger than its neighbors might not be useful. In the presence of noise, you could a) smooth the data, b) have a min threshold between the peak point and the closest neighbor, and/or c) require a peak point be a peak in some window in which points on the left increase while points on the right decrease (with some heuristic for constant sequences). Note that (c) is different from MassDefect's maximum over some range of values method. Other methods exist.
– Robert Jacobson
Nov 24 at 19:00
@RobertJacobson You're absolutely right. My algorithm is pretty weak. I'm not entirely clear on what kind of peaks the poster is looking for.
– MassDefect
Nov 25 at 22:53
If your data has noise, defining a peak as a point that is larger than its neighbors might not be useful. In the presence of noise, you could a) smooth the data, b) have a min threshold between the peak point and the closest neighbor, and/or c) require a peak point be a peak in some window in which points on the left increase while points on the right decrease (with some heuristic for constant sequences). Note that (c) is different from MassDefect's maximum over some range of values method. Other methods exist.
– Robert Jacobson
Nov 24 at 19:00
If your data has noise, defining a peak as a point that is larger than its neighbors might not be useful. In the presence of noise, you could a) smooth the data, b) have a min threshold between the peak point and the closest neighbor, and/or c) require a peak point be a peak in some window in which points on the left increase while points on the right decrease (with some heuristic for constant sequences). Note that (c) is different from MassDefect's maximum over some range of values method. Other methods exist.
– Robert Jacobson
Nov 24 at 19:00
@RobertJacobson You're absolutely right. My algorithm is pretty weak. I'm not entirely clear on what kind of peaks the poster is looking for.
– MassDefect
Nov 25 at 22:53
@RobertJacobson You're absolutely right. My algorithm is pretty weak. I'm not entirely clear on what kind of peaks the poster is looking for.
– MassDefect
Nov 25 at 22:53
add a comment |
up vote
4
down vote
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.
add a comment |
up vote
4
down vote
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.
add a comment |
up vote
4
down vote
up vote
4
down vote
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.
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.
edited Nov 24 at 20:01
answered Nov 24 at 19:36
C. E.
49k395197
49k395197
add a comment |
add a comment |
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f186630%2ffind-a-local-maximum-in-list-without-using-built-in-function-findpeaks%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown