python scipy.signal.peak_widths -->绝对高度?(fft -3dB 阻尼) [英] python scipy.signal.peak_widths --> absolute heigth? (fft -3dB damping)

问题描述

如您所见，每个峰的宽度以相同的恒定偏移量 1 计算.通过使用 peak_prominences 提供的原始 left_bases 和 right_bases 我们正在限制最大测量宽度(例如，见 299 和 533 处的峰值).如果您想取消该限制，您必须自己创建这些数组.

https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.peak_widths.html

I think the linked function can only calculate the peak widths at a relative height. Does anyone know if there is a function that calculates the width at a fixed value (peak_amplitude - x) for all peaks?

Currently I am trying to change the original inner function "_peak_widths". Fail already with the cimport. Understand the source code here only partially. I added in the code where I would make a modification.

 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 

CHANGE HERE TO x[peak] - 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

解决方案

In case this is still relevant to you, you can use scipy.signal.peak_widths "as is" to achieve what you want by passing in modified prominence_data. Based on your own 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)

As stated in peak_widths's documentation the height at which the width is measured is calculated as h_eval = h_peak - prominence * relative_height

We can control the latter two variables through the parameters prominence_data and rel_height. So instead of passing in the calculated prominence which differs for each peak we can create an array where all values are the same and use that to create an absolute height:

# 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()

As you can see the width is evaluated for each peak at the same constant offset of 1. By using the original left_bases and right_bases as provided by peak_prominences we are limiting the maximal measured width (e.g. see peaks at 299 and 533). If you want to remove that limitation you must create these arrays yourself.

这篇关于python scipy.signal.peak_widths -->绝对高度?(fft -3dB 阻尼)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！