% Demonstration of the linearization of the power transform method by
% taking the nth power of data and then the 1/n power of the measured
% areas. Simulates two overlapping peaks with specified amplitudes. Peak
% areas are measured by perpendicular drop, using the half-way point to
% mark the boundary between the peaks. You can change the parameters in
% lines 14-22. Requires gaussian.m, halfwidth.m, val2ind.m, and plotit.m
% downloaded from my web site: http://tinyurl.com/cey8rwh

%
format compact
format short g

% User changable constants
n=10; % <<<< CHANGE POWER HERE
k1=2; % <<<< "Unknown concentration" of component 1
k2=4; % <<<< "Unknown concentration" of component 2
width=3.5; % Initial peak width (FWHM).
sigma=width/(2*log(2)); % peak width (sigma)
position1=4; % x-axis position of peak 1
position2=8.5; % x-axis position of peak 2
Resolution=(position2-position1)./(4*sigma); % Original 4-sigma resolution
noise=.01; % Random noise standard deviation (optional)

% Main processing loop
figure(1)
clf
disp('Calibration results:')
disp('        Power  Resolution   Slope 1    Slope 2         R2')
x=-2:.01:15;
clf
for k=1:1:20 % k is trial counter for changing standard mixtures.
   % Choose between two ways of preparing and running the standards
%   c1(k)=1+rand(); % Randomized order of standard mixtures
%   c2(k)=1+rand(); % Randomized order of standard mixtures
    c1(k)=k/20; % Fixed ratio of standard mixtures
    c2(k)=k/20; % Fixed ratio of standard mixtures
    % Component 1 gives 2x signal strength for the same concentration
    y1=c1(k).*2.*gaussian(x,position1,width); % Peak 1 (variable height)
    y2=c2(k).*gaussian(x,position2,width);% Peak 2 (variable height)
    y=y1+y2+noise*randn(size(x)); % Overlapping peak 1 and 2 signal
    yp=y.^n; % Raise to power n
    % Compute resolution after power raising
    Resolutionp=(position2-position1)./(2.8853*halfwidth(x,y1.^n,position1));
    hold on
    Peak1MaxIndex=val2ind(x,position1); % Index of maximum peak 1
    Peak2MaxIndex=val2ind(x,position2); % Index of maximum peak 2
    valleyy=min(yp(Peak1MaxIndex:Peak2MaxIndex)); % Y value at valley
    valleyx=x(val2ind(yp,valleyy)); % X value at valley
    
    % Two alternatives ways to compute boundary between peaks
    % 1. Look for actual valley between peaks
    % ValleyIndex=val2ind(x,valleyx); % Index of valley point
    % 2. Compute half-way point between peaks, if peak positions are known
    ValleyIndex=val2ind(x,position1+(position2-position1)/2);
    
    figure(1)
    subplot(2,1,1)
    hold on
    plot(x,y);
    xlabel('x (time)')
    ylabel('y')
    hold off
    title(['Signals of the standard solutions.  Power = ' num2str(1) '    Resolution = ' num2str(Resolution) ])
    subplot(2,1,2)
    hold on
    plot(x,yp);
    xlabel('x (time)')
    ylabel('y^p')
    hold off
    title(['Power-raised signals of the standard solutions.  Power = ' num2str(n) '    Resolution = ' num2str(Resolutionp) ])
    
   ap1(k)=sum(yp(1:ValleyIndex))*(x(2)-x(1)); %  perpendicular drop areas of peak 1
    ap2(k)=sum(yp(ValleyIndex:length(x)))*(x(2)-x(1)); %  perpendicular drop areas of peak 2
    apr1(k)=ap1(k).^(1/n); % nth root of perpendicular drop areas of peak 1
    apr2(k)=ap2(k).^(1/n); % nth root of perpendicular drop areas of peak 2
end
xlabel('x (time)')
ylabel('y^p')
%     hold off

% Plot calibration curve for component 1 and fit to straight line
figure(2)
clf
[coef1,RSquared1]=plotit(c1,apr1,1,'or','b-');
xlabel('k, concentration of standard solution with fixed ratio of components')
ylabel('1/n power of area of signal raised to n power')
title(['Calibration curve 1 of the Power Method. Power = ' num2str(n) '   Resolution = ' num2str(Resolution) ])
% disp(['Peak 1   ' num2str([n Resolutionp coef(1) RSquared]) ])

% Plot calibration curve for component 2 and fit to straight line
figure(3)
clf
[coef2,RSquared2]=plotit(c2,apr2,1,'or','b-');
xlabel('k, concentration of standard solution with fixed ratio of components')
ylabel('1/n power of area of signal raised to n power')
title(['Calibration curve 2 of the Power Method. Power = ' num2str(n) '   Resolution = ' num2str(Resolution) ])
disp(['         ' num2str([n Resolutionp coef1(1) coef2(1) mean([RSquared1(1) RSquared2(1)]) ]) ])

% Measurement of unknown sample with specified concentrations k1 and k2

% Generate signal with both components, plus noise, and raise to power n.
s1=k1*2.*gaussian(x,position1,width); % Peak 1 (variable height)
s2=k2.*gaussian(x,position2,width);% Peak 2 (variable height)
s=s1+s2+noise*randn(size(x)); % Overlapping peak 1 and 2 signal
sp=s.^n; % Raise to power n

% Plot signal of sample mixture.
figure(4)
plot(x,sp)
xlabel('x (time)')
ylabel('y^p')
title(['Signal with component concentrations k1 = ' num2str(k1) ' and k2 = ' num2str(k2) ])

% Measure the areas of the two peaks and take nth root.
as1=sum(sp(1:ValleyIndex))*(x(2)-x(1)); %  perpendicular drop areas of peak 1
as2=sum(sp(ValleyIndex:length(x)))*(x(2)-x(1)); %  perpendicular drop areas of peak 2
asr1=as1.^(1/n); % nth root of perpendicular drop areas of peak 1
asr2=as2.^(1/n); % nth root of perpendicular drop areas of peak 2

% Solution of calibration equation for concentration of unknown: Cx = (Sx - intercept)/slope
slope1=coef1(1);
slope2=coef2(1);
intercept1=coef1(2);
intercept2=coef2(2);
disp('Measured values and errors for components 1 and 2 in the unknown sample')
Measured1=(asr1-intercept1)/slope1;
PercentError1=100*(Measured1-k1)/k1;
Measured2=(asr2-intercept2)/slope2;
PercentError2=100*(Measured2-k2)/k2;
disp('          Concentration  % error')
disp(['Component 1: ' num2str(Measured1)  '   ' num2str(PercentError1) '%' ])
disp(['Component 2: ' num2str(Measured2)  '   ' num2str(PercentError2) '%' ])
disp(['Calibration curve slope ratio: ' num2str(slope1/slope2) ])
disp(' ')