% Method of one-step least-squares fit to a single gaussian.
% Generates a noisy Gaussian function, converts the
%  y-axis to a log scale, fits a parabola
% (quadratic) to the (x,log(y)) data, then calculates
% the position, width, and height of the original
% Gaussian from the three coefficients of the
% quadratic fit.  This works only if the Gaussian has
% no baseline offset (that is, goes to zero far off the
% peak) and if there are no negative values in y.

format compact
% Simulation section
% Generates a noisy Gaussian function
x=[25:31:200];
Height=100;
Width=100;
Position=100;
Noise=1;
y=Height.*gaussian(x,Position,Width)+Noise.*randn(size(x));

% Analysis section;
z=log(abs(y));
coef=polyfit(x,z,2);
x2=[-50:300];
fit=polyval(coef,x2);
a=coef(3);
b=coef(2);
c=coef(1);
Height
MeasuredHeight=exp(a-c*(b/(2*c))^2);
MeasuredPosition=-b/(2*c);
MeasuredWidth=2.35703/(sqrt(2)*sqrt(-c));
HeightResults=[Height MeasuredHeight]
PositionResults=[Position MeasuredPosition]
WidthResults=[Width MeasuredWidth]

% Plotting section
subplot(2,2,1) % Upper left segment
plot(x,y,'ro')
axis([-50,300,0,120]);
ylabel('y')
title('Original data')

subplot(2,2,2) % Upper right segment
plot(x,z,'ro')
axis([-50,300,2,5]);
ylabel('ln(y)')
title('Ln of original data')

subplot(2,2,3) % Lower left segment
plot(x,z,'ro',x2,fit)
axis([-50,300,2,5]);
ylabel('ln(y)')
title('Ln of original data with quadratic fit')
xlabel(['a = '  num2str(coef(3))  '    b = '   num2str(coef(2))  '    c = '  num2str(coef(3))] )

subplot(2,2,4) % Lower right segment
plot(x,y,'Ro',x2,Height.*gaussian(x2,MeasuredPosition,MeasuredWidth),'b')
axis([-50,300,0,120]);
xlabel(['Height = '  num2str(MeasuredHeight)  '    Position = '   num2str(MeasuredPosition)  '    Width = '  num2str(MeasuredWidth)] )
ylabel('y')
title('Original data with computed Gaussian fit')
'     Height   Position   Width'
PercentErrors=[100*(MeasuredHeight-Height)/Height 100*(MeasuredPosition-Position)/Position 100*(MeasuredWidth-Width)/Width]