peak_widths w.r.t to x axis

Question

I am using scipy.signal to calculate the width of the different peaks. I have 100 values wrt to different time points. I am using following code to calculate the peak, then width. The problem is it is not considering the time on x axis while calculating the width.

peaks_control, _ = find_peaks(x_control, height=2100)
time_control     = time[:100]
width_control    = peak_widths(x_control, peaks_control, rel_height=0.9)

The output of width_control is

array([12.84785714, 13.21299534, 13.4502381 , 12.71311143]),
array([2042.5, 2048.8, 2057.4, 2065. ]),
array([ 5.795 ,28.29469697, 51.245 , 74.17150396]),
array([18.64285714, 41.50769231, 64.6952381 , 86.88461538]))

I am using following to use time on x axis and show the signals, which is correct

plt.plot(time_control, x_control)
plt.plot(time_control[peaks_control], x_control[peaks_control], "x")
#plt.plot(np.zeros_like(x_control), "--", color="gray")
#plt.xlim(time_control.tolist())
plt.title('Control')
plt.xlabel('Time (s)')
plt.ylabel('RFU')
plt.show()

I am using following code to show the width also, but not able to put the actual time on x axis.

plt.plot(x_control)
plt.plot(peaks_control, x_control[peaks_control], "x")
plt.hlines(*width_control[1:], color="C3")

plt.title('Control')
plt.xlabel('Time (s)')
plt.ylabel('RFU')
plt.show()

Answer 1

I had the same problem just now, so here's my solution (there are probably more elegant solutions but this worked for me):

peak_widths() returns the widths (in samples), height at which the widths were calculated, and the interpolated positions of left and right intersection points of a horizontal line at the respective evaluation height (also in samples).

To convert those values from samples back to our x-axis we can use scipy.interpolate.interp1():

from scipy.signal import find_peaks, peak_widths
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
import numpy as np

def index_to_xdata(xdata, indices):
    "interpolate the values from signal.peak_widths to xdata"
    ind = np.arange(len(xdata))
    f = interp1d(ind,xdata)
    return f(indices)

x = np.linspace(0, 1, 10)
y = np.sin(4*x-0.2)

peaks, _ = find_peaks(y)
widths, width_heights, left_ips, right_ips = peak_widths(y, peaks)

widths = index_to_xdata(x, widths)
left_ips = index_to_xdata(x, left_ips)
right_ips = index_to_xdata(x, right_ips)

plt.plot(x,y)
plt.plot(x[peaks], y[peaks], "x")
plt.hlines(width_heights, left_ips, right_ips, color='r')
plt.xlabel('x values')
plt.ylabel('y values')
plt.show()

image of plot