I have been trying to figure out the full width half maximum (FWHM) of the the blue peak (see image). The green peak and the magenta peak combined make up the blue peak. I have been using the following equation to find the FWHM of the green and magenta peaks:
fwhm = 2*np.sqrt(2*(math.log(2)))*sd where sd = standard deviation. I created the green and magenta peaks and I know the standard deviation which is why I can use that equation.
I created the green and magenta peaks using the following code:
def make_norm_dist(self, x, mean, sd): import numpy as np norm =  for i in range(x.size): norm += [1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x[i] - mean)**2/(2*sd**2))] return np.array(norm)
If I did not know the blue peak was made up of two peaks and I only had the blue peak in my data, how would I find the FWHM?
I have been using this code to find the peak top:
peak_top = 0.0e-1000 for i in x_axis: if i > peak_top: peak_top = i
I could divide the
peak_top by 2 to find the half height and then try and find y-values corresponding to the half height, but then I would run into trouble if there are no x-values exactly matching the half height.
I am pretty sure there is a more elegant solution to the one I am trying.
You can use spline to fit the [blue curve - peak/2], and then find it's roots:
import numpy as np from scipy.interpolate import UnivariateSpline def make_norm_dist(x, mean, sd): return 1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x - mean)**2/(2*sd**2)) x = np.linspace(10, 110, 1000) green = make_norm_dist(x, 50, 10) pink = make_norm_dist(x, 60, 10) blue = green + pink # create a spline of x and blue-np.max(blue)/2 spline = UnivariateSpline(x, blue-np.max(blue)/2, s=0) r1, r2 = spline.roots() # find the roots import pylab as pl pl.plot(x, blue) pl.axvspan(r1, r2, facecolor='g', alpha=0.5) pl.show()
Here is the result:
This worked for me in iPython (quick and dirty, can be reduced to 3 lines):
def FWHM(X,Y): half_max = max(Y) / 2. #find when function crosses line half_max (when sign of diff flips) #take the 'derivative' of signum(half_max - Y) d = sign(half_max - array(Y[0:-1])) - sign(half_max - array(Y[1:])) #plot(X[0:len(d)],d) #if you are interested #find the left and right most indexes left_idx = find(d > 0) right_idx = find(d < 0)[-1] return X[right_idx] - X[left_idx] #return the difference (full width)
Some additions can be made to make the resolution more accurate, but in the limit that there are many samples along the X axis and the data is not too noisy, this works great.
Even when the data are not Gaussian and a little noisy, it worked for me (I just take the first and last time half max crosses the data).