CODEWITHSUNDEEP
pythonsignalssignal-processingdata-processing
Alex Walczak
Alex Walczak

Reputation: 1396

In a dataset with multiple peaks, how do I find FWHM of each peak (Python)?

Here's my denoised data:enter image description here

I've calculated my peaks, but if I define the peak width as full width at half maximum (FWHM) (while assuming zero is defined as the smallest point in data between ~25 to ~375 on x axis), is there a Numpy/Scipy way to calculate their width? I'm okay with coding my own function, but I'm not sure how to begin an implementation since there are multiple peaks in my data.

Upvotes: 3

Views: 1461

Answers (1)

Vlad
Vlad

Reputation: 1

Here's my solution:

'''Calculates fwhm of each peak in signal'''

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

# ==== import data ==== #
x, y = np.genfromtxt('data.txt', skip_header=True, unpack=True)

# ==== fwhm calc ==== #
peaks, _ = find_peaks(y, distance=20) # peaks are separated by at least 20px
params = dict(rel_height=0.5, wlen=35)
'''wlen parameter is recommended to use for periodic signals; 
it should be slighly bigger than distance between peaks
its meaning is area around peak's vertex where algorithm searches for peaks' bases'''
fwhm_0 = peak_widths(y, peaks, **params) # output columns are widths, width_heights, left_ips, right_ips
'''** operator to unpack the dictionary'''
fwhm_widths = fwhm_0[0]

# ==== plot ==== #
plt.plot(y)
plt.plot(peaks, y[peaks], "x")
plt.hlines(*fwhm_0[1:], color="C2")
'''An asterisk means unpacking. It spreads the iterable into several arguments'''
plt.show()

# ==== remove the outliers ==== #
zs = scipy.stats.zscore(fwhm_widths)
fwhm = fwhm_widths[np.abs(zs) < 3]

# ==== output ==== #
M = 1 # scale factor, units/pixels
Delta = M * statistics.mean(np.diff(peaks)).item()
'''.item() converts from numpy int64 to simple float '''
print("Mean distance between peaks is :", Delta, "units")
fwhm = M * fwhm
fwhm_mean = statistics.mean(fwhm)
print("Mean FWHM is :", fwhm_mean, "units")
fwhm_stdev = statistics.stdev(fwhm)
print("Standart deviation is :", fwhm_stdev, "units")

Upvotes: 0

Related Questions