# Is there a stable numerical algorithm for FWHM that isn't 2.35*sigma?

This is a question that should have a simple answer, but which I can find no proper discussion of in the literature or on the internet.

I start from the assumption that I have a noisy numerical signal with a peak in it - for example amplitude as a function of $x$, i.e. $A(x)$ - and I wish to determine its full-width at half maximum (FWHM). I am aware of two basic algorithms that may be used on numerical data:

1. Assume it's Gaussian, determine the r.m.s. value, and multiply by 2.35 ( $2\sqrt{2 \log{2}}$).

2. Determine the FWHM using a peak-finding algorithm that locates the peak $A_{max} = A(x_{peak}$), then locates the first positions either side where $A(x)$ falls to $1/2 A_{max}$.

Method 1 is numerically robust as you don't have to bin (i.e. smooth) the data that determines $A(x)$. However, if the peak is significantly non-Gaussian in shape you get a systematically wrong answer.

Method 2 is direct, but there are plenty of numerical cases where there are multiple points in $x$ either side of the peak where $A(x)$ crosses $1/2 A_{max}$. And of course the number of crossing points depends sensitively on whether I smooth A(x) or not.

My question is therefore whether there is a mathematically justifiable algorithm for determining FWHM numerically that doesn't assume the peak is Gaussian,

-
I don't quite understand you method 1. Do you mean, fit a gaussian peak to the data, determine its RMS value...? If you determine the RMS of the experimential data, you will get completely useless results if there's a lot of noise, because this will of cause also be taken into account. – leftaroundabout Aug 27 '11 at 19:54
Good question and answer, but not physics. – Colin K Aug 27 '11 at 21:22
I'm of a mind to migrate this to Statistics.SE. Any objections? – dmckee Aug 27 '11 at 22:41
I think it belongs in physics. This is the sort of thing (experimental) physicists do all day; this is where physicists would search for an answer. They speak a different language in CrossValidated. – nibot Aug 28 '11 at 19:46
I personally find this kind of question a refreshing change from the usual soft questions we get in physics.SE. – nibot Aug 28 '11 at 19:47

Typically one knows the functional form of something when they are interested in the FWHM. In this case, you can use least-squares fitting to fit the function to the data and extract the parameters. A very common functional form of a resonance is the Lorentzian:

$$f(A,\gamma,x_0;x) = A\frac{1}{1 + \gamma^2\left(x-x_0\right)^2}$$ where $1/\gamma$ gives the HWHM. Sometimes you will have to tailor this to your data a little bit; for instance, while this function goes to zero as $x\to\pm\infty$, your real data will probably have some artificial offset, which you would add a term to account for. Also, of course, you will want to either exclude or simultaneously fit any other nearby resonances.

Here is some code that demonstrates this in Matlab:

% define the functional form of the resonance:
lorentzian = @(a, x) a(3)./(1 + a(2)^2*(x - a(1)).^2);

% make a signal and noise:
x = linspace(0, 100, 101);
signal = lorentzian([70 0.1 1], x);
noise = 0.1*randn(size(signal));
y = signal + noise;

% make initial guess of the parameters:
[m, ii] = max(y);
a0 = [x(ii) m 1];

% do a least-squares fitting of the lorentzian to the measured data
a = lsqcurvefit(lorentzian, a0, x, y);

% plot the results
plot(x, y, 'o');
hold all
plot(x, signal, '--');
plot(x, lorentzian(a, x), '-');
hold off

legend('data', 'true signal',  'best fit to data', ...
'Location', 'Best');


Despite a fairly large amount of scatter in the data points, the least squares fit recovers the parameters of the curve quite well:

                          x0       gamma        A
true value                70         0.1        1
estimated value           70.33      0.0956     0.97

-
Great answer. In general, you certainly need to assume something about the functional form, or else there's no sensible way to answer the question. After all, without some assumption, it's logically possible that the true shape wiggles in all sorts of complicated ways. If you don't want to assume a fully-specified functional form, I suppose you could fit to some broad class of smooth functions. If you think that the true shape is (in some sense) close to Gaussian, for instance, you might want to look into the Edgeworth expansion. – Ted Bunn Aug 28 '11 at 17:02