Getting the width and area of peaks from Scipy Signal object

Question

How do I get the peak objects with the properties such as position, peak aarea, peak width etc from Scipy Signal function using cwt get peaks method:

def CWT(trace):
x = []
y = []
for i in range(len(trace)):
    x.append(trace[i].Position)
    y.append(trace[i].Intensity)
x = np.asarray(x)
y = np.asarray(y)
return signal.find_peaks_cwt(x,y)

This just returns an array?

Answer 1

First, it looks like you are using find_peaks_cwt incorrectly. Its two positional parameters are not x and y coordinates of data points. The first parameter is y-values. The x-values are not taken at all, they are assumed to be 0,1,2,.... The second parameter is a list of peak widths that you are interested in;

1-D array of widths to use for calculating the CWT matrix. In general, this range should cover the expected width of peaks of interest.

There is no reason for width parameter to be of the same size as the data array. In my example below, the data has 500 values, but the widths I use are 30...99.

Second, this method only finds the position of peaks (the array you get has the indexes of peaks). There is no analysis of their widths and areas. You will either have to look elsewhere (blog post Peak Detection in the Python World lists some alternatives, though none of them return the data you want), or come up with your own method of estimating those things.

My attempt is below. It does the following:

1. Cuts the signal by midpoints between peaks
2. For each piece, uses the median of values in it as the baseline
3. Declares the peak to consist of all values that are greater than 0.5*(peak value + baseline), i.e., midway between median and maximum.
4. Finds where the peak begins and where it ends. (The width is just the difference of these)
5. Declares the area of the peak to be the sum of (y - baseline) over the interval found in step 4.

Complete example:

t = np.linspace(0, 4.2, 500)
y = np.sin(t**2) + np.random.normal(0, 0.03, size=t.shape)  # simulated noisy signal
peaks = find_peaks_cwt(y, np.arange(30, 100, 10))

cuts = (peaks[1:] + peaks[:-1])//2    # where to cut the signal 
cuts = np.insert(cuts, [0, cuts.size], [0, t.size])
peak_begins = np.zeros_like(peaks)
peak_ends = np.zeros_like(peaks)
areas = np.zeros(peaks.shape)
for i in range(peaks.size):
  peak_value = y[peaks[i]]
  y_cut = y[cuts[i]:cuts[i+1]]           # piece of signal with 1 peak
  baseline = np.median(y_cut)
  large = np.where(y_cut > 0.5*(peak_value + baseline))[0]
  peak_begins[i] = large.min() + cuts[i]
  peak_ends[i] = large.max() + cuts[i]
  areas[i] = np.sum(y[peak_begins[i]:peak_ends[i]] - baseline)

The arrays areas, peak_begins and peak_ends are of interest here. The widths are [84 47 36], indicating the peaks get thinner (recall these are in index units, the width is the number of data points in the peak). I use this data to color the peaks in red:

widths = peak_ends - peak_begins
print(widths, areas)
plt.plot(t, y)
for i in range(peaks.size):
  plt.plot(t[peak_begins[i]:peak_ends[i]], y[peak_begins[i]:peak_ends[i]], 'r')
plt.show()