function [MeasuredSTD ResultsMatrix]=GaussFitMC2(TrueHeight,TruePosition,TrueWidth,NumPoints,Noise,NumTrials)
% Monte Carlo simulation of the measurement of the peak height, position,
% and width of a noisy x,y Gaussian peak, comparing the gaussfit parabolic 
% fit to the fitgaussian iterative fit (included as subfunctions). 
% Plots the peak height, position, and width of the gaussfit results 
% vs the height, position, and width of the fitgaussian results, 
% plus a randomly selected single data set and Gaussian fit.
% The input arguments are:
% TrueHeight = True peak height of the Gaussian peak.
% TruePosition = True x-axis value at the peak maximum.
% TrueWidth =  True half-width (FWHM) of the Gaussian peak.
% NumPoints = Number of points taken for the least-squares fit
% Noise = standard deviaiton of (normally-distributed) random noise
% NumTrials = Number of trials used to compute standard deviations
% Optional Outputs:
% MeasuredSTD = computed standard deviation of Measured [gaussfitHeight
% gaussfitPosition gaussfitWidth fitgaussianHeight fitgaussianPosition fitgaussianWidth]
% ResultsMatrix = Measured [gaussfitHeight gaussfitPosition gaussfitWidth fitgaussianHeight fitgaussianPosition fitgaussianWidth] of each trial
% Example: GaussFitMC2(100,100,100,50,1,1000)
% ans =
%     0.2135    0.1403    0.4756    0.2103    0.1299    0.4590
global PEAKHEIGHTS
startx=TruePosition-TrueWidth/2;
endx=TruePosition+TrueWidth/2;
deltax=(endx-startx)/(NumPoints-1);
x=startx:deltax:endx;
ResultsMatrix=zeros(6,NumTrials);
for trial=1:NumTrials,
    NoiseArray=Noise*randn(size(x));
    y=TrueHeight*gaussian([startx:deltax:endx],TruePosition,TrueWidth)+NoiseArray;
    [Height, Position, Width]=gaussfit(x,y);
    start=[Position Width];    % generic initial guesses for peak position and width
    options = optimset('TolX',0.1);  % Determines how close the model must fit the data
    parameter=fminsearch(@(lambda)(fitgaussian(lambda,x,y)),start);
    ResultsMatrix(:,trial)=[Height, Position, Width, PEAKHEIGHTS, parameter(1), parameter(2)];
end
FullDataSetSTD=std(ResultsMatrix');
% Display results
MeasuredSTD=FullDataSetSTD;
subplot(2,2,1)
plotit(ResultsMatrix(1,:),ResultsMatrix(4,:),1);
title('Peak Heights');
xlabel('Parabolic fit')
ylabel('Interative fit')
axis([98.5 101.5 98.5 101.5])
subplot(2,2,2)
plotit(ResultsMatrix(2,:),ResultsMatrix(5,:),1);
title('Peak Positions');
xlabel('Parabolic fit')
ylabel('Interative fit')
axis([98.5 101.5 98.5 101.5])
subplot(2,2,3)
plotit(ResultsMatrix(3,:),ResultsMatrix(6,:),1);
title('Peak Widths');
xlabel('Parabolic fit')
ylabel('Interative fit')
axis([98.5 101.5 98.5 101.5])
subplot(2,2,4)
plot(x,y,'o',linspace(min(x),max(x)),Height*gaussian(linspace(min(x),max(x)),Position,Width))
title('Single data set');
xlabel('x')
ylabel('y')

function [Height, Position, Width]=gaussfit(x,y)
% Converts y-axis to a log scale, fits a parabola
% (quadratic) to the (x,ln(y)) data, then calculates
% the position, width, and height of the
% Gaussian from the three coefficients of the
% quadratic fit.  This works only if the data have
% no baseline offset (that is, trends to zero far off the
% peak) and if there are no zeros or negative values in y.
% Example: [Height, Position, Width]=gaussfit([1 2 3],[1 2 1]) 
% returns Height = 2, Position = 2, Width = 2
z=log(y);
coef=polyfit(x,z,2);
a=coef(3);
b=coef(2);
c=coef(1);
Height=exp(a-c*(b/(2*c))^2);
Position=-b/(2*c);
Width=2.35482/(sqrt(2)*sqrt(-c));

function err = fitgaussian(lambda,t,y)
% Fitting function for a Gaussian band signal.
global PEAKHEIGHTS
numpeaks=round(length(lambda)/2);
A = zeros(length(t),numpeaks);
for j = 1:numpeaks,
    A(:,j) = gaussian(t,lambda(2*j-1),lambda(2*j))';
end 
PEAKHEIGHTS = abs(A\y');
z = A*PEAKHEIGHTS;
err = norm(z-y');

function g = gaussian(x,pos,wid)
%  gaussian(x,pos,wid) = gaussian peak centered on pos, half-width=wid
%  x may be scalar, vector, or matrix, pos and wid both scalar
%  T. C. O'Haver, 1988
% Examples: gaussian([0 1 2],1,2) gives result [0.5000    1.0000    0.5000]
% plot(gaussian([1:100],50,20)) displays gaussian band centered at 50 with width 20.
g = exp(-((x-pos)./(0.6005615.*wid)) .^2);