Home > wafo > misc > findpeaks.m

# findpeaks

## PURPOSE

find peaks of vector or matrix possibly rainflow filtered

## SYNOPSIS

[indx, indy] = findpeaks(S,Np,min_h,min_p)

## DESCRIPTION

``` FINDPEAKS find peaks of vector or matrix possibly rainflow filtered

CALL: [ix iy] = findpeaks(S,Np,min_h,min_p)

ix     = row subscripts to the peaks of S if iy present
otherwise linear index to peaks of S
iy     = column subscripts to the peaks of S
S      = matrix or vector
Np     = The Np highest peaks are found (if exist). (default 2)
min_h  = The threshold in the rainflowfilter (default 0.05*range(S(:))).
A zero value will return all the peaks of S.
min_p  = 0..1, Only the peaks that are higher than
min_p*max(max(S))  min_p*(the largest peak in S)
are returned (default  0).
Example:
x=(0:0.01:10); S=x.^2+10*sin(3*x)+0.5*sin(50*x); clf, plot(x,S)
ix = findpeaks(S',8,5,0.3) % Find highest 8 peaks that are not
% less that 0.3*"global max" and have
% rainflow amplitude larger than 5.

## CROSS-REFERENCE INFORMATION

This function calls:
 dat2tp Extracts turning points from data, range Calculates the difference between the maximum and minimum values. error Display message and abort function. ind2sub Multiple subscripts from linear index. sub2ind Linear index from multiple subscripts.
This function is called by:
 wafofig7 Joint distribution (pdf) of crest wavelength, Lc, and crest amplitude, Ac wafofig8 Joint distribution (pdf) of crest wavelength, Lc, and crest amplitude, Ac for extremal waves wspecplot Plot a spectral density

## SOURCE CODE

```001 function [indx, indy] = findpeaks(S,Np,min_h,min_p)
002 %FINDPEAKS find peaks of vector or matrix possibly rainflow filtered
003 %
004 % CALL: [ix iy] = findpeaks(S,Np,min_h,min_p)
005 %
006 %    ix     = row subscripts to the peaks of S if iy present
007 %             otherwise linear index to peaks of S
008 %    iy     = column subscripts to the peaks of S
009 %    S      = matrix or vector
010 %    Np     = The Np highest peaks are found (if exist). (default 2)
011 %    min_h  = The threshold in the rainflowfilter (default 0.05*range(S(:))).
012 %             A zero value will return all the peaks of S.
013 %    min_p  = 0..1, Only the peaks that are higher than
014 %                   min_p*max(max(S))  min_p*(the largest peak in S)
015 %                   are returned (default  0).
016 % Example:
017 %  x=(0:0.01:10); S=x.^2+10*sin(3*x)+0.5*sin(50*x); clf, plot(x,S)
018 %  ix = findpeaks(S',8,5,0.3) % Find highest 8 peaks that are not
019 %                             % less that 0.3*"global max" and have
020 %                             % rainflow amplitude larger than 5.
021 %
023
024
025 % Tested on: matlab 5.3
026 % History:
027 % revised pab Feb2004
028 % -changed default value for min_h
030 % revised pab 12.10.1999
031 %   fixed a bug
032 % by Per A. Brodtkorb  25.09.1999
033
034 error(nargchk(1,4,nargin))
035 if (nargin <4) | isempty(min_p)
036   min_p=0;
037 end
038 if (nargin <3) | isempty(min_h)
039   min_h = 0.05*range(S(:));
040 end
041
042 if (nargin <2) | isempty(Np)
043  Np = 2;
044 end
045 Ssiz=size(S);
046 if min(Ssiz)<2
047   ndim=1;
048   S=S(:)';
049   m=max(Ssiz);
050   n=1;
051 else
052   ndim=2;
053   n=Ssiz(1);
054   m=Ssiz(2);
055 end
056
057 x=(1:m).';
058
059 % Finding turningpoints of the spectrum
060 % Returning only those with rainflowcycle heights greater than h_min
061 indP=[]; % indices to peaks
062 ind=[];
063 for iy=1:n % find all peaks
064   TuP =dat2tp([x S(iy,:)'],min_h); % find turning points
065
066   if ~isempty(TuP)
067     ind=TuP(2:2:end,1); % extract indices to maxima only
068   else % did not find any , try maximum
069     [y ind]=max(S(iy,:));
070   end
071   if ndim>1,
072     switch iy,
073      case 1,ind2=find((S(iy,ind)>S(iy+1,ind)) );
074      case n,ind2=find((S(iy,ind)>S(iy-1,ind)));
075      otherwise
076       ind2=find((S(iy,ind)>S(iy-1,ind)) & (S(iy,ind)>S(iy+1,ind)));
077     end
078     if ~isempty(ind2)
079        indP=[indP;ind(ind2) iy*ones(length(ind2),1)];
080     end
081   end
082
083 end
084
085 if isempty(ind)|(isempty(indP)&(ndim>1))
086   indx=[];
087   indy=[];
088   return
089 end
090 if ndim>1
091   % indP(1) = col subscript, indP(2) = row subscript
092   ind=sub2ind([n m],indP(:,2),indP(:,1));
093 end
094
095 [Y ind2] = sort(S(ind));
096 ind2 = flipud(ind2(:)) ;
097
098 % keeping only the Np most significant peak frequencies.
099 ind=ind(ind2(1:min(Np,length(ind))));
100 if (min_p >0 ), % Keeping only peaks larger than min_p percent relative
101                 % to the maximum peak
102   [Y ind2]=find(S(ind) >min_p*max(S(ind)));
103   ind=ind(ind2);
104 end
105
106 if nargout>1
107   [indx indy]=ind2sub(Ssiz,ind);
108 else
109   indx=ind;
110 end
111```

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003