13

Function argrelextrema from scipy.signal does not detect flat extrema. Example:

import numpy as np
from scipy.signal import argrelextrema
data = np.array([ 0, 1, 2, 1, 0, 1, 3, 3, 1, 0 ])
argrelextrema(data, np.greater)
(array([2]),)

the first max (2) is detected, the second max (3, 3) is not detected.

Any workaround for this behaviour? Thanks.

Marcus
  • 159
  • 1
  • 1
  • 8
  • Depending on the granularity of the data, a simple practical work-around can be to add a sort of "deterministic noise" to resolve ties in finding the min/max. I.e., internally using `data + np.linspace(-1e-5, +1e-5, len(data))`. Using linspace in an increasing/decreasing way allows to control in which direction ties should be resolved. – bluenote10 Dec 30 '22 at 10:05

2 Answers2

26

Short answer: Probably argrelextrema will not be flexible enough for your task. Consider writing your own function matching your needs.

Longer answer: Are you bound to use argrelextrema? If yes, then you can play around with the comparator and the order arguments of argrelextrema (see the reference).

For your easy example, it would be enough to chose np.greater_equal as comparator. 

>>> data = np.array([ 0, 1, 2, 1, 0, 1, 3, 3, 1, 0 ])
>>> print(argrelextrema(data, np.greater_equal,order=1))
(array([2, 6, 7]),)

Note however that in this way

>>> data = np.array([ 0, 1, 2, 1, 0, 1, 3, 3, 4, 1, 0 ])
>>> print(argrelextrema(data, np.greater_equal,order=1))
(array([2, 6, 8]),)

behaves differently that you would probably like, finding the first 3 and the 4 as maxima, since argrelextrema now sees everything as a maximum that is greater or equal to its two nearest neighbors. You can now use the order argument to decide to how many neighbors this comparison must hold - choosing order=2 would change my upper example to only find 4 as a maximum.

>>> print(argrelextrema(data, np.greater_equal,order=2))
(array([2, 8]),)

There is, however, a downside to this - let's change the data once more:

>>> data = np.array([ 0, 1, 2, 1, 0, 1, 3, 3, 4, 1, 5 ])
>>> print(argrelextrema(data, np.greater_equal,order=2))
(array([ 2, 10]),)

Adding another peak as a last value keeps you from finding your peak at 4, as argrelextrema is now seeing a second-neighbor that is greater than 4 (which can be useful for noisy data, but not necessarily the behavior expected in all cases).

Using argrelextrema, you will always be limited to binary operations between a fixed number of neighbors. Note, however, that all argrelextrema is doing in your example above is to return n, if data[n] > data[n-1] and data[n] > data[n+1]. You could easily implement this yourself, and then refine the rules, for example by checking the second neighbor in case that the first neighbor has the same value.

For the sake of completeness, there seems to be a more elaborate function in scipy.signal, find_peaks_cwt. I have however no experience using it and can therefore not give you more details about it.

pascal
  • 456
  • 4
  • 12
  • 1
    For other peak detection algorithms, you can also check https://github.com/MonsieurV/py-findpeaks `argrelextrema` is indeed not the easier choice when you want to tinker with algorithm behavior. – Yoan Tournade Feb 10 '18 at 13:09
0

I'm really surprised that no one figured out an answer to this. All you need to do is preprocess the array to remove duplicates that are located next to each other and you can run argrelextrema like so:

import numpy as np
from scipy.signal import argrelextrema
data = np.array([ 0, 1, 2, 1, 0, 1, 3, 3, 1, 0 ])

filter_table = [False] + list(np.equal(data[:-1], data[1:]))
data = np.array([x for idx, x in enumerate(data) if not filter_table[idx]])
argrelextrema(data, np.greater)
Ian
  • 53
  • 9
  • A practical problem of this solution is that typically the indices returned from `argrelextrema` are supposed to keep the semantics of the original indices. If you filter duplicates the semantics of the indices will change, and the result may become useless in many use cases. – bluenote10 Dec 30 '22 at 09:56
  • This is a fair point, but in the problem I was working on when I came across this it wasn't the case. – Ian Dec 30 '22 at 14:45