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
Extracts turning points from data, | |
Calculates the difference between the maximum and minimum values. | |
Display message and abort function. | |
Multiple subscripts from linear index. | |
Linear index from multiple subscripts. |
Joint distribution (pdf) of crest wavelength, Lc, and crest amplitude, Ac | |
Joint distribution (pdf) of crest wavelength, Lc, and crest amplitude, Ac for extremal waves | |
Plot a spectral density |
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
Comments or corrections to the WAFO group