help-octave
[Top][All Lists]

## Re: Finding peaks/max in a graph

 From: bertrand roessli Subject: Re: Finding peaks/max in a graph Date: Mon, 5 Apr 2004 12:26:56 +0200 (CEST)

```You might try the following set of functions.

Hope it helps (and works),

BR

function [peak,pmax,sy,ny,sdy,ndy,sddy,nddy,pmin,pmaxsav] = findpeaks(x,ys,thr,
pos,width,last)
% findpeaks : Quick search of peaks in a signal.
%Syntax: [peakmat,pmax,sy,ny,sdy,ndy,sddy,nddy,pmin] =
findpeaks(x,y,{threshold=1, pos,width,mode ='silent'}) or findpeaks(y)
%
% This function searches peaks into an (x,y) signal, using a specified
% 'threshold' in noise units, usually around 1 ; few peaks for a high value.
% 'peakmat' is an array of rows for each maximum identified
%          [ index left_width right_width  max_pos Intensity Width ].
% 'pmax' is a Dirac type peaks signal.
% Averaged signals (including derivatives) and noises are also returned.
% opt : pos, width = part of 'y' to be scanned for peaks around pos.
%       mode is 'silent', 'normal' or 'verbose' (do plot).

% Description: quick search of peaks in a signal.

% Part of 'Spectral tools'. E.Farhi. 07/96 rev 12/97
% uses : diff.m
%        smooth.m
%        interp.m
%        vect2row.m

if ~exist('last')
last = 'normal';
end

if isempty(last)
last = 'normal';
end

if strcmp(last,'silent')
tmp = 0;
elseif strcmp(last,'normal')
tmp = 1;
else
tmp = 2;
end

if (nargin < 1)
error('usage : [peaks,pmax,sy,ny,sdy,ndy,sddy,nddy] =
findpeaks(x,y,{threshold=1}) or findpeaks(y)');
end

if (tmp>0)
fprintf(1,'Peak Searching\n');
end

if (nargin == 1)
ys=x;
x=1:length(x);
end

if ~exist('thr')
thr = [];
end
if isempty(thr)
thr = 1;
end

if ~exist('pos') | ~exist('width')
pos = [];
width = [];
end
if isempty(pos) |  isempty(width)
nd = length(x);
pos = x(ceil(nd/2));
width = abs(x(nd)-x(1));
end

if thr == 0, thr = 0.0001; end

y=ys(:);

nd = length(y);
dx = x(2) - x(1);

%if any(diff(x) ~= dx) % step is not constant, need to interpolate
%       nx = linspace(x(1), x(nd), round((x(nd)-x(1))/dx));
%       y = interp1(x,y,nx);
%       x=nx;
%       plot(x,y)
%end

t=find( (x>=(pos-width)) & (x<=(pos+width)) );          % index zone around pos
t = t(:);
y = y(t); y=y(:);
x = x(t); x=x(:);
nd = length(y);

sy = smooth (y);
dy = diff (sy); dy(nd) = dy(nd-1);
sdy = smooth (dy);
ddy = diff (sdy); ddy(nd) = ddy(nd-1);
sddy = smooth (ddy);

t=(abs(y) <= median(abs(y)));           % noises computation
t = t(:);
nt = sum(t);
ny = std(t.*y - t.*sy);
t=(abs(dy) <= median(abs(dy)));
ndy = std(t.*dy - t.*sdy);
t=(abs(ddy) <= median(abs(ddy)));
nddy = std(t.*ddy - t.*sddy);

t = dy(1:(nd-1)).*dy(2:nd);     % for 1st derivative sign changes.
t(nd) = sdy(nd)^2;
t = t(:);

pmax = (abs(dy) <= ndy/thr) | (t <= (ndy/thr)^2);       % extrema

pmin = ( pmax & (sddy > nddy*thr));
pmax = ( pmax & (sddy < -nddy*thr));

% now alterning maxima and minima ...

pmaxsav = pmax;
peaklist=find(pmax);

if (tmp > 0)
fprintf(1,'* analysing extrema (2nd order interpolation), %i possible
peaks.',length(peaklist));
fprintf(1,'\nnumber        indexes     max_pos     width   intensity
(exp int)');
end

index =1;               % index in original peak list (from pmax)
ip = 1;         % index in final peak table
peak = [ 0 0 0 0 0 ];
pmax = 0*y;

while (index <= length(peaklist))        % scanning all peaks in pmax
peakindex = peaklist(index);
t = find(pmin(1:peakindex));    % search min before
if (~isempty(t))
minbefore = t(length(t));
else
minbefore = 1;
end
t = find(pmin(peakindex:nd));   % search min after
if (~isempty(t))
minafter = t(1) + peakindex -1;
else
minafter = nd;
end
t = minbefore:minafter;
[dummy,t] = max(y(t));  % find real max beween 2 min.
if (isempty(t))
index = index+1;                        % no max found
% trying next in pmax
else
max_pos = t(1) + minbefore -1;                  % max found
t = max(y(minbefore),y(minafter));
if ((y(max_pos) - t) <= ny)     % this is a noise peak.
if (y(minbefore) > y(minafter)) % equivalent minima
pmin(minbefore) = 0;
else                    % higher minimum removed
pmin(minafter) = 0;
end
if (minafter == nd)
break;                  % finish
else
if (minbefore == 1)
index = index +1;
end  % trying again

end
else                            % real peak
if (~isempty(find((peak(:,1) - max_pos) == 0)))
index = index + 1;              % already found
else
peak(ip,1) = max_pos;
t = min(y(minbefore),y(minafter));
t = 0.1*y(max_pos) + 0.9*t;     % 10 percent hight for peak
halfafter = find(y(max_pos:minafter) < t);
if (isempty(halfafter))
halfafter = minafter;
else
halfafter = halfafter(1) + max_pos -1;
end
halfbefore = find(y(minbefore:max_pos) < t);
if (isempty(halfbefore))
halfbefore = minbefore;
else
halfbefore = halfbefore(length(halfbefore)) + minbefore
-1 ;
end
peak(ip, 2) = max_pos - halfbefore;
peak(ip, 3) = halfafter - max_pos;
t = halfbefore:halfafter;
%disp([ max_pos length(t) ])
% might be done with polyfit but this is better...
[ smoothed, p] =
interpsp(x(t),y(t),x(max_pos),max(ceil(length(t)/2),3),2);
%               p = polyfit(x(t),y(t),2);

S = x(max_pos);
W = abs(x(halfafter) - x(halfbefore));
In = y(max_pos);

if (length(p) == 3)
a=p(1); b=p(2); c=p(3); % 2nd order interpolation
delta = (b*b - 4*a*c);
if ((a < 0) & (delta >0))
S=-b/2/a;
W=abs(sqrt(delta)/a);
In=a*S*S+b*S+c;
else
t = minbefore:minafter;
[ smoothed, p] =
interpsp(x(t),y(t),x(max_pos),max(ceil(length(t)/2),3),2);
%                               p = polyfit(x(t),y(t),2);

if (length(p) == 3)
a=p(1); b=p(2); c=p(3); % 2nd order
interpolation
delta = (b*b - 4*a*c);
delta = -1;
if ((a < 0) & (delta >0))
S=-b/2/a;
W=abs(sqrt(delta)/a);
In=a*S*S+b*S+c;
else
p = [];
end
end
end
else                    % interpolation not performed
p = [];
end
end
if S > max(x) | S < min(x) | W > abs(max(x)-min(x))
p = [];
end
if isempty(p)                   % interpolation not performed
S = x(max_pos);
W = abs(x(halfafter) - x(halfbefore));
In = y(max_pos);
end

t = [ minbefore minafter ];
wx = find(x > S+W);
if ~isempty(wx), t = [ t wx(1)]; end
wx = find(x >= S-W);
if ~isempty(wx), t = [ t wx(1) ]; end
t = abs(t - max_pos) * 3;
t = max(t);
t = (max_pos -t):(max_pos+t);
t = t(find(t>=1 & t<= nd));
t = min(y(t)); % compared with nearer min peak
peak(ip,4) = S;
peak(ip,5) = In - t;
peak(ip,6) = W;
if (tmp > 0)
fprintf(1,'\nPeak %3i : (%3i-%3i-%3i) S=%7.2f W=%7.2f
In=%7.2f (%7.2f)', ip, peak(ip, 2), peak(ip, 1), peak(ip, 3), S, W, In,
y(max_pos));
end
ip = ip + 1;
index = index + 1;
pmax(max_pos) = 1;
end % from if (~isempty(find((peak(:,1) - max_pos) == 0)))
end
end
end
if tmp
fprintf(1,'\n');
end

if (tmp > 1)
fprintf(1,'\n%i peaks identified\n', ip-1);
clf;
title('Peaks position in signal');
plot(x,y,x,pmax.*y);
end

if (size(ys,2)==1)
sy = sy';
pmax = pmax';
sdy = sdy';
sddy = sddy';
end

function [y, p] = interpsp(xd,yds, x, nb_of_pts, order, last)
% interpsp : Lagrange 1D interpolation/smoothing.
%Syntax: [itrp, p] = interpsp(x,y, {new_x=x, nb_of_pts, order=1,mode='silent'})
or interpsp(y)
%
% This function returns the interpolation/smoothing of (x,y) by Lagrange
% method at 'new_x' values (new_x=x if not precised).
% optional 'nb_of_pts': number of points used for each interpolation computation
% optional 'order': order of polynomial interpolation with nb_of_points
% 'itrp' is the interpolated data. 'p' is last computed polynome.

% Description: Lagrange 1D interpolation/smoothing.

% Part of 'Spectral tools'. E.Farhi. 07/96 rev 11/97
% uses : vect2row.m

% Argument processing -----------------------------------------------

if ~exist('last')
last = 'silent';
end

if isempty(last)
last = 'silent';
end

if strcmp(last,'silent')
tmp = 0;
else
tmp = 1;
end

if (tmp>0)
fprintf(1,'Smoothing/Interpolation\n');
end

if (nargin < 1)
error('usage: [interpolated,p] = interpsp(x,y, {new_x=x, nb_of_pts,
order=1}) or interp(y)');
end

if (nargin == 1)
yd=xd;
xd=1:length(yd);
end

xd = vect2row(xd);
yd = vect2row(yds);

xdmax = max(abs(xd));
nd=length(xd);

if ~exist('order')
order = [];
end
if isempty(order)
order = 1;
end

if ~exist('nb_of_pts')
nb_of_pts = [];
end
if isempty(nb_of_pts)
nb_of_pts = order+1;
end

if ~exist('x')
x = [];
end
if isempty(x)
x = xd;                 % initial X axis
if (tmp>0)
disp('new_x = x');
end
end

if (xd(1) > xd(nd))
xd = xd(nd:-1:1);
yd = yd(nd:-1:1);
x = x(length(x):-1:1);
xddec = 1;
else
xddec = 0;
end

xd = xd  / xdmax;       % xd set to 0:1
x = vect2row(x);
x = x  / xdmax;

if (length(yd) ~= nd)
error('x and y vectors must have same number of elements.');
end;
nb_of_pts = min(nb_of_pts, nd);
order = min(order, nb_of_pts-1);
if (nb_of_pts < 2)
p = [];
y = yds;
return;
end

if (tmp>0)
fprintf(1,'Nb of points %i, Order %i. ',nb_of_pts,order);
end
nb_of_pts = nb_of_pts -1;
fn = floor(nb_of_pts/2);
cn = ceil(nb_of_pts/2);

lx = length(x);

if (tmp>0)
fprintf(1,'Running ');
end
t0 = clock;
X = eye(order+1);
Y = zeros(1,order+1);

if (lx == nd)
if(x == xd)                   % computing nearest indexes in xd yd to be used
id_tab = 1:nd;                                  % no change on
x axis.
else
dif = x-xd;
if all( dif == dif(1) )                 % x is just xd shifted by
constant value
step = x(2) - x(1);
id_tab = (1:nd) + round( step/dif(1) );
end
end
end
if (~exist('id_tab'))                           % already computed ?
x = sort(x);
xd = sort(xd);
yds = sort(yds);
for index = 1:lx                        % general case.
[dummy,id] = min(abs(x(index) - xd ));  % closest index in x
if (length(id))
id_tab(index) = id(1);
else
id_tab(index) = nd;
end
end
end

% ----------------------- main iteration loop -----------------------

if (order >= 1)

if (tmp>0)
fprintf(1,'(10 %%) #');
end

idprec = 0;
for index = 1:lx                        % index in new_x
id = id_tab(index);             % corresponding index in x data.
id = max(1, min(id, nd));

if (id ~= idprec)               % need to re-compute polynome.
idprec = id;
if (xd(id) < x(index))
id = max(1, id-fn):min(nd,id+cn);       % indexes in
xd,yd for extracting data
else
id = max(1, id-cn):min(nd,id+fn);
end
if id == 1
id = 1:min(nb_of_pts+1,nd);
elseif id == nd
id = max(1,nd - nb_of_pts):nd;
end
xs = xd(id);
ys = yd(id);

% ------------------- smoothing/interpolation -----------------------

for j=1:(order+1)
X(j,j) = sum(xs.^(2*j-2));
Y(j) = sum(ys.*(xs.^(j-1)));
for k=(j+1):(order+1)
X(k,j) = sum(xs.^(k+j-2));
X(j,k) = X(k,j);
end;
end;

A = X\Y';                       % Least Square Approx.

end;

% ----------------------- evaluates polynome -------------------------

k=0;
for j=1:order+1
k = k +A(j)*x(index)^(j-1);
end;
y(index) = k;

if (rem(index,ceil(lx/10)) == 0)
if (tmp>0)
fprintf(1,'#');
end
end
end

else                    % order = 0 -> averaging signal ...
ys = zeros(1,nd+nb_of_pts);
xs = ones(1,nd+nb_of_pts)*nb_of_pts;
for index = 1:nb_of_pts
ys(index:(nd+index-1)) = ys(index:(nd+index-1)) + yd(1:nd);
xs(index) = index;
xs(nd+index) = nb_of_pts-index+1;
if (tmp>0)
fprintf(1,'#');
end
end
y = ys ./ xs;
y = y(cn:nd+cn);
y = y(id_tab);
A = y(length(y));
end

if (tmp>0)
fprintf(1,' - OK %5.1f s.\n',etime (clock, t0));
end

if (nargout >= 2)       % for standard MatLab/Octave notation.
X = 1 / xdmax;
for index = 1:order+1
A(index) = A(index) * X^(index-1);
end
p = fliplr(vect2row(A));
if (tmp>0)
disp('Interpolation polynome (Last) :');
disp(p);
end;
end

if (size(yds,2)==1)
y = y';
end

if (xddec == 1)
y = y(length(y):-1:1);
end

function [ny,c] = smooth(yd,N,M)
% smooth : Data smoothing by Savitzky-Golay method
%Syntax: [smoothed_y, coefs] = smooth(y,{N,M})
%
% Smoothes the y signal by M-th order Savitzky-Golay method with N points.
% This algorithm assumes that corresponding x axis is evenly spaced.

% Description:  data smoothing by Savitzky-Golay method

% Part of 'Spectral tools'. E.Farhi. 07/96
% From : Numerical recipes in C. p 650

% Argument processing -----------------------------------------------
% uses : vect2row.m

if (nargin < 1)
error('usage: [smoothed_y, coefs] = smooth(y,{N,M=2})');
end

nd=length(yd);

if (nargin < 3)
M = 2;
end

if (nargin <= 1)
N = ceil(max(2*M,nd/50));
end

if (nd <= 2*N)
error('not enough points in data');
end

%  Savitzky-Golay coefficients -----------------------------------------------

N = ceil(N/2);

A = zeros (2*N +1, M+1);
c = zeros(1,2*N+1);

n=(-N):N;
for j=0:M
A(:,j+1) = n'.^j;           % Aij = i^j
end

B = pinv(A'*A);
B = B(1,1:(M+1));

for n=1:(2*N+1)
c(n) = A(n,:) * B';   % these are Savitzky-Golay coefficients
end

% Smoothing ------------------------------------------------------------------

ny = 0*yd;

for n=1:(2*N+1)
ny((N+1):(nd-N)) = ny((N+1):(nd-N)) + c(n) * yd(n:(nd-2*N-1+n));
end
for n=1:N
ny(n)=yd(n);
ny(nd-n+1) = yd(nd-n+1);
end

function v=vect2row(vin)
% vect2row : Make row
%Syntax: V=vect2row(Vin)
%
% 'V' is 'Vin' reshaped as a row vector.
% For matrix input, output is as Nrows < Ncolumns.

% Description: make row

% Spectral tools.

[nr,nc] = size (vin);

if ( nr > nc)
v=vin';
else
v=vin;
end

%end

On Mon, 5 Apr 2004, edA-qa mort-ora-y wrote:

> I need to find the peaks in a graph (sampled data).  By peaks I mean
> each point at which the slope changes from positive to negative (or 0
> slope, but these are usually cusps).
>
> As this is sampled data (resulting from an FFT actually), I need to find
> only a certain number of the peaks (the maximum peaks, since the noise
> in the sample yields a high number of real peaks).
>
> I wrote a function which does this, but it is extremely slow (I'm
> working with sample sets of 100K in size).  Does anybody else know how I
> can achieve this any faster, or which standard functions.
>
> Purpose: To find the most predominant frequencies in a waveform (sample
> fed into FFT, then into peaks finding algorithm)
>
> Attached: The script I wrote
>
>

--

Dr. Bertrand Roessli
Laboratory for Neutron Scattering
ETHZ and Paul-Scherrer Institute
CH-5232 Villigen PSI

Tel.:+41 56 310 44 01
Fax :+41 56 310 29 39

-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------

```