python scipy.signal.peak_widths --> absolute heigth? (fft -3dB damping)

Question

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

lagru · Accepted Answer · 2019-09-12T19:11:17.297

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.

score 0 · Answer 2 · answered Dec 14 '18 at 13:36

I just removed the c content. Thats my solution:

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

python scipy.signal.peak_widths --> absolute heigth? (fft -3dB damping)

CHANGE HERE TO x[peak] - 3

2 Answers2

Linked