% A simple self-contained demonstration of the findpeaksnr function (line 54)
% applied to noisy synthetic data set consisting of a random number of narrow 
% peaks.  Each time you run this, a different set of peaks is generated.
% Calls the fundpeaksnr nrfunction, which must be in the Matlab path.
% See http://terpconnect.umd.edu/~toh/spectrum/Smoothing.html and 
% http://terpconnect.umd.edu/~toh/spectrum/PeakFindingandMeasurement.htm
% Tom O'Haver (toh@umd.edu). Version 3.1 February 2013
% You can change the signal characteristics in lines 9-16
format short g
format compact
increment=2;
x=[1:increment:8000];

% For each simulated peak, compute the amplitude, position, and width
pos=[200:200:7800];   % Positions of the peaks (Change if desired)
amp=round(5.*randn(1,length(pos)));  % Amplitudes of the peaks  (Change if desired)
wid=50.*ones(size(pos));   % Widths of the peaks (Change if desired)
Noise=.2; % Amount of random noise 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)>.1,  % Keep only those peaks above a certain amplitude
      % Create a series of peaks of different x-positions
      A(k,:)=exp(-((x-pos(k))./(0.6005612.*wid(k))).^2); % Gaussian 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.0646.*amp(k)*wid(k)]; 
      p=p+1;
  end; 
end 
z=amp*A;  % Multiplies each row by the corresponding amplitude and adds them up
sqrtnoise=sqrt(abs(z)).*randn(size(z));
propnoise=abs(z).*randn(size(z));
sz=std(sqrtnoise);
pz=std(sqrtnoise);
nz=Noise.*sqrtnoise./sz;
y=z+nz;  % Adds square-root signal-dependent random noise
% y=y+.1.*randn(size(z));  % Adds constant random noise
% y=y+5.*gaussian(x,0,4000); % Optionally adds a broad background signal
% y=y+gaussian(x,0,4000); % Optionally adds a broad background signal
% demodata=[x' y']; % Assembles x and y vectors into data matrix

figure(1);plot(x,y,'r')  % Graph the signal in red
title('Detected peaks are numbered. Peak table is printed in Command Window')

% Initial values of variable parameters
WidthPoints=mean(wid)/increment; % Average number of points in half-width of peaks
SlopeThreshold=0.5*WidthPoints^-2; % Formula for estimating value of SlopeThreshold
AmpThreshold=0.02*max(y);
SmoothWidth=round(WidthPoints);  % SmoothWidth should be roughly equal the peak width (in points)
FitWidth=round(WidthPoints); % FitWidth should be roughly equal to the peak widths (in points)

% Lavel the x-axis with the parameter values
xlabel(['SlopeThresh. = ' num2str(SlopeThreshold) '    AmpThresh. = ' num2str(AmpThreshold) '    SmoothWidth = ' num2str(SmoothWidth) '    FitWidth = ' num2str(FitWidth) ])

% Find the peaks
tic;
Measuredpeaks=findpeaksnr(x,y,SlopeThreshold,AmpThreshold,SmoothWidth,FitWidth,3);
ElapsedTime=toc;
PeaksPerSecond=length(Measuredpeaks)/ElapsedTime;

% Display results
disp('---------------------------------------------------------')
disp(['SlopeThreshold = ' num2str(SlopeThreshold) ] )
disp(['AmpThreshold = ' num2str(AmpThreshold) ] )
disp(['SmoothWidth = ' num2str(SmoothWidth) ] )
disp(['FitWidth = ' num2str(FitWidth) ] )
disp(['Speed = ' num2str(round(PeaksPerSecond)) ' Peaks Per Second' ] )
disp('          Peak #     Position      Height      Width        SNR')
Measuredpeaks  % Display table of peaks
figure(1);text(Measuredpeaks(:, 2),Measuredpeaks(:, 3),num2str(Measuredpeaks(:,1)))  % Number the peaks found on the graph
clear PercentErrors
if length(ActualPeaks)==length(Measuredpeaks),
    PercentErrors(:,1:4)=100.*(ActualPeaks(:,1:4)-Measuredpeaks(:,1:4))./ActualPeaks(:,1:4);
    PercentErrors(:,1)=Measuredpeaks(:,1)
    AverageAbsolutePercentErrors=mean(abs(100.*(ActualPeaks(:,1:4)-Measuredpeaks(:,1:4))./ActualPeaks(:,1:4)))
end