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

Question

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,

Answer 1

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