% Script to create a realistic model of an experimental multi-peak signal,
% based on the measured characteristics of an original signal.
clear Peaks
format compact
format short g
% Estimated from the original experimental data:  
x=.11:.101:450; % Independent variable of original data
Baseline=180+(290-180)./290.*x; % Slightly tilted but straight baseline 
noise=21;% Estimated from the random noise in original data.
% 
% Define P as a table of peaks in the original data, from findpeaksG,
% iPeak, ipf, etc.
%         Position      Height       Width         Area
P=[1       119.73       700.16       2.5385       2640.27
   2       206.31       3060.6       1.3842       4510.78
   3       214.78       1260.1       1.5837       2120.62
   4       252.48       16629        1.1265       19943 
   5       301.5        18945        1.0348       20871 
   6       360.2        504.17       1.0252       550.17 ];
% Create a model for each peak (Gaussian in this case, but could be any
% shape) based on the height, position, and width in the above matrix P.
% The width in the P table is ignored.
for peak=1:length(P) % peak=peak number (1 through 6 in this case)
    Height=P(peak,3);
    Position=P(peak,2);
    Width=P(peak,4);
    % Create a Gaussian for each peak and add to a row of matrix "Peaks"
    Peaks(peak,:)=Height.*gaussian(x,Position,Width); 
end            
y=sum(Peaks); % Add up all the separate peak into a multi-peak vector y
y=y+Baseline+noise.*randn(size(x)); % Add noise and baseine
clf
plot(x,y)

% Use this simulated signal to test the peak area measurement accuracy for
% a specified integration window width and baseline correction range. The
% small peaks give very poor area accuracy, both because the signal-to-
% noise ratio is poor but also because the baseline connection is less
% accurate. In this example, the baseline is estimated for the local region
% of each peak and is assumed to straight and tilted, which is accurate
% enough, because the limiting factor in this case is the random noise (as
% you can prove by setting noise=0 in line 9).
%
% Create a data point number integration range for each peak
ww=90; % Integration window width, expressed in number of data points
BasePoints=5; % Number of points at edges of window to count as baseline
disp('Peak area measurement by the peakfit function')
% Column heading for area comparison table:
disp('        peak #     Position    True area   Measured area   % error') 
for peak=1:length(P)
  TrueArea(peak)=trapz(x,Peaks(peak,:)); % Area without baseline or noise
  PeakPoint=val2ind(x,P(peak,2)); % Peak maxima expressed in data point
  prange=PeakPoint-ww/2:PeakPoint+ww/2; % Measurement range for each peak
  % Plot each peak over the measurement range in a separate figure window
  figure(peak+1) 
  clf
  plot(x(prange),y(prange))
  grid
  % Peak measurement for known peak shapes, using my peakfit.m function. 
  % https://terpconnect.umd.edu/~toh/spectrum/InteractivePeakFitter.htm
  NumPeaks=1; % Number of peaks to fit in each segment
  PeakShape=1; % Peak shape (1=Gaussian, 2-Lorentzian, 3=logistic distribution, etc.)
  BaselineMode=1; % Straight tilted baseline correction mode
  [FitResults(peak,:),GOF(peak,:)]=peakfit([x(prange);y(prange)],0,0,NumPeaks,PeakShape,0,0,0,BaselineMode);
  area(peak)=FitResults(peak,5); % Peak area is in the 5th column of FitResults
  title(['Peak ' num2str(peak)])
  PercentAreaError(peak)=100.*(TrueArea(peak)-area(peak))./TrueArea(peak);
  xlabel(['Area Error = ' num2str(PercentAreaError(peak)) '%' ])
  disp([P(peak,1)'  P(peak,2)'  TrueArea(peak) area(peak) PercentAreaError(peak)])
end
MeanAbsoluteError=mean(abs(PercentAreaError))
disp('Individual peak fits are displayed in Figures 2 through 7')
figure(1)
clf
plot(x,y)
xlabel('Independent variable: x') 
ylabel('Dependent variable, y')
title('Simulated signal, designed to match a real experimental signal')
% In this particular case, most of the error in area measurement is due to
% noise, not baseline estimation. You can prove this by setting noise=0 in
% line 9.