% Demonstration script for iPeak function. It generates a test signal
% consisting of several peaks, adds random noise, runs iPeak, then prints out 
% the actual peak positions, heights, widths and areas. 
% Each time you run it you get a different set of peaks. You can easily
% evaluate the accuracy of the measurements because the actual peak
% parameter values in this simulation are always integers.
% Must have ipeak in the path
%   T. C. O'Haver, February 2013, revised July 2019.
increment=1;
x=[1:increment:4000];

% For each simulated peak, compute the amplitude, position, and width
amp=round(5.*randn(1,38));  % Amplitudes of the peaks  (Change if desired)
pos=[200:100:3900];   % Positions of the peaks (Change if desired)
wid=65.*ones(size(pos));   % Widths of the peaks (Change if desired)

% Optionally add noise and background (Change if desired) 
Noise=.01; % Amount of random noise added to the signal. (Change if desired) 
background=0; % Strength of bnackground added to the signal. (Change if desired)

% A = matrix containing one of the unit-amplidude peak in each of its rows
A = zeros(length(pos),length(x));
ActualPeaks=[0 0 0 0 0];
p=1;
for k=1:length(pos)
  if amp(k)>.3  % Keep only those peaks above a certain amplitude
      % Create a series of peaks of different x-positions
      A(k,:)=gaussian(x,pos(k),wid(k)); % Gaussian peaks
      % A(k,:)=rectanglepeak(x,pos(k),wid(k)); % rectangle peaks
      % A(k,:)=ones(size(x))./(1+((x-pos(k))./(0.5.*wid(k))).^2);  % Lorentzian peaks
      % Assembles actual parameters into ActualPeaks matrix: each row = 1
      % peak; columns are Peak #, Position, Height, Width, Area
      ActualPeaks(p,:) = [p pos(k) amp(k) wid(k) 1.065.*amp(k).*wid(k)]; 
      p=p+1;
  end 
end
z=amp*A;  % Multiplies each row by the corresponding amplitude and adds them up
y=z+Noise.*randn(size(z));  % Optionally adds random noise
y=y+background.*gaussian(x,-2000,4000); % Optionally adds a broad background signal
demodata=[x' y']; % Assembles x and y vectors into data matrix

% Now call iPeak, with specified values of AmpT, SlopeT, SmoothW, and FitW.
% (You can change theses values if desired).
 EstimatedPeaks=ipeak(demodata,0,0.5,0.00005,9,27,2499.5,999,0);

% Compare MeasuredPeaks to ActualPeaks
format short g
disp('-----------------------------------------------------------------------------')
disp('         Peak #    Position      Height       Width        Area        ')
disp(ActualPeaks)
disp('         Peak #    Position      Height       Width        Area        FitError')

disp(EstimatedPeaks)
try
    disp('Peak Area Measurement (Shift-U key)')
disp('          Peak      Position    PeakMax    Peak-valley    Perp drop   Tan skim')
M=measurepeaks(x,y,0.00005,0.5,9,27,1);
disp(M)
ActualAreas=ActualPeaks(:,5);
MeasuredAreas=M(:,5);
PerpDropErrors=100*(ActualAreas-MeasuredAreas)./ActualAreas;
AverageAbsPercentError=mean(abs(PerpDropErrors))
ipeak(demodata,0,0.5,0.00005,9,27,2499.5,999,0);
catch me
end