python - 改进算法以找到一维 numpy 数组的局部最大值

Question

我试图找到函数 f(x) = (sin(x)/x)^2 的局部最大值。对于近似解，我初始化了两个变量 x 和 y 并首先绘制了一个图形以进行可视化表示。

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

这显示 graph .

然后我尝试创建一个算法来找到最大值:

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

当我尝试 returnMaxima(2, x, y) 时，我得到了 [[-4.4961240310077528, 0.04719010162459622],
[4.4961240310077528, 0.04719010162459622]]。

这是不正确的，因为它跳过了 x = 0 处的局部最大值。我怀疑这是因为 y[i-1] 和 y[i+1] 值与 y[i] 处的最大值相邻x=0 大约等于 y[i] 导致代码

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

不考虑这一点。这是因为当我将 x = np.linspace(-20.0, 20.0, num=int((20+20)/0.01)) 更改为 x = np.linspace(-20.0 , 20.0, num=int((20+20)/0.1)) 即 x 中的较大步长，正确找到了 x=0 处的局部最大值。然而，即使我将上面代码中的 > 符号更改为 >=，x=0 处的最大值仍未计算在内。

为什么会这样？我应该如何改进我的代码以获得正确的结果？谢谢!

Answer 1

你最好使用类似 scipy.signal.find_peaks_cwt 的东西.像这样的东西:

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

结果: