% Demonstration script (for peakfit.m, version 1.4). Generates an
% overlapping peak signal, adds noise, fits it with peakfit.m , then
% computes the accuracy and precision of peak paramater measurements. 
% You can changeany of the initial values in lines 13-30. T. C. O'Haver 
% (toh@umd.edu) Ver. 6: April, 2011, for use with peakfit ver. 1.4 (2011)

clear
format long g
format compact
warning off all

% >>> Initial settings of user-adjustable parameters <<<
n=300;  % maximum x-value
increment=1; % Difference between adjacent x values
x=[1:increment:n];  % Creates x-vector
amp=[1 2 3];  % Amplitude  s of the peaks (one entry for each peak)
pos=[110 140 180];   % Positions of the peaks (in increasing order)
wid=[25 25 25];   % Widths of the peaks (full width at half-maximum)
noise=.05; % Amplitude of random noise added to signal
NumRepeats=2; % Number of repeat samples with independent noise
NumTrials=10; % Number of trail fits per sample
PeakShape=1; % Gaussian=1, Lorentzian=2, Logistic=3, Pearson=4
% ("extra" determines the shape of the Pearson function: as extra
% goes from 1 -> large, shape goes from Lorentzian -> Gaussian).  
% Exponentially-broadened Gaussian = 5; Equal-widths Gaussians = 6;
% Equal-widths Lorentzians = 7. 
extra=2; 
center=150; % x-value in center of fitted range
window=200; % range of x-values in fitted range
startvector=[100 30 150 30 200 30]; % First guess vector [pos1 wid1 pos2 wid2 ...] Must be two values for each peak

% Generate simulated signal peaks
% A = matrix containing one of the unit-amplidude functions in each of its rows
clear A
A=zeros(length(pos),length(x));
for k=1:length(pos),
  switch PeakShape
  case 1
    A(k,:)=exp(-((x-pos(k))./(0.6006.*wid(k))) .^2);  % Gaussian function
  case 2
    A(k,:)=ones(size(x))./(1+((x-pos(k))./(0.5.*wid(k))).^2); % Lorentzian function
  case 3
    L = exp(-((x-pos(k))/(.477.*wid(k))) .^2);
    A(k,:)= (2.*L)./(1+L);
  case 4
      % pearson function (shape determined by value of extra in line 21
      A(k,:)=ones(size(x))./(1+((x-pos(k))./((0.5.^(2/extra)).*4.62.*wid(k))).^2).^extra;
  case 6
    A(k,:)=exp(-((x-pos(k))./(0.6006.*wid(k))) .^2);  % Gaussian function
  case 7
    A(k,:)=ones(size(x))./(1+((x-pos(k))./(0.5.*wid(k))).^2);  % Laorentzian function
  end % switch
  area(k)=amp(k).*trapz(x,A(k,:)); % Computes peak areas via trapaziod method
end % for
TrueY=amp*A;  % Multiplies each row by the corresponding amplitude and adds them up
% To test the autozero baseline subtraction, un-comment next line.
% TrueY=TrueY+ones(size(x))./(1+((x+100)./(300)).^2); % Adds sloping curved baseline

% Arrange the peak parameters of the simulated signal into a matrix
ColumnLabels='           Peak#   Position     Height     Width       Area';
NumPeaks=length(pos);
for m=1:NumPeaks,
  if m==1,  
     ActualParameters=[m pos(m) amp(m) wid(m) area(m) ]; 
  else
     ActualParameters=[ActualParameters; [m pos(m) amp(m) wid(m) area(m) ]];
  end
end

% Generate and fit NumRepeats separate samples of the simulated signal with noise.
AllFitResults=zeros(NumPeaks,5,NumRepeats);
AllFittingErrors=zeros(1,NumRepeats);
AllPercentParameterErrors=zeros(NumPeaks,5,NumRepeats);
for repeat=1:NumRepeats,
   %  Add normally-distributed random noise to simulated signal
   y = TrueY + noise .* randn(size(x));
   % Execute the peakfit script
   [FitResults,LowestError]=peakfit([x' y'],center,window,NumPeaks,PeakShape,extra,NumTrials,startvector); 
   FitResults=sortrows(FitResults,2); % Sort results by increasing peak positions (column 2)
   AllFitResults(:,:,repeat)=FitResults;
   AllFittingErrors(repeat)=LowestError;
   AllPercentParameterErrors(:,:,repeat)=100.*(FitResults-ActualParameters)./ActualParameters;
end

% Calculate the mean and relative standard deviaiton of peak parameters
% measured by curve fitting. \
if NumRepeats>1, 
   for peak=1:NumPeaks,
     MeanParameters(peak,:)=mean(squeeze(AllFitResults(peak,:,:))');
     RSDparameters(peak,:)=100.*std(squeeze(AllFitResults(peak,:,:))')./mean(squeeze(AllFitResults(peak,:,:))');
   end
else
    MeanParameters=AllFitResults;
end
  disp('-------------------------------------------------------')
    switch PeakShape
    case 1
        ShapeString='1 Gaussian';
    case 2
        ShapeString='2 Lorentzian';    
    case 3
        ShapeString='3 Logistic';
    case 4
        ShapeString='4 Pearson';
    case 5
        ShapeString='5 ExpGaussian';
    case 6
        ShapeString='6 Equal width Gaussians';
    case 7
        ShapeString='7 Equal width Lorentzians';
      otherwise
  end % switch peakshape
  disp([ 'Peak shape: ' ShapeString ])
  disp(['Number of peaks: ' num2str(NumPeaks)])
  disp(['Noise amplitude: ' num2str(noise)])
  disp(['Number of data points in signal: ' num2str(length(x))])
  disp(['Number of repeat signals: ' num2str(NumRepeats)])
  disp(['Number of trial fits per sample signal: ' num2str(NumTrials)])
  
  disp('Actual parameters:');
  disp('Position      Height       Width       Area');
  disp(num2str(ActualParameters(:,2:5)))
  disp(['Average percent fitting error: ' num2str(mean(AllFittingErrors))])
  
  disp('Average measured parameters:');
  disp('Position      Height       Width       Area');
   if NumRepeats>1,
     disp(num2str(MeanParameters(:,2:5)))
   else
     disp(num2str(AllFitResults(:,2:5))) 
   end % NumRepeats>1,
  RSDparameters(:,1)=ActualParameters(:,1); % Replace first column with peak numbers
  if NumRepeats>1,
    disp('Percent RSDs of measured parameters:');
    disp('Position      Height       Width       Area'); 
    disp(num2str(RSDparameters(:,2:5)))
    meanRSD=mean(RSDparameters);
    meanRSD=meanRSD(2:5);
  disp('Average %RSD of measured parameters for all peaks: '); disp(num2str(meanRSD))
  end % if NumRepeats>1
     if NumRepeats>1,
     PercentParameterErrors=100.*(MeanParameters-ActualParameters)./ActualParameters;
   else
     PercentParameterErrors=100.*(MeanParameters-ActualParameters)./ActualParameters; 
   end % NumRepeats>1,
  PercentParameterErrors(:,1)=ActualParameters(:,1); % Replace first column with peak numbers
  
  disp('Percent errors of measured parameters:');
  disp('Position      Height       Width       Area');
  disp(num2str(PercentParameterErrors(:,2:5)))
  AverageParameterErrors=mean(abs(PercentParameterErrors(:,2:5)));
  disp('Average Percent Parameter Error for all peaks: '); disp(num2str(AverageParameterErrors))