How to fit a gaussian to data in matlab/octave?

Question

I have a set of frequency data with peaks to which I need to fit a Gaussian curve and then get the full width half maximum from. The FWHM part I can do, I already have a code for that but I'm having trouble writing code to fit the Gaussian.

Does anyone know of any functions that'll do this for me or would be able to point me in the right direction? (I can do least squares fitting for lines and polynomials but I can't get it to work for gaussians)

Also it would be helpful if it was compatible with both Octave and Matlab as I have Octave at the moment but don't get access to Matlab until next week.

Any help would be greatly appreciated!

Answer 1

Fitting a single 1D Gaussian directly is a non-linear fitting problem. You'll find ready-made implementations here, or here, or here for 2D, or here (if you have the statistics toolbox) (have you heard of Google? :)

Anyway, there might be a simpler solution. If you know for sure your data y will be well-described by a Gaussian, and is reasonably well-distributed over your entire x-range, you can linearize the problem (these are equations, not statements):

   y = 1/(σ·√(2π)) · exp( -½ ( (x-μ)/σ )² )
ln y = ln( 1/(σ·√(2π)) ) - ½ ( (x-μ)/σ )²
     = Px² + Qx + R         

where the substitutions

P = -1/(2σ²)
Q = +2μ/(2σ²)    
R = ln( 1/(σ·√(2π)) ) - ½(μ/σ)²

have been made. Now, solve for the linear system Ax=b with (these are Matlab statements):

% design matrix for least squares fit
xdata = xdata(:);
A = [xdata.^2,  xdata,  ones(size(xdata))]; 

% log of your data 
b = log(y(:));                  

% least-squares solution for x
x = A\b;                    

The vector x you found this way will equal

x == [P Q R]

which you then have to reverse-engineer to find the mean μ and the standard-deviation σ:

mu    = -x(2)/x(1)/2;
sigma = sqrt( -1/2/x(1) );

Which you can cross-check with x(3) == R (there should only be small differences).