% Comparison of iterative fitting of noisy overlapping peaks
% Requires peakfit.m function in the path, download from
% https://terpconnect.umd.edu/~toh/spectrum/peakfit.m
warning off MATLAB:polyfit:RepeatedPointsOrRescale
format compact

figure(1);clf

% Simulate data set
increment=1;
x=[1:increment:400];
% For each simulated peak, enter the amplitude, position, and width below
amp=[3 4 5];  % Amplitudes of the peaks
pos=[100 200 300];   % Positions of the peaks
wid=[20 40 70];   % Widths of the peaks
Noise=.2;
% A = matrix containing one of the unit-amplitude peak in each of its rows
A = zeros(length(pos),length(x));
for k=1:length(pos)
    if amp(k)>0, A(k,:)=gaussian(x,pos(k),wid(k)); end
end
z=amp*A;  % Multiplies each row by the corresponding amplitude and adds them up
y=z+Noise.*randn(size(z));

% Fit using hand-coded fitshape2 function
start=[pos(1) wid(1) pos(2) wid(2) pos(3) wid(3)];
disp('Hand-coded function:')
tic
[Positions,Heights,Widths]=fitshape2(x,y,start);
toc

% Assemble best-fit model from the peak parameters
for k=1:length(pos)
    if amp(k)>0, A(k,:)=gaussian(x,Positions(k),Widths(k)); end
end
z=Heights*A;  % Multiplies each row by the corresponding amplitude and adds them up
% Plot original data and best fit in lower panel
subplot(2,1,1)
plot(x,y,'.',x,z)
title('Hand-coded fitshape2.m function, with default input arguments. (No first guesses supplied).')

% Fit using ChatGPT-created iterativefitGaussians function
n=length(amp);
initPositions=pos;
initWidths=wid;
initHeights=[1 1 1];
disp('ChatGPT-coded function:')
tic
[positions,widths,heights]=iterativefitGaussians(x,y,n,initPositions,initWidths);
toc
% Assemble best-fit model from the peak parameters
for k=1:length(pos)
    if amp(k)>0, A(k,:)=gaussian(x,positions(k),(widths(k))./0.6); end
end
z=heights*A;  % Multiplies each row by the corresponding amplitude and adds them up
% Plot original data and best fit in upper panel
subplot(2,1,2)
plot(x,y,'.',x,z)
title('Coded by ChatGPT, with specified first guesses for peak positions and widths.')

%----------------------------- local functions -------------------

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.60056120439323.*wid)) .^2);
end

function [Positions,Heights,Widths]=fitshape2(x,y,start)
global PEAKHEIGHTS
% Call the fminsearch function with specified fitting function
FitResults=fminsearch(@(lambda)(fitfunction(lambda,x,y)),start);
NumPeaks=round(length(start)./2);

title('fitshape2')
for m=1:NumPeaks
    Positions(m)=FitResults(2*m-1);
    Widths(m)=FitResults(2*m);
end
Heights=PEAKHEIGHTS';
end

function err = fitfunction(lambda,t,y)
%  Fitting function for a signal consisting of overlapping peaks.
A = zeros(length(t),round(length(lambda)/2));
for j = 1:length(lambda)/2
    A(:,j) = shapefunction(t,lambda(2*j-1),lambda(2*j))';
end
PEAKHEIGHTS = A\y';
z = A*PEAKHEIGHTS;
err = norm(z-y');
end

function g = shapefunction(x,a,b)
%  shapefunction(x,pos,wid) = peak peak centered on pos, half-width=wid
%  x may be scalar, vector, or matrix, pos and wid both scalar
%  T. C. O'Haver, 1988
% Replace the last line with the expression for your shape,
% similar to the gaussian definition originally there.
g = exp(-((x-a)./(0.6005612.*b)) .^2); % Change this to change shape
end

function [positions, widths, heights] = iterativefitGaussians(x, y, n, initPositions, initWidths)
% ChatGPT: Fit n Gaussian functions to the data given in x and y
%   [positions,widths,heights]=fitGaussians(x,y,n,initPositions,initWidths)
%   fits the sum of n Gaussian functions to the data given in x and y. The
%   initial estimates of the positions and widths of the Gaussians are given
%   in the vectors initPositions and initWidths. The function returns the
%   position, width, and height of all the Gaussians.

% Initialize variables
positions = initPositions;
widths = initWidths;
heights = zeros(1, n);

% Set optimization options
options = optimoptions('lsqcurvefit', 'Display', 'off');

% Fit the Gaussians
for i = 1:n
    % Define the function to fit
    fitFcn = @(p, x) p(1) .* exp(-((x - p(2))./p(3)).^2);
    
    % Define the initial values of the parameters
    initParams = [y(i), positions(i), widths(i)];
    
    % Define the parameter bounds
    lb = [0, min(x), 0];
    ub = [inf, max(x), inf];
    
    % Fit the Gaussian
    params = lsqcurvefit(fitFcn, initParams, x, y, lb, ub, options);
    
    % Extract the position, width, and height of the Gaussian
    positions(i) = params(2);
    widths(i) = params(3);
    heights(i) = params(1);
end
end



