Re: Finding peaks/max in a graph

[bertrand roessli]

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
-------------------------------------------------------------