Code covered by the BSD License

Peak finding and measurement

Tom O'Haver (view profile)

20 Jul 2006 (Updated )

Function to locate and measure the positive peaks and valleys in noisy data sets.

ipeakdemo1
```function ipeakdemo1
% Demonstration script for iPeak function. It generates a test signal
% consisting of several peaks, adds random noise, runs ipeak, then prints out
% the actual and the measured 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.
%   T. C. O'Haver, September 2011
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=40.*ones(size(pos));   % Widths of the peaks (Change if desired)
Noise=0.5; % 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)>.2,  % Keep only those peaks above a certain amplitude
% Create a series of peaks of different x-positions
A(k,:)=exp(-((x-pos(k))./(0.6005615.*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
y=z+Noise.*randn(size(z));  % Optionally adds random noise
demodata=[x' y']; % Assembles x and y vectors into data matrix

% Initial values of variable peak detection 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.05*max(y);
SmoothWidth=round(WidthPoints/2)+2;  % SmoothWidth should be roughly equal to 1/2 the peak width (in points)
FitWidth=round(WidthPoints); % FitWidth should be roughly equal to 1/2 the peak widths(in points)

% Now call iPeak, with specified values of AmpT, SlopeT, SmoothW, and FitW,
% with AUTOZERO=1 (ON) and the top window showing the first peak. Once
% iPeak is running, you can change any of these using
% keystroke commands (press "K" to see a list of commands).
MeasuredPeaks=ipeak(demodata,0,AmpThreshold,SlopeThreshold,SmoothWidth,FitWidth,ActualPeaks(1,2),120,1);

% Compare MeasuredPeaks to ActualPeaks
disp('-----------------------------------------------------------------')
disp(['Signal to noise ratio of smallest peak = ' num2str(min(ActualPeaks(:,3))./Noise)])
disp('         Peak #    Position      Height       Width        Area')
ActualPeaks
MeasuredPeaks(:,1:5)

SizeResults=size(MeasuredPeaks);
merror=zeros(SizeResults(1),5);
for n=1:SizeResults(1),
indexn=val2ind(ActualPeaks(1:5,2),MeasuredPeaks(n,2));
merror(n,:)=100.*abs(ActualPeaks(indexn,:)-MeasuredPeaks(n,1:5))./ActualPeaks(indexn,:);
merror(n,1)=MeasuredPeaks(n,1);
% MeasuredPeaks(n,1)=indexn;
end
AveragePercentError=mean(merror)
disp('Demonstration of autozero background correction, for separated, narrow ')
disp('peaks on a large baseline. Hint: Turn on the Autozero mode (T key),')
disp('adjust the zoom setting so that the peaks are shown one at a time in')
disp('the upper window, then press the P key to display the peak table.')

disp('End of demo.')
% ----------------------------------------------------------------------
function g = gaussian(x,pos,wid)
%  gaussian(X,pos,wid) = gaussian peak centered on pos, half-width=wid
%  X may be scalar, vector, or matrix, pos and wid both scalar
% Examples: gaussian([0 1 2],1,2) gives result [0.5000    1.0000    0.5000]
% plot(gaussian([1:100],50,20)) displays gaussian band centered at 50 with width 20.
g = exp(-((x-pos)./(0.6005615.*wid)).^2);

```