[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
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).
% Author: EF <address@hidden>
% 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.
% Author: EF <address@hidden>
% 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.
% Author: EF <address@hidden>
% 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.
% Author: EF <address@hidden>
% Description: make row
% E.Farhi. 12/95 (address@hidden)
% 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
-------------------------------------------------------------