如何平滑图形以找到所有局部最大值？

Question

我想找到一组给定点的局部最大值。 对于前。 -

xpoints_ball = np.array([2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,189,190,191,192,193,194,195,197,198,199,200,201,202,203,204,205,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,224,234,235,236,237,238,239,240,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298])
ypoints_ball = np.array([238,241,243,243,246,251,254,258,259,263,270,276,280,287,288,293,301,312,319,330,331,339,351,360,373,386,385,398,407,428,394,453,454,469,486,504,520,537,538,555,575,593,615,634,633,655,676,697,724,733,734,732,732,730,729,729,728,729,730,730,733,736,736,739,742,746,751,756,756,760,767,774,781,788,790,796,807,773,742,715,718,683,660,636,615,594,595,572,556,538,524,510,511,494,469,455,448,448,436,429,418,411,403,403,394,389,385,381,375,376,372,369,364,363,360,360,357,357,356,355,354,352,340,330,318,310,311,301,292,285,276,268,304,264,256,251,246,241,240,235,231,226,220,227,225,234,237,229,240,241,244,249,250,251,255,257,262,264,268,268,274,277,280,286,291,291,294,301,309,313,318,319,325,331,338,354,354,361,367,379,386,393,403,411,421,430,441,441,454,461,474,494,496,509,522,534,545,561,560,573,587,603,616,632,629,643,662,688,721,721,719,720,718,718,722,721,723,723,727,729,732,735,739,739,742,747,752,729,692,692,667,639,611,591,565,569,547,529,512,491,478,463,449,439,429,415,418,409,400,391,384,380,381,374,367,363,361,360,360,358,354,354,352,352,352,352,352,352,347,334,338])

plt.plot(xpoints_ball, ypoints_ball)
plt.show()

为了平滑图形，我使用了三次插值 -

f_ball = interp1d(xpoints_ball, ypoints_ball, kind='cubic' , fill_value="extrapolate")
plt.plot(xpoints_ball, f_ball(xpoints_ball))
plt.show()

现在，为了找到局部最大值，我使用 Scipy 库的 peak 函数 -

coordinates = f_ball(xpoints_ball)
peaks, _ = find_peaks(coordinates)
plt.plot(coordinates)
plt.plot(peaks, coordinates[peaks], "x")

这个结果给了我很多不是局部最大值的误报值。 如何修改上面的代码，使局部最大值的结果是这样的 -

Answer 1

您可以使用例如平滑您的数据高斯滤波并使用sigma参数对其进行控制

from scipy.ndimage import gaussian_filter
from scipy.signal import find_peaks

xpoints_ball = np.array([2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,189,190,191,192,193,194,195,197,198,199,200,201,202,203,204,205,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,224,234,235,236,237,238,239,240,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298])
ypoints_ball = np.array([238,241,243,243,246,251,254,258,259,263,270,276,280,287,288,293,301,312,319,330,331,339,351,360,373,386,385,398,407,428,394,453,454,469,486,504,520,537,538,555,575,593,615,634,633,655,676,697,724,733,734,732,732,730,729,729,728,729,730,730,733,736,736,739,742,746,751,756,756,760,767,774,781,788,790,796,807,773,742,715,718,683,660,636,615,594,595,572,556,538,524,510,511,494,469,455,448,448,436,429,418,411,403,403,394,389,385,381,375,376,372,369,364,363,360,360,357,357,356,355,354,352,340,330,318,310,311,301,292,285,276,268,304,264,256,251,246,241,240,235,231,226,220,227,225,234,237,229,240,241,244,249,250,251,255,257,262,264,268,268,274,277,280,286,291,291,294,301,309,313,318,319,325,331,338,354,354,361,367,379,386,393,403,411,421,430,441,441,454,461,474,494,496,509,522,534,545,561,560,573,587,603,616,632,629,643,662,688,721,721,719,720,718,718,722,721,723,723,727,729,732,735,739,739,742,747,752,729,692,692,667,639,611,591,565,569,547,529,512,491,478,463,449,439,429,415,418,409,400,391,384,380,381,374,367,363,361,360,360,358,354,354,352,352,352,352,352,352,347,334,338])


ypoints_ball_smooth=gaussian_filter(ypoints_ball, sigma=0.8)

plt.plot(xpoints_ball, ypoints_ball)
plt.plot(xpoints_ball, ypoints_ball_smooth,"--")

peaks, _ = find_peaks(ypoints_ball_smooth)
plot(xpoints_ball[peaks],ypoints_ball_smooth[peaks],"ro")

Answer 2

试试这个功能

import numpy as np

def smooth(array, window_len=10, window="hanning"):
    """np.array, int, str => np.array
    Smooth the data using a window with requested size.
    
    This method is based on the convolution of a scaled window with the signal.
    The signal is prepared by introducing reflected copies of the signal 
    (with the window size) at both ends so that transient parts are minimized
    in the beginning and end part of the output signal.
    
    input:
        x: the input signal 
        window_len: the dimension of the smoothing window; should be an odd integer
        window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
            flat window will produce a moving average smoothing.
    output:
        the smoothed signal
        
    example:
        t = linspace(-2,2,0.1)
        x = sin(t)+randn(len(t))*0.1
        y = smooth(x)
    
    see also: 
    
    numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
    scipy.signal.lfilter
    """

    if array.ndim != 1:
        raise ValueError("smooth only accepts 1 dimension arrays.")

    if array.size < window_len:
        raise ValueError("Input vector needs to be bigger than window size.")

    if window_len < 3:
        return x

    if not window in ["flat", "hanning", "hamming", "bartlett", "blackman"]:
        raise ValueError(
            "Window must be one of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
        )

    s = np.r_[array[window_len - 1 : 0 : -1], array, array[-2 : -window_len - 1 : -1]]
    
    if window == "flat":  # moving average
        w = np.ones(window_len, "d")
    else:
        w = eval("np." + window + "(window_len)")

    y = np.convolve(w / w.sum(), s, mode="valid")

    y = y[int(window_len / 2 - 1) : -int(window_len / 2)]
    offset = len(y) - len(array)  # in case input and output are not of the same lenght
    assert len(array) == len(y[offset:])
    return y[offset:]