4

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

PS 我无法上传错误的峰值检测波形图的图像。我希望我的解释足够充分......代码如下

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
4

1 回答 1

3

将凹槽与真实峰区分开来的能力意味着您在峰之间具有基本间距。换句话说，有一个您希望运行峰值检测搜索的最小频率分辨率。如果您放大一个比本底噪声窄得多的信号，您会观察到似乎每隔几个样本就会出现“峰值”的锯齿形。

听起来你想做的是以下内容：

  1. 平滑信号。
  2. 找到“真正的”峰。

或者更准确地说，

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

第 1 步：低通滤波

做第 1 步，我推荐你使用scipy 提供的信号处理工具。

我改编了这个食谱条目，它展示了如何使用 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 步：寻找峰值

对于第 2 步，我建议您使用scipy 提供的小波变换峰值查找器。您必须提供过滤后的信号和从最小到最大可能峰宽运行的向量作为输入。该向量将用作小波变换的基础。

#------------------------------------------------
# Step 2: Find the peaks
#------------------------------------------------

figure(4)
plot(t[N-1:]-delay, filtered_x[N-1:], 'g', linewidth=1)

peakind = find_peaks_cwt(filtered_x, arange(3,20))
scatter([t[i] - delay for i in peakind], [filtered_x[i] for i in peakind], color="red")
for i in peakind:
    print t[i] + delay

xlabel('t')
grid(True)

show()

在此处输入图像描述

于 2013-10-14T03:45:56.087 回答