Python 中的峰值检测算法
Peak detection algorithm in Python
我正在 Python 中实施峰值检测算法,该算法仅检测那些高于阈值幅度的峰值。我不想使用内置函数,因为我还必须将此模拟扩展到硬件实现。
from math import sin,isnan
from pylab import *
def peakdet(v, delta,thresh,x):
delta=abs(delta)
maxtab = []
mintab = []
v = asarray(v)
mn, mx = v[0], v[0]
mnpos, mxpos = NaN, NaN
lookformax = True
for i in arange(len(v)):
this = v[i]
if abs(this)>thresh:
if this > mx:
mx = this
mxpos = x[i]
if this < mn:
mn = this
mnpos = x[i]
if lookformax:
if (this < mx-delta):
if (mx>abs(thresh)) and not isnan(mxpos):
maxtab.append((mxpos, mx))
mn = this
mnpos = x[i]
lookformax = False
else:
if (this > mn+delta):
if (mn<-abs(thresh)) and not isnan(mnpos):
mintab.append((mnpos, mn))
mx = this
mxpos = x[i]
lookformax = True
return array(maxtab), array(mintab)
#Input Signal
t=array(range(100))
series=0.3*sin(t)+0.7*cos(2*t)-0.5*sin(1.2*t)
thresh=0.95 #Threshold value
delta=0.0 #
a=zeros(len(t)) #
a[:]=thresh #
maxtab, mintab = peakdet(series,delta,thresh,t)
#Plotting output
scatter(array(maxtab)[:,0], array(maxtab)[:,1], color='red')
scatter(array(mintab)[:,0], array(mintab)[:,1], color='blue')
xlim([0,t[-1]])
title('Peak Detector')
grid(True)
plot(t,a,color='green',linestyle='--',dashes=(5,3))
plot(t,-a,color='green',linestyle='--',dashes=(5,3))
annotate('Threshold',xy=(t[-1],thresh),fontsize=9)
plot(t,series,'k')
show()
这个程序的问题是它无法检测到一些峰值,即使它们高于阈值。
这是我得到的输出:
我看到其他帖子有峰值检测问题,但找不到任何解决方案。请帮助并提出更正建议。
这些代码
if lookformax:
if (this < mx-delta):
if (mx>abs(thresh)) and not isnan(mxpos):
maxtab.append((mxpos, mx))
mn = this
mnpos = x[i]
lookformax = False
else:
if (this > mn+delta):
if (mn<-abs(thresh)) and not isnan(mnpos):
mintab.append((mnpos, mn))
mx = this
mxpos = x[i]
lookformax = True
仅运行条件下
if abs(this)>thresh:
所以你只能在阈值上方的下一个点小于它时找到一个峰值。
放出条件
您的函数使用了很多参数。您可以将问题分解为几个步骤:
- 首先检测所有高于阈值的点。将这些点添加到
maxthresh
和 minthresh
列表中。
- 遍历
maxthresh
列表,如果该点之前的y值小于该点,且该点之后的y值小于该点,则该点为峰。
- 遍历
minthresh
列表,如果该点之前的y值大于该点,且该点之后的y值大于该点,则该点为峰。
代码实现:
from math import sin
from matplotlib import pylab
from pylab import *
def peakdet(v, thresh):
maxthresh = []
minthresh = []
peaks = []
valleys = []
for x, y in v:
if y > thresh:
maxthresh.append((x, y))
elif y < -thresh:
minthresh.append((x, y))
for x, y in maxthresh:
try:
if (v[x - 1][1] < y) & (v[x + 1][1] < y):
peaks.append((x, y))
except Exception:
pass
for x, y in minthresh:
try:
if (v[x - 1][1] > y) & (v[x + 1][1] > y):
valleys.append((x, y))
except Exception:
pass
return peaks, valleys
测试代码:
# input signal
t = array(range(100))
series = 0.3 * sin(t) + 0.7 * cos(2 * t) - 0.5 * sin(1.2 * t)
arr = [*zip(t, series)] # create a list of tuples where the tuples represent the (x, y) values of the function
thresh = 0.95
peaks, valleys = peakdet(arr, thresh)
scatter([x for x, y in peaks], [y for x, y in peaks], color = 'red')
scatter([x for x, y in valleys], [y for x, y in valleys], color = 'blue')
plot(t, 100 * [thresh], color='green', linestyle='--', dashes=(5, 3))
plot(t, 100 * [-thresh], color='green', linestyle='--', dashes=(5, 3))
plot(t, series, 'k')
show()
额外测试以确保在多个点高于阈值时检测到峰值:
# input signal
t = array(range(100))
series = 6.3 * sin(t) + 4.7 * cos(2 * t) - 3.5 * sin(1.2 * t)
arr = [*zip(t, series)]
thresh = 0.95
peaks, valleys = peakdet(arr, thresh)
scatter([x for x, y in peaks], [y for x, y in peaks], color = 'red')
scatter([x for x, y in valleys], [y for x, y in valleys], color = 'blue')
plot(t, 100 * [thresh], color='green', linestyle='--', dashes=(5, 3))
plot(t, 100 * [-thresh], color='green', linestyle='--', dashes=(5, 3))
plot(t, series, 'k')
show()
所以,这里你有一个 numpythonic 解决方案(这比明确地做一个循环要好得多)。
我使用滚动函数将位置中的数字+1 或-1 移动。另外一个"peak"被定义为局部最大值,前后数都小于中心值。
完整代码为:
import numpy as np
import matplotlib.pyplot as plt
# input signal
x = np.arange(1,100,1)
y = 0.3 * np.sin(x) + 0.7 * np.cos(2 * x) - 0.5 * np.sin(1.2 * x)
threshold = 0.95
# max
maxi = np.where(np.where([(y - np.roll(y,1) > 0) & (y - np.roll(y,-1) > 0)],y, 0)> threshold, y,np.nan)
# min
mini = np.where(np.where([(y - np.roll(y,1) < 0) & (y - np.roll(y,-1) < 0)],y, 0)< -threshold, y,np.nan)
如果你绘制它,你会得到:
来自 scipy.signal
的 find_peaks
的解决方案
from scipy.signal import find_peaks
import numpy as np
import matplotlib.pyplot as plt
# Input signal
t = np.arange(100)
series = 0.3*np.sin(t)+0.7*np.cos(2*t)-0.5*np.sin(1.2*t)
# Threshold value (for height of peaks and valleys)
thresh = 0.95
# Find indices of peaks
peak_idx, _ = find_peaks(series, height=thresh)
# Find indices of valleys (from inverting the signal)
valley_idx, _ = find_peaks(-series, height=thresh)
# Plot signal
plt.plot(t, series)
# Plot threshold
plt.plot([min(t), max(t)], [thresh, thresh], '--')
plt.plot([min(t), max(t)], [-thresh, -thresh], '--')
# Plot peaks (red) and valleys (blue)
plt.plot(t[peak_idx], series[peak_idx], 'r.')
plt.plot(t[valley_idx], series[valley_idx], 'b.')
plt.show()
结果图如下所示。
注意find_peaks
有一个参数height
,也就是我们这里所说的thresh
。它还有一个名为 threshold
的参数,它正在做其他事情。
我正在 Python 中实施峰值检测算法,该算法仅检测那些高于阈值幅度的峰值。我不想使用内置函数,因为我还必须将此模拟扩展到硬件实现。
from math import sin,isnan
from pylab import *
def peakdet(v, delta,thresh,x):
delta=abs(delta)
maxtab = []
mintab = []
v = asarray(v)
mn, mx = v[0], v[0]
mnpos, mxpos = NaN, NaN
lookformax = True
for i in arange(len(v)):
this = v[i]
if abs(this)>thresh:
if this > mx:
mx = this
mxpos = x[i]
if this < mn:
mn = this
mnpos = x[i]
if lookformax:
if (this < mx-delta):
if (mx>abs(thresh)) and not isnan(mxpos):
maxtab.append((mxpos, mx))
mn = this
mnpos = x[i]
lookformax = False
else:
if (this > mn+delta):
if (mn<-abs(thresh)) and not isnan(mnpos):
mintab.append((mnpos, mn))
mx = this
mxpos = x[i]
lookformax = True
return array(maxtab), array(mintab)
#Input Signal
t=array(range(100))
series=0.3*sin(t)+0.7*cos(2*t)-0.5*sin(1.2*t)
thresh=0.95 #Threshold value
delta=0.0 #
a=zeros(len(t)) #
a[:]=thresh #
maxtab, mintab = peakdet(series,delta,thresh,t)
#Plotting output
scatter(array(maxtab)[:,0], array(maxtab)[:,1], color='red')
scatter(array(mintab)[:,0], array(mintab)[:,1], color='blue')
xlim([0,t[-1]])
title('Peak Detector')
grid(True)
plot(t,a,color='green',linestyle='--',dashes=(5,3))
plot(t,-a,color='green',linestyle='--',dashes=(5,3))
annotate('Threshold',xy=(t[-1],thresh),fontsize=9)
plot(t,series,'k')
show()
这个程序的问题是它无法检测到一些峰值,即使它们高于阈值。 这是我得到的输出:
我看到其他帖子有峰值检测问题,但找不到任何解决方案。请帮助并提出更正建议。
这些代码
if lookformax:
if (this < mx-delta):
if (mx>abs(thresh)) and not isnan(mxpos):
maxtab.append((mxpos, mx))
mn = this
mnpos = x[i]
lookformax = False
else:
if (this > mn+delta):
if (mn<-abs(thresh)) and not isnan(mnpos):
mintab.append((mnpos, mn))
mx = this
mxpos = x[i]
lookformax = True
仅运行条件下
if abs(this)>thresh:
所以你只能在阈值上方的下一个点小于它时找到一个峰值。
放出条件
您的函数使用了很多参数。您可以将问题分解为几个步骤:
- 首先检测所有高于阈值的点。将这些点添加到
maxthresh
和minthresh
列表中。 - 遍历
maxthresh
列表,如果该点之前的y值小于该点,且该点之后的y值小于该点,则该点为峰。 - 遍历
minthresh
列表,如果该点之前的y值大于该点,且该点之后的y值大于该点,则该点为峰。
代码实现:
from math import sin
from matplotlib import pylab
from pylab import *
def peakdet(v, thresh):
maxthresh = []
minthresh = []
peaks = []
valleys = []
for x, y in v:
if y > thresh:
maxthresh.append((x, y))
elif y < -thresh:
minthresh.append((x, y))
for x, y in maxthresh:
try:
if (v[x - 1][1] < y) & (v[x + 1][1] < y):
peaks.append((x, y))
except Exception:
pass
for x, y in minthresh:
try:
if (v[x - 1][1] > y) & (v[x + 1][1] > y):
valleys.append((x, y))
except Exception:
pass
return peaks, valleys
测试代码:
# input signal
t = array(range(100))
series = 0.3 * sin(t) + 0.7 * cos(2 * t) - 0.5 * sin(1.2 * t)
arr = [*zip(t, series)] # create a list of tuples where the tuples represent the (x, y) values of the function
thresh = 0.95
peaks, valleys = peakdet(arr, thresh)
scatter([x for x, y in peaks], [y for x, y in peaks], color = 'red')
scatter([x for x, y in valleys], [y for x, y in valleys], color = 'blue')
plot(t, 100 * [thresh], color='green', linestyle='--', dashes=(5, 3))
plot(t, 100 * [-thresh], color='green', linestyle='--', dashes=(5, 3))
plot(t, series, 'k')
show()
额外测试以确保在多个点高于阈值时检测到峰值:
# input signal
t = array(range(100))
series = 6.3 * sin(t) + 4.7 * cos(2 * t) - 3.5 * sin(1.2 * t)
arr = [*zip(t, series)]
thresh = 0.95
peaks, valleys = peakdet(arr, thresh)
scatter([x for x, y in peaks], [y for x, y in peaks], color = 'red')
scatter([x for x, y in valleys], [y for x, y in valleys], color = 'blue')
plot(t, 100 * [thresh], color='green', linestyle='--', dashes=(5, 3))
plot(t, 100 * [-thresh], color='green', linestyle='--', dashes=(5, 3))
plot(t, series, 'k')
show()
所以,这里你有一个 numpythonic 解决方案(这比明确地做一个循环要好得多)。
我使用滚动函数将位置中的数字+1 或-1 移动。另外一个"peak"被定义为局部最大值,前后数都小于中心值。
完整代码为:
import numpy as np
import matplotlib.pyplot as plt
# input signal
x = np.arange(1,100,1)
y = 0.3 * np.sin(x) + 0.7 * np.cos(2 * x) - 0.5 * np.sin(1.2 * x)
threshold = 0.95
# max
maxi = np.where(np.where([(y - np.roll(y,1) > 0) & (y - np.roll(y,-1) > 0)],y, 0)> threshold, y,np.nan)
# min
mini = np.where(np.where([(y - np.roll(y,1) < 0) & (y - np.roll(y,-1) < 0)],y, 0)< -threshold, y,np.nan)
如果你绘制它,你会得到:
来自 scipy.signal
的 find_peaks
的解决方案
from scipy.signal import find_peaks
import numpy as np
import matplotlib.pyplot as plt
# Input signal
t = np.arange(100)
series = 0.3*np.sin(t)+0.7*np.cos(2*t)-0.5*np.sin(1.2*t)
# Threshold value (for height of peaks and valleys)
thresh = 0.95
# Find indices of peaks
peak_idx, _ = find_peaks(series, height=thresh)
# Find indices of valleys (from inverting the signal)
valley_idx, _ = find_peaks(-series, height=thresh)
# Plot signal
plt.plot(t, series)
# Plot threshold
plt.plot([min(t), max(t)], [thresh, thresh], '--')
plt.plot([min(t), max(t)], [-thresh, -thresh], '--')
# Plot peaks (red) and valleys (blue)
plt.plot(t[peak_idx], series[peak_idx], 'r.')
plt.plot(t[valley_idx], series[valley_idx], 'b.')
plt.show()
结果图如下所示。
注意find_peaks
有一个参数height
,也就是我们这里所说的thresh
。它还有一个名为 threshold
的参数,它正在做其他事情。