改进的波检测算法

Question

我是python和编程的新手，在我的样条插值图上研究波峰检测算法。我使用了这个链接上给出的代码：https://gist.github.com/endolith/250860。我必须使该算法适用于任何类型的波形，即低振幅和高振幅，基线不对齐等。目标是计算曲线图中的波数。但是我的峰值检测计算出“无效”的峰值，所以给出了错误的答案。所谓的“无效”峰值，我的意思是，如果在波峰处有两个彼此接近的缺口，程序检测到2个峰值，即2个波，而实际上它只有1个波。我试过改变在链接上给出的峰值检测函数中定义的“delta”参数，但这并不能解决我正在研究的泛化目标。请对算法或我应该使用的任何其他方法提出任何改进建议。欢迎任何帮助。提前谢谢。在

另外，我无法上传错误峰值检测波图的图像。我希望我的解释足够。。。 代码如下

wave = f1(xnew)/(max(f1(xnew))) ##interpolated wave
maxtab1, mintab1 = peakdet(wave,.005) 
maxi = max(maxtab1[:,1])
for i in range(len(maxtab1)):
    if maxtab1[i,1] > (.55 * maxi) : ## Thresholding
        maxtab.append((maxtab1[i,0],maxtab1[i,1]))
arr_maxtab = np.array(maxtab)
dist = 1500 ## Threshold defined for minimum distance between two notches to be considered as two separate waves
mtdiff = diff(arr_maxtabrr[:,0])
final_maxtab = []
new_arr = []
for i in range(len(mtdiff)):
if mtdiff[i] < dist :
            new_arr.append((arr_maxtab[i+1,0],arr_maxtab[i+1,1])) 
for i in range(len(arr_maxtab)):
if  not arr_maxtab[i,0] in new_arr[:,0]:
    final_maxtab.append((arr_maxtab[i,0], arr_maxtab[i,1]))
else:
final_maxtab = maxtab

Answer 1

区分缺口和真实峰值的能力意味着峰值之间有一个基本间距。换言之，有一个最低频率分辨率，你想运行你的峰值检测搜索。如果你放大一个比噪声底限窄得多的信号，你会观察到似乎每隔几个采样点就会出现“峰值”的锯齿形。在

听起来你想做的是：

1. 把信号调平。在
2. 找到真正的山峰。在

或者更准确地说

1. 使信号通过低通滤波器。在
2. 在你可接受的峰值宽度范围内找到具有足够信噪比的峰值。在

第一步：低通滤波

要执行步骤1，我建议您使用scipy提供的signal processing tools。在

我改编了this cookbook entry，它展示了如何使用FIR滤波器来使用scipy对信号进行低通。在

from numpy import cos, sin, pi, absolute, arange, arange
from scipy.signal import kaiserord, lfilter, firwin, freqz, find_peaks_cwt
from pylab import figure, clf, plot, xlabel, ylabel, xlim, ylim, title, grid, axes, show, scatter

#                        
# Create a signal. Using wave for the signal name.
#                        

sample_rate = 100.0
nsamples = 400
t = arange(nsamples) / sample_rate
wave = cos(2*pi*0.5*t) + 0.2*sin(2*pi*2.5*t+0.1) + \
        0.2*sin(2*pi*15.3*t) + 0.1*sin(2*pi*16.7*t + 0.1) + \
            0.1*sin(2*pi*23.45*t+.8)


#                        
# Create a FIR filter and apply it to wave.
#                        

# The Nyquist rate of the signal.
nyq_rate = sample_rate / 2.0

# The desired width of the transition from pass to stop,
# relative to the Nyquist rate.  We'll design the filter
# with a 5 Hz transition width.
width = 5.0/nyq_rate

# The desired attenuation in the stop band, in dB.
ripple_db = 60.0

# Compute the order and Kaiser parameter for the FIR filter.
N, beta = kaiserord(ripple_db, width)

# The cutoff frequency of the filter.
cutoff_hz = 10.0

# Use firwin with a Kaiser window to create a lowpass FIR filter.
taps = firwin(N, cutoff_hz/nyq_rate, window=('kaiser', beta))

# Use lfilter to filter x with the FIR filter.
filtered_x = lfilter(taps, 1.0, wave)

第二步：寻峰

对于步骤2，我建议您使用scipy提供的wavelet transformation peak finder。您必须提供滤波信号和从最小峰值宽度到最大峰值宽度的矢量作为输入。该矢量将作为小波变换的基础。在

^{pr2}$