python scipy.signal.peak_widths -->绝对高度？(fft -3dB阻尼)
python scipy.signal.peak_widths -->绝对高度？(fft -3dB阻尼)
提问于 2018-12-14 19:15:54
https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.peak_widths.html
我认为关联函数只能计算出相对高度处的峰值宽度。有没有人知道是否有一个函数可以计算所有峰值的宽度(peak_amplitude - x)？
目前我正在尝试修改原来的内部函数"_peak_widths“。使用cimport时已失败。仅部分理解此处的源代码。我在代码中添加了我要做修改的地方。
with nogil:
for p in range(peaks.shape[0]):
i_min = left_bases[p]
i_max = right_bases[p]
peak = peaks[p]
# Validate bounds and order
if not 0 <= i_min <= peak <= i_max < x.shape[0]:
with gil:
raise ValueError("prominence data is invalid for peak {}"
.format(peak))
height = width_heights[p] = x[peak] - prominences[p] * rel_height 在此处更改为xpeak 3
# Find intersection point on left side
i = peak
while i_min < i and height < x[i]:
i -= 1
left_ip = <np.float64_t>i
if x[i] < height:
# Interpolate if true intersection height is between samples
left_ip += (height - x[i]) / (x[i + 1] - x[i])
# Find intersection point on right side
i = peak
while i < i_max and height < x[i]:
i += 1
right_ip = <np.float64_t>i
if x[i] < height:
# Interpolate if true intersection height is between samples
right_ip -= (height - x[i]) / (x[i - 1] - x[i])
widths[p] = right_ip - left_ip
if widths[p] == 0:
show_warning = True
left_ips[p] = left_ip
right_ips[p] = right_ip回答 2
Stack Overflow用户
回答已采纳
发布于 2019-09-12 05:06:51
如果这仍然与您相关，您可以通过传入修改后的prominence_data，“按原样”使用scipy.signal.peak_widths来实现您想要的结果。基于您自己的answer
import numpy as np
from scipy.signal import find_peaks, peak_prominences, peak_widths
# Create sample data
x = np.linspace(0, 6 * np.pi, 1000)
x = np.sin(x) + 0.6 * np.sin(2.6 * x)
# Find peaks
peaks, _ = find_peaks(x)
prominences, left_bases, right_bases = peak_prominences(x, peaks)正如peak_widths文档中所述，测量宽度的高度是以h_eval = h_peak - prominence * relative_height计算的
我们可以通过参数prominence_data和rel_height控制后两个变量。因此，我们可以创建一个所有值都相同的数组，并使用它来创建绝对高度，而不是传递计算出的每个峰值不同的prominence：
# Create constant offset as a replacement for prominences
offset = np.ones_like(prominences)
# Calculate widths at x[peaks] - offset * rel_height
widths, h_eval, left_ips, right_ips = peak_widths(
x, peaks,
rel_height=1,
prominence_data=(offset, left_bases, right_bases)
)
# Check that h_eval is 1 everywhere
np.testing.assert_equal(x[peaks] - h_eval, 1)
# Visualize result
import matplotlib.pyplot as plt
plt.plot(x)
plt.plot(peaks, x[peaks], "x")
plt.hlines(h_eval, left_ips, right_ips, color="C2")
plt.show()
正如您所看到的，在相同的常量偏移量1处评估每个峰值的宽度。通过使用peak_prominences提供的原始left_bases和right_bases，我们限制了最大测量宽度(例如，参见299和533处的峰值)。如果想要消除该限制，则必须自己创建这些数组。
Stack Overflow用户
发布于 2018-12-14 21:36:26
我刚删除了c的内容。这就是我的解决方案：
def gauss(x, p): # p[0]==mean, p[1]==stdev
return 1.0/(p[1]*np.sqrt(2*np.pi))*np.exp(-(x-p[0])**2/(2*p[1]**2))
def _peak_widths(x,peaks,prop,val=3):
i_min = prop['left_bases']
i_max = prop['right_bases']
peak = peaks[0]
# Validate bounds and order
height = x[peak] - val
# Find intersection point on left side
i = peak
while i_min < i and height < x[i]:
i -= 1
left_ip = i
if x[i] < height:
# Interpolate if true intersection height is between samples
left_ip += (height - x[i]) / (x[i + 1] - x[i])
# Find intersection point on right side
i = peak
while i < i_max and height < x[i]:
i += 1
right_ip = i
if x[i] < height:
# Interpolate if true intersection height is between samples
right_ip -= (height - x[i]) / (x[i - 1] - x[i])
widths = right_ip - left_ip
left_ips = left_ip
right_ips = right_ip
return [height, widths, int(left_ips), int(right_ips)]
if __name__ == '__main__':
# Create some sample data
known_param = np.array([2.0, 0.07])
xmin,xmax = -1.0, 5.0
N = 1000
X = np.linspace(xmin,xmax,N)
Y = gauss(X, known_param)
fig, ax= plt.subplots()
ax.plot(X,Y)
#find peaks
peaks, prop = signal.find_peaks(Y, prominence = 3.1)
ax.scatter(X[peaks],Y[peaks], color='r')
#calculate peak width
y, widths, x1, x2 = _peak_widths(Y,peaks, prop)
print(f'width = { X[x1] - X[x2]}')
l = mlines.Line2D([X[x1],X[x2]], [y,y], color='r')
ax.add_line(l)
plt.show()页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/53778703
复制相关文章
相似问题
腾讯云开发者
