Re: Finding peaks/max in a graph

[David Bateman]

A for-loop is not the way to go about this. Consider the code fragment

n = 9
a = [1 2 3 4 5 4 3 2 1]';
peaks = find([a(2:n,1) - a(1:n-1,1) < 0; 1] & [1; a(1:n-1,1) - a(2:n,1) < 0]);

that will find all indexes of the peaks in the data in a single statement

D.

Daprès edA-qa mort-ora-y <address@hidden> (le 05/04/2004):
> 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
> 
> -- 
> edA-qa mort-ora-y
> Idea Architect
> http://disemia.com/
> 

> function p = find_peaks( data, num )
> 
>       cur = data(1);
>       mn = cur;
>       dir = 0;
> 
>       peaks = [];
> 
>       for i = 2:length(data)
>               s = sign( data(i)-cur );
>               if( s != dir )
>                       dir = s;
> 
>                       if( dir < 0 )
>                               peaks   = [peaks; i, data(i)];
>                       endif
>               endif
> 
>               cur = data(i);
>               mn = min(cur,mn);
>       endfor
> 
>       l = rows(peaks);
>       if( num > l )
>               num = l;
>       endif
> 
>       # Ensure we are using a matrix, otherwise sort treats it as a vector
>       # we should populate this with the minimum value actually... (TBD)
>       if( l == 1 )
>               peaks = [peaks; 0, mn];
>               l = 2;
>       endif
> 
>       [ spv, spi ] = sort( peaks );
> 
>       p = [];
>       
>       for i = 0:(num-1)
>               p = [p; peaks( spi(l-i,2), 1 ), peaks( spi(l-i,2), 2 ) ];
>       endfor
>       
> endfunction

-- 
David Bateman                                address@hidden
Motorola CRM                                 +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax) 
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary



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