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.
 
  See also  dat2tp

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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 %
022 % See also  dat2tp
023 
024 
025 % Tested on: matlab 5.3
026 % History:
027 % revised pab Feb2004
028 % -changed default value for min_h
029 % -added nargchk  
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