% Measuring a weak peak next to a much larger peak of same basic shape.
% Demonstration of increasing measurement accuracy by position and width
% constraints. Hint: Spread out the figure windows so they can all be seen
% separately and don't overlap. You can change the constants in lines 11 -
% 20 or the peak shape in line 26. Try SmallHeight = .1 and .05; Noise = 0
% and .01; xshift = 0 and 1. Requires gaussian.m, peakfit.m and findpeaks.m
% in the path.
format compact
format short g
clear
increment=.05; % x-axis increment (smaller => more data points)
P1=4; % Peak position of first (larger) interfering peak
P2=5; % Peak position of second (smaller) measured peak
LargeHeight=1; % Height of peak 1 (the uninteresting or interfering peak)
SmallHeight=.05; % Height of peak 2 (the measured peak)
width=1.66; % Full width at half maximum of both peaks
Noise=.01; % Fractional random white noise
xshift=0; % Average unexpected shift in x-axis, causing both peaks to shift
NumTrials=5; % Number of trial fits per signal from different starting points
NumSignals=20; % Number of repeated measurements with independent noise
SignalToNoiseRatio=SmallHeight/Noise % SNR of smaller peak
x=0:increment:10;n=length(x);
for trial=1:NumSignals,
    % Compute signal vector y for this measurement (change shape in line 26)
    shift=xshift.*randn(); % Shift in position of peaks for this measurement
    y=SmallHeight.*gaussian(x,P2+shift,width)+LargeHeight.*gaussian(x,P1+shift,width)+Noise.*randn(size(x));
    
    figure(1)
    F0=peakfit([x;y],0,0,2,1,0,NumTrials,0,0);
    F0=sortrows(F0,3);
    subplot(2,1,1)
    title('Unconstrained iterative fit (variable positions and widths), using peakfit.m: ')
    
    figure(2)
    F1=peakfit([x;y],0,0,2,6,0,NumTrials);
    F1=sortrows(F1,3);
    subplot(2,1,1)
    title('Equal-width iterative fit (variable positions and equal widths), using peakfit.m')
    
    figure(3)
    F2=peakfit([x;y],0,0,2,16,0,NumTrials,0,0,[P1 P2]);
    F2=sortrows(F2,3);
    subplot(2,1,1)
    title('Fixed positions iterative fit (variable widths), using peakfit.m')
 
    figure(4)
    CLSheights=cls(x,y,2,1,[4 5],[width width],0);
    clsp1=LargeHeight.*exp(-(x-P1).^2);
    clsp2=SmallHeight.*exp(-(x-P2).^2);
    clsM=clsp1+clsp2;
    subplot(2,1,1)
    plot(x,y,'b.',x,clsp1,'g',x,clsp2,'g',x,clsM,'r')
    axis([0 10 0 1.2])
    FittingError=100*norm(y-clsM)./(sqrt(n)*max(y));
    xlabel(['Peakheights = '  num2str(CLSheights) '        Fitting error = ' num2str(FittingError) '%' ])
    title('Fixed positions and widths, using CLS.m')
    subplot(2,1,2)
    plot(x,y-clsM,'m.')
    xlabel('Residual plot')
    
    figure(5)
    % Find peak position of this measurement (P(2))
    SlopeThreshold=.003;
    AmpThreshold=.1;
    SmoothWidth=3;
    FitWidth=11;
    Peaktable=findpeaksG(x,y,SlopeThreshold,AmpThreshold,SmoothWidth,FitWidth,3); 
    CLSheights2=cls(x,y,2,1,[Peaktable(2) Peaktable(2)+1],[width width],0);
    % Assume peak 1 has this position
    clsp1=LargeHeight.*exp(-(x-Peaktable(2)).^2);
    % Assume peak 2 is also shifted the same amount in this measurement
    clsp2=CLSheights2(2).*exp(-(x-Peaktable(2)-(P2-P1)).^2); 
    clsM=clsp1+clsp2;
    subplot(2,1,1)
    plot(x,y,'b.',x,clsp1,'g',x,clsp2,'g',x,clsM,'r')
    axis([0 10 0 1.2])
    FittingError=100*norm(y-clsM)./(sqrt(n)*max(y));
    xlabel(['Peakheights = '  num2str(CLSheights2) '        Fitting error = ' num2str(FittingError) '%' ])
    title('Position 1 calculated from findpeaks, fixed widths and distance between peaks, using CLS.m')
    subplot(2,1,2)
    plot(x,y-clsM,'m.')
    xlabel('Residual plot')
    
     figure(6)
     start=[Peaktable(2) width Peaktable(2)+(P2-P1) width];
     F5=peakfit([x;y],0,0,2,6,0,NumTrials,start);
     F5=sortrows(F5,3);
     subplot(2,1,1)
     title('Equal-width iterative fit (first guess derived from findpeaks), using peakfit.m')
    
    e0(trial)=100.*(SmallHeight-F0(1,3))./SmallHeight;
    e1(trial)=100.*(SmallHeight-F1(1,3))./SmallHeight;
    e2(trial)=100.*(SmallHeight-F2(1,3))./SmallHeight;
    e3(trial)=100.*(SmallHeight-CLSheights(2))./SmallHeight;
    e4(trial)=100.*(SmallHeight-CLSheights2(2))./SmallHeight;
    e5(trial)=100.*(SmallHeight-F5(1,3))./SmallHeight;
end
Unconstrained=sqrt(norm(e0));
EqualWidths=sqrt(norm(e1));
FixedPositions=sqrt(norm(e2));
FixedPositionsAndWidths=sqrt(norm(e3));
findpeaksP=sqrt(norm(e4));
findpeaksP2=sqrt(norm(e5));
disp('      Unconstr.     EqualW       FixedP      FixedP&W    findpeaksP   findpeaksP2')
disp([Unconstrained EqualWidths FixedPositions FixedPositionsAndWidths findpeaksP findpeaksP2])