scipy 如何计算下降到0的峰?Python查找峰

sd2nnvve  于 2022-11-10  发布在  Python
关注(0)|答案(1)|浏览(204)

我使用Scipy的find_peaks来计算时间序列中的峰值数量。
我需要计算峰的数量,要求它从0开始,然后福尔斯到0。从右边开始的第二个峰(用一条垂直线表示)在这里被计算,但它不应该被计算,因为它在最后一个峰之前没有下降到0。有没有办法在find_peaks中指定它?

peaks1 = find_peaks(array, height=(1,1.5),prominence=1)
peaks1_5 = find_peaks(array, height=(1.5,2),prominence=1.5)
peaks2 = find_peaks(array, height=2,prominence=2)

fig, ax = plt.subplots(figsize=(30, 10), dpi=80)
plt.plot(spi_neg['date'],spi["SPI-12"])
[plt.axvline(spi_neg.date.iloc[p],c='red',linewidth=0.3) for p in peaks1[0]]
[plt.axvline(spi_neg.date.iloc[p],c='green',linewidth=0.3) for p in peaks1_5[0]]
[plt.axvline(spi_neg.date.iloc[p],c='purple',linewidth=0.3) for p in peaks2[0]]
plt.axhline(2,linestyle='dashed',linewidth=1)
plt.axhline(1.5,linestyle='dashed',linewidth=1)
plt.axhline(1,linestyle='dashed',linewidth=1)

Peaks chart

eqfvzcg8

eqfvzcg81#

一个带问题的运行代码将是有帮助的,并且一个更精确的可数峰的定义也是有帮助的;- )
首先,我们生成一些数据

import numpy as np
import matplotlib.pyplot as plt

# ---- generate data

mp = 200
freq = 20
t = np.linspace(0,freq*np.pi,mp)
signal = np.sin(t)
noise = np.random.rand(mp)
X = 0.5*signal + noise

# ---- scale X

def scale01(a):
    return (a-a.min())/(a.max()-a.min())

X = scale01(X) - 0.5
X = np.maximum(X,0.0)

# ---- grafics

with plt.style.context('ggplot'):       
    fig = plt.figure(figsize=(15,3))
    plt.plot(t, X)
    plt.plot(t, X, 'o')

现在我们识别零湖和非零岛

a = np.array(np.where(X<=0))[0]   # extract the indices with X<=0
b = np.array(np.where(X>0) )[0]   # extract the indices with X>0

with plt.style.context('ggplot'):       
    fig = plt.figure(figsize=(15,3))
    plt.plot(t[b], X[b], 'or', label=">0")
    plt.vlines(t[b], 0, X[b], colors='k')
    plt.plot(t[a], -X[a], 'og', label="<=0")
    plt.legend(); plt.show()

接下来我们用numpy数组填充列表中的非零岛。每个numpy数组包含一个非零岛。

X_ = X[b]
m = len(X_)
list_y = list()
list_Y = list()

for j in range(1,m):
    if  b[j]-b[j-1]>1 :
        list_Y.append(list_y)
        list_y = list()  
        #print("------------------------------------------------------ new list")
    #print(j, b[j], X_[j])
    list_y.append(X_[j])
list_Y.append(list_y)

print("list_Y"); 
n = len(list_Y)
for j in range(n):
    print(list_Y[j])

对于列表中的每个numpy数组,您可以根据您的定义(我无法完全捕获)评估峰

list_Y
[0.22062475371241386, 0.29207471279008657, 0.35072832015294586, 0.1251594602284437, 0.24379282278250836, 0.06896727821692716]
[0.06271739133976328]
[0.2689504650818615, 0.011887999386713255, 0.055442917743508624, 0.2876317343278316, 0.24084993011027578, 0.12097014134978235]
[0.1907699022464584]
[0.08249052680941726]
[0.10205561805376617]
[0.18903867830269638, 0.26990334850384257, 0.5, 0.3288200602696131, 0.05036869827824486, 0.040381419904307436]
[0.08618838642790339]
[0.0053279353208096625, 0.3468189863146819, 0.05644254569326557, 0.3985674171686334, 0.14897985190026097, 0.0025548308606182513, 0.32765453143333545, 0.3328107320769136, 0.1838328219774621, 0.21123652127176762]
[0.18870251894797663]
[0.13453490446867422, 0.25258744200608363, 0.4981866504733391, 0.35180043079867795, 0.08425183513691303, 0.3376976620831299, 0.22348609066402825]
[0.0716155758184146]
[0.052227024152749935, 0.08639499278421903]
[0.1581304564482665, 0.2273016493144655, 0.26721741895716056, 0.33665669827299305, 0.19255497112246478, 0.16227876457894175]
[0.10236622631923908, 0.06039140456773806, 0.053391261130168344]
[0.21170561257978093, 0.11669466945342633, 0.2479665749659119, 0.25792206298341824, 0.19579440295962314, 0.15210847528158666, 0.23531008247873408]
[0.05340116678342899]
[0.2088166123161308, 0.26031072203571415, 0.2786317264092839, 0.289871721166855, 0.25460661866030165, 0.3214937091565473, 0.36293451974436275]
[0.04525610202919361, 0.1740374143631349, 0.17258947174880612]
[0.14217066607610684, 0.03435965315335088, 0.09996473411377804, 0.48290831305140514, 0.09407783896892297]
[0.03224632110920911, 0.08787466747977346, 0.20032938280871493, 0.23646809723694695, 0.13380244841935984, 0.05305696510866664, 0.2657761536751757, 0.34514204517200975]
[0.17123014194168007, 0.2397521290598289]

相关问题