
% Script that measures the areas of two simulated overlapping 
% Gaussian peaks by several methods described in 
% https://terpconnect.umd.edu/~toh/spectrum/Integration.html
format compact
% User-adjustable parameters
increment=.01;
a1=1; % +
% Height of peak 1
a2=.5; % Height of peak 2
m1=40; % Position of peak 1
m2=57; % Position of peak 2
s1=6; % Width (sigma) of peak 1
s2=6; % Width (sigma) of peak 2

% Calculate simulated overlapping Gaussian peak pair
x=1:increment:100;
e=2.718;
g1=a1.*e.^(-(x - m1).^2/(2.*s1.^2))/(sqrt(2.*pi).*s1);
g1=g1./max(g1);
g2=a2.*e.^(-(x - m2).^2/(2.*s2.^2))/(sqrt(2.*pi).*s2);
g2=g2./max(g2);
g2=a2.*g2;
gs=g1+g2;
Area1=trapz(x,g1); % True area of peak 1
Area2=trapz(x,g2);% True area of peak 2
minindex=val2ind(x,10); % index value of start of first peak
maxindex=val2ind(x,90); % index value of end of second peak
w1=halfwidth(x,g1); % FWHM of peak 1
w2=halfwidth(x,g2); % FWHM of peak 2
Resolution=1.18*((m2-m1)./(w1+w2)); % Resolution definned by https://www.shimadzu.com/an/hplc/support/lib/lctalk/resol-1.html

% Perpendicular drop based on minimum valley between peaks
index1=val2ind(x,m1);
index2=val2ind(x,m2);
valleyy=min(gs(index1:index2));
valleyindex=val2ind(gs,valleyy);
valleyx=x(valleyindex);
PDMeasArea1=trapz(x(minindex:valleyindex),gs(minindex:valleyindex));
PDMeasArea2=trapz(x(valleyindex:maxindex),gs(valleyindex:maxindex));
PercentValley=100*valleyy/max(gs);
PDPercentError1=100.*(Area1-PDMeasArea1)/Area1;
PDPercentError2=100.*(Area2-PDMeasArea2)/Area2;

% Perpendicular drop based on half-way point between peaks
midindex=val2ind(x,(m1+m2)/2); % index value of halfway between peaks
HalfwayArea1=trapz(x(minindex:midindex),gs(minindex:midindex));
HalfwayArea2=trapz(x(midindex:maxindex),gs(midindex:maxindex));
HalfwayPercentError1=100.*(Area1-HalfwayArea1)/Area1;
HalfwayPercentError2=100.*(Area2-HalfwayArea2)/Area2;

% Normalization method
% Two-point straight line between the two peaks
x1=m1;
x2=m2;
y1=a1;
y1=gs(index1);
y2=a2;
y2=gs(index2);
cline=(x-x1).*(y2-y1)./(x2-x1)+y1;
for z=1:length(x)
    if cline(z)>y1,cline(z)=y1;end
    if cline(z)<y2,cline(z)=y2;end        
end
mg1=g1./cline;
mg2=g2./cline;
mgs=mg1+mg2;
minindex=val2ind(x,10); % index value of start of first peak
maxindex=val2ind(x,90); % index value of end of second peak
midindex=val2ind(x,(m1+m2)/2); % index value of halfway between peaks
index1=val2ind(x,m1);
index2=val2ind(x,m2);
valleymy=min(mgs(index1:index2));
valleymindex=val2ind(mgs,valleymy);
valleymx=x(valleymindex);
figure(1)
plot(x,gs,x,cline,':',x,gs./cline,'g--',[valleymx valleymx],[0 valleymy],'k-.','LineWidth',2)
title('Overlapping peak area measurement by geometrical height normalization')
legend('Observed signal','cline','(g1+g2)/cline','valley point')
axis([min(x) max(x) 0 1.1*a1])
grid
MeasArea1=trapz(x(1:valleymindex),mgs(1:valleymindex).*cline(1:valleymindex));
MeasArea2=trapz(x(valleymindex:index2),mgs(valleymindex:index2).*cline(valleymindex:index2))+trapz(x(index2:length(x)),mgs(index2:length(x)).*cline(index2:length(x)));
PercentMeasError1=100.*(MeasArea1-Area1)/Area1;
PercentMeasError2=100.*(MeasArea2-Area2)/Area2;

% Subtraction method
SignalMax=max(gs);
Maxindex=val2ind(gs,SignalMax);
Maxx=x(Maxindex);
% Estimate isolated peak 1 based on assumption of symmetry
rg1=[g1(1:Maxindex-1) fliplr(g1(1:Maxindex))];
rg1(2*Maxindex:length(x))=0;
figure(2)
clf
xx=x(1:length(rg1));
% Estimate isolated peak 2 by subtracting estimated peak 1 from total
% signal
rg2=gs(1:length(x))-rg1(1:length(x));
plot(x(1:length(x)),gs(1:length(x)),'b',xx,rg1,'r-.',xx,rg2,'g-.','LineWidth',2)
legend('Observed signal','Estimated peak 1','Calculated peak 2')
title('Overlapping peak area measurement by subtracting estimated first peak')

% Print table comparing the accuracy of peak area measurement.
disp(' ')
disp('Percent peak area errors by method')
disp(['Resolution: ' num2str(Resolution)  '    Heights:  ' num2str(a1) '/' num2str(a2)  '    Widths:  ' num2str(s1) '/' num2str(s2) ]);
disp(' ')
disp('Perpendicular drop, valley point')
disp([ 'Peak 1: ' num2str(PDPercentError1) '%     Peak 2: ' num2str(PDPercentError2) '% '])
disp(' ')
disp('Perpendicular drop, half-way point')
disp([ 'Peak 1: ' num2str(HalfwayPercentError1) '%     Peak 2: ' num2str(HalfwayPercentError2) '% '])
disp(' ')
disp('Geometrical peak height normalization')
disp([ 'Peak 1: ' num2str(PercentMeasError1) '%     Peak 2: ' num2str(PercentMeasError2) '% '])
disp(' ')
disp('subtracting estimated first peak')
disp([ 'Peak 1: ' num2str(PercentMeasError1) '%     Peak 2: ' num2str(PercentMeasError2) '% '])