% Comparison of four methods of overlapping peak pair area measurement for
% ecponentially broadened asymmetrical peaks:
% 1. Perpendicular drop (based on valley between peaks).
% 2. Perpendicular drop (based on and half-way points between peaks)
% 3. Peak height equalization and correction
% 4. Subtraction of estimated larger peak
% Requires val2ind.m, halfwidth.m, fastsmooth.m, ExpBroaden.m, and plotit.m functions
%
format compact
clear
increment=.001;
x=1:increment:100;
TimeConstant=.3; % Exponential broadening time constant, x units. Try .2 - 3
noise=0.001; % Try values from 0.1 or 0.001 (10% to 0.1%)
Peak1Variation=0; % fractional random variation of peak 1 amplitude
HeightsVector=[.99 .95 .9 .85 .8 .75 .7 .65 .6 .55 .5 .45 .4 .35 .3 .25 .2];
disp('Percent errors of peak 2 area measurement as a function of relative height')
disp('    Height2    Valley PD    Halfway PD    Equalization    Subtraction')
% Main loop
for trial=1:length(HeightsVector)
a1=1+Peak1Variation.*randn();
a2=HeightsVector(trial);
e=2.718;
m1=30; % Position of peak 1
m2=45; % Position of peak 2
s1=5; % width (sigma) of peak 1
s2=5;% width (sigma) of peak 2
separation=m2-m1; % Estimate of peak separation
g1=a1.*e.^(-(x - m1).^2/(2.*s1.^2))/(sqrt(2.*pi).*s1);
g1=ExpBroaden(g1',-TimeConstant./increment)';
g1=g1./max(g1);
w1=halfwidth(x,g1); % FWHM of peak 1
g2=a2.*e.^(-(x - m2).^2/(2.*s2.^2))/(sqrt(2.*pi).*s2);
g2=ExpBroaden(g2',-TimeConstant./increment)';
g2=g2./max(g2);
w2=halfwidth(x,g2); % FWHM of peak 2
g2=a2.*g2;
Resolution=1.18*((m2-m1)./(w1+w2)); % Resolution definned by https://www.shimadzu.com/an/hplc/support/lib/lctalk/resol-1.html
gs=g1+g2;
SmoothWidth=.1*s1./increment;
% gs=ExpBroaden(gs',-TimeConstant./increment)'; % Add exponential broadening
gs=gs+noise.*randn(size(x));
gs=fastsmooth(gs,SmoothWidth,3);

Area1=trapz(x,g1);
Area2(trial)=trapz(x,g2);

% Determine shifted peak positions after broadening, added noise, and smoothing
minindex=val2ind(x,10); % estimated index value of start of first peak
betweenindex=val2ind(x,40); % Somewhere between the peaks
maxindex=val2ind(x,90); % estimated index value of end of second peak
m1=max(gs(minindex:betweenindex));
m1index=val2ind(gs,m1);
m1=x(m1index);
m2=m1+separation; % Compensates for shift in m2 due to exponential broadening

% 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(trial)=trapz(x(valleyindex:maxindex),gs(valleyindex:maxindex));
PercentValley=100*valleyy/max(gs);
PDPercentError1=100.*(Area1-PDMeasArea1)/Area1;
PDPercentError2(trial)=100.*(Area2(trial)-PDMeasArea2(trial))/Area2(trial);

% 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(trial)=trapz(x(midindex:maxindex),gs(midindex:maxindex));
HalfwayPercentError1=100.*(Area1-HalfwayArea1)/Area1;
HalfwayPercentError2(trial)=100.*(Area2(trial)-HalfwayArea2(trial))/Area2(trial);

% Two-point straight line between the two peaks
x1=m1;
x2=m2;
y1=gs(index1);
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

mgs=gs./cline;
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,[valleymx(1) valleymx(1)],[0 valleymy],'k--','LineWidth',2)
title('Overlapping peak area measurement by geometrical height equalization')
legend('Observed signal','cline','Equalization','valley point')
xlabel( [' Resolution = ' num2str(Resolution) '   TimeConstant = ' num2str(TimeConstant) '     Noise = ' num2str(100.*noise/a1) '%' ])
axis([min(x) max(x) 0 1.1*a1])
grid

EMeasArea1=trapz(x(1:valleymindex),mgs(1:valleymindex).*cline(1:valleymindex));
EMeasArea2(trial)=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.*(EMeasArea1-Area1)/Area1;
PercentMeasError2(trial)=100.*(EMeasArea2(trial)-Area2(trial))/Area2(trial);

% Measure areas by reconstruction peak 1 by reflection and subtracting
% it from total signal to reveal peak 2.
SignalMax=max(gs(minindex:betweenindex));
Maxindex=val2ind(gs,SignalMax);
Maxx=x(Maxindex);
% Estimate isolated peak 1 based on assumption of symmetry
rg1=[gs(1:Maxindex-1) fliplr(gs(1:Maxindex))];
rg1(2*Maxindex:length(x))=0;

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));
figure(2)
clf
plot(x,gs,xx,rg1,'-.',xx,rg2,'-.','LineWidth',2)
legend('Observed signal','Estimated peak 1','Calculated peak 2')
title('Overlapping peak area measurement by subtracting estimated first peak')
xlabel( [' Resolution = ' num2str(Resolution) '   TimeConstant = ' num2str(TimeConstant) '     Noise = ' num2str(100.*noise/a1) '%' ])
axis([min(x) max(x) 0 1.1*a1])
grid

SignalMax2=max(rg2);
Maxindex2=val2ind(rg2,SignalMax2);
Maxx2=x(Maxindex2);

SMeasArea1=trapz(xx,rg1);
SMeasArea2(trial)=trapz(xx,rg2);
PercentAreaError1=100.*(Area1-SMeasArea1)./Area1;
PercentAreaError2(trial)=100.*(Area2(trial)-SMeasArea2(trial))./Area2(trial);

Heights(trial)=a2;
PD1(trial)=abs(PDPercentError2(trial));
PD2(trial)=abs(HalfwayPercentError2(trial));
GPHE(trial)=abs(PercentMeasError2(trial));
SEFP(trial)=abs(PercentAreaError2(trial));

disp([ '    ' num2str(a2) '       ' num2str(PDPercentError2(trial)) '%      '  num2str(HalfwayPercentError2(trial)) '%      '  num2str(PercentMeasError2(trial)) '%      ' num2str(PercentAreaError2(trial)) '% '])
end

figure(3)
clf
plot(HeightsVector,PD1,HeightsVector,PD2,HeightsVector,GPHE,HeightsVector,SEFP,'LineWidth',2)
grid
axis([min(HeightsVector) max(HeightsVector) 0 10])
xlabel( [' Resolution = ' num2str(Resolution) '   TimeConstant = ' num2str(TimeConstant) '     Noise = ' num2str(100.*noise/a1) '%' ])
ylabel('Absolute percent area measurement error')
title('Effect of the relative height of the smaller peak')
legend('Perp. drop valley','Perp drop halfway','Height equalization','Larger peak subtraction')

figure(4)
[coef1,RSquared1]=plotit(HeightsVector,PDMeasArea2,1,'o');
title('True vs measured peak 2 area by valley perpenducular drop method. ')
ylabel('Measured area of peak 2')
xlabel( [' Resolution = ' num2str(Resolution) '   TimeConstant = ' num2str(TimeConstant) '     Noise = ' num2str(100.*noise/a1) '%' ])

figure(5)
[coef2,RSquared2]=plotit(HeightsVector,HalfwayArea2,1,'o');
title('True vs measured peak 2 area by halfway perpenducular drop method. ')
ylabel('Measured area of peak 2')
xlabel( [' Resolution = ' num2str(Resolution) '   TimeConstant = ' num2str(TimeConstant) '     Noise = ' num2str(100.*noise/a1) '%' ])

figure(6)
[coef3,RSquared3]=plotit(HeightsVector,EMeasArea2,1,'o');
title('True vs measured peak 2 area by equalization method. ')
ylabel('Measured area of peak 2')
xlabel( [' Resolution = ' num2str(Resolution) '   TimeConstant = ' num2str(TimeConstant) '     Noise = ' num2str(100.*noise/a1) '%' ])

figure(7)
[coef4,RSquared4]=plotit(HeightsVector,SMeasArea2,1,'o');
title('True vs measured peak 2 area by subtraction method. ')
ylabel('Measured area of peak 2')
xlabel( [' Resolution = ' num2str(Resolution) '   TimeConstant = ' num2str(TimeConstant) '     Noise = ' num2str(100.*noise/a1) '%' ])

