Improving algorithm to find local maxima of 1D numpy array

Question

I am trying to find the local maxima of the function f(x) = (sin(x)/x)^2. For approximate solutions, I initialised two variables x and y and first plotted a graph to have a visual representation.

x = np.linspace(-20.0, 20.0, num=int((20+20)/0.01)) 
y = np.power(np.sin(x)/x, 2)
plt.plot(x, y, 'r.', markersize= 1)
plt.show()

This shows graph.

I then tried to create an algorithm to find the Maxima:

def returnMaxima(num, x, y):
    """
    number, np.array, np.array -> list
    num: number of maxima needed | x: x 1D array | y: y 1D array
    returns [[x1,y1], [x2,y2]...] in descending order of y
    """
    allMaximaPoints = [] # stores all Maxima points
    reqMaximaPoints = [] # stores num Maxima points
    for i in range(y.size): 
        # for first y value
        if i == 0: 
            if y[i] > y[i+1]:
                allMaximaPoints += [[x[i], y[i]], ]
        # for last y value
        elif i == y.size - 1:
            if y[i] > y[i-1]:
                allMaximaPoints += [[x[i], y[i]], ]
        # for rest y values
        else: 
            if y[i] > y[i-1] and y[i] > y[i+1]:
                allMaximaPoints += [[x[i], y[i]], ]
    # extract largest maximas from allMaximaPoints
    while num > 0: 
        reqMaximaPoints += [max(allMaximaPoints, key=lambda item:item[1]),]
        del allMaximaPoints[allMaximaPoints.index(max(allMaximaPoints, key=lambda item:item[1]))]
        num -= 1
    return reqMaximaPoints

When I tried returnMaxima(2, x, y) I get [[-4.4961240310077528, 0.04719010162459622], [4.4961240310077528, 0.04719010162459622]].

This is incorrect as it skipped the local maxima at x = 0. I suspect it is because the y[i-1] and y[i+1] values adjacent to the maxima at y[i] at x=0 is approximately equal to y[i] causing the code

else:
    if y[i] > y[i-1] and y[i] > y[i+1]:
        allMaximaPoints += [[x[i], y[i]], ]

to not account for that point. This is because when I changedx = np.linspace(-20.0, 20.0, num=int((20+20)/0.01)) to say x = np.linspace(-20.0, 20.0, num=int((20+20)/0.1)) i.e. larger steps in x, the local maxima at x=0 is correctly found. However, even if I changed the > signs in the above code to >=, that maxima at x=0 is still not counted.

Why is that so? How should I improve my code to get the correct results? Thanks!

The best library I found for maxima is peakutils: pypi.python.org/pypi/PeakUtils . It allows you to select sensitivity and threshold to make it robust to noise. — Gerges
– Gerges, Commented Feb 5, 2018 at 21:26

Alok Singhal · Accepted Answer · 2018-02-05 23:01:39Z

1

You might be better off using something like scipy.signal.find_peaks_cwt. Something like:

indices = scipy.signal.find_peaks_cwt(y, [1, 2, 3, 4], noise_perc=50)
plt.plot(x, y)
plt.plot(x[indices], y[indices], 'r.')
plt.show()

Results in:

answered Feb 5, 2018 at 23:01

Alok Singhal

96.8k21 gold badges130 silver badges158 bronze badges

Sign up to request clarification or add additional context in comments.

Comments

dtauxe · Accepted Answer · 2018-02-05 21:23:33Z

0

Short answer: x is never 0 (and you probably don't want it to be either).

We can test this by running

In [53]: 0 in x
Out[53]: False

This is actually fortunate as x=0 would cause a runtime warning and produce a "nan" (Not a Number) value due to the division by 0. The problem is that the function is undefined at 0.

Otherwise, you are essentially correct about why it doesn't work. If you swap the > for >= you'll see that the output is different:

In [55]: returnMaxima(2, x, y)
Out[55]:
[[-0.0050012503125778096, 0.99999166252624228],
[0.0050012503125778096, 0.99999166252624228]]

This is actually the "correct" output, as these are the largest values and are at the values of x closest to 0.

answered Feb 5, 2018 at 21:23

dtauxe

15311 bronze badges

Collectives™ on Stack Overflow

Improving algorithm to find local maxima of 1D numpy array

2 Answers 2

Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related