MM2LC Extracts level-crossing spectrum from min2Max cycles. CALL: [lc, alpha] = mm2lc(mM,def,plotflag,sa); lc = two column matrix with levels and number of upcrossings, i.e., level-crossing spectrum. alpha = irregularity factor, approximately Tmaxima/Tz mM = min2Max cycles (possibly rainflow filtered). def = 1, only upcrossings. 2, upcrossings and maxima (default). 3, upcrossings, minima, and maxima. 4, upcrossings and minima. plotflag = 0, no plotting 1, plot the number of upcrossings overplotted with Rice formula for the crossing intensity for a Gaussian process (default). sa = standard deviation of the process (Default estimates it from the number of upcrossings) Example: x = load('sea.dat'); tp = dat2tp(x); mM = tp2mm(tp); lc = mm2lc(mM); See also lcplot
Plots level-crossing spectrum (lc) | |
Clear variables and functions from memory. | |
Difference and approximate derivative. | |
Display message and abort function. | |
Set unique. |
001 function [lc,alpha] = mm2lc(mm,def,plotflag,sa) 002 %MM2LC Extracts level-crossing spectrum from min2Max cycles. 003 % 004 % CALL: [lc, alpha] = mm2lc(mM,def,plotflag,sa); 005 % 006 % lc = two column matrix with levels and number of upcrossings, 007 % i.e., level-crossing spectrum. 008 % alpha = irregularity factor, approximately Tmaxima/Tz 009 % mM = min2Max cycles (possibly rainflow filtered). 010 % 011 % def = 1, only upcrossings. 012 % 2, upcrossings and maxima (default). 013 % 3, upcrossings, minima, and maxima. 014 % 4, upcrossings and minima. 015 % 016 %plotflag = 0, no plotting 017 % 1, plot the number of upcrossings overplotted 018 % with Rice formula for the crossing intensity 019 % for a Gaussian process (default). 020 % 021 % sa = standard deviation of the process 022 % (Default estimates it from the number of upcrossings) 023 % 024 % Example: 025 % x = load('sea.dat'); tp = dat2tp(x); mM = tp2mm(tp); 026 % lc = mm2lc(mM); 027 % 028 % See also lcplot 029 030 031 % Tested on: matlab 5.3 032 % History: 033 % revised pab Feb2004 034 % - added alpha 035 % revised pab 25.04.2001 036 % -speeded up the for loop even further. 037 % revised pab 19.04.2001 038 % - fixed a bug: a forgotten transpose. 039 % revised pab 30.12.2000 040 % - vectorized the for loop to speed up things 041 % revised by pab 11.08.99 042 % changed name from mm2cross to mm2lc 043 % revised by Per A. Brodtkorb 01.10.98 044 % added: overplot the crossingspectrum with Rice formula for crossing 045 % intensity for a Gaussian process 046 047 048 049 error(nargchk(1,4,nargin)) 050 051 % Default values 052 if nargin<4, sa=[]; end % unknown stdev is default 053 if nargin<3|isempty(plotflag), plotflag=1; end % default plot final result 054 if nargin<2|isempty(def), def=2; end % default upcrossings & maxima 055 056 057 if ((def<1) | (def>4)) 058 error('def must be one of (1,2,3,4).') 059 end 060 061 index = find(mm(:,1) <= mm(:,2)); 062 063 if isempty(index) 064 error('Error in input mM.') 065 end 066 067 cc = mm(index,:); clear index 068 ncc = length(cc); 069 070 minima = [cc(:,1) ones(ncc,1) zeros(ncc,1) ones(ncc,1)]; 071 maxima = [cc(:,2) -ones(ncc,1) ones(ncc,1) zeros(ncc,1)]; 072 073 extremes = [maxima; minima]; 074 [tmp, ind] = sort(extremes(:,1)); 075 extremes = extremes(ind,:); 076 077 if 1, 078 % pab 30.12.2000 079 % indices to matching entries. 080 ind = (diff(tmp) == 0); 081 % Create position mapping vectors 082 tmp = [1;~ind]; 083 iy = cumsum(tmp); 084 ix = [find(~ind);2*ncc]; 085 086 ind = find(ind).'; % added transpose (pab 19.04.2001) 087 % Alternative call for finding ix: 088 %ix = (1:2*ncc).'; 089 %ix(ind) = []; 090 else 091 %Alternatively, ix,iy and ind may be found by the following: (slower) 092 [tmp, ix, iy] = unique(extremes(:,1)); 093 ind = find(diff(iy)==0)'; 094 end 095 096 clear tmp 097 extr = extremes(ix,:); % Keep only unique crossing levels 098 nx = size(extr,1); 099 if any(ind) 100 if 1, 101 % pab 25.04.2001 speeded up the for loop: 102 l = diff([0; ix]); % run lengths (ie, number of crossings) 103 ind1 = find([1 diff(ind)>1]); 104 for iz = ind1 105 jy1 = ind(iz); 106 jx = iy(jy1); 107 jy2 = jy1+l(jx)-2; 108 extr(jx,2:4) = extr(jx,2:4) + sum(extremes(jy1:jy2,2:4),1); 109 end 110 else 111 % Old call: kept just in case (slow) 112 for iz = ind 113 extr(iy(iz),2:4) = extr(iy(iz),2:4) + extremes(iz,2:4); 114 end 115 end 116 end 117 clear extremes 118 119 switch def 120 case 1,% Only upcrossings 121 lc=[extr(1:nx,1) cumsum(extr(1:nx,2)) - extr(1:nx,4)]; 122 case 2,% Upcrossings + maxima 123 lc=[extr(1:nx,1) cumsum(extr(1:nx,2)) + extr(1:nx,3)-extr(1:nx,4)]; 124 case 3,% Upcrossings + minima + maxima 125 lc=[extr(1:nx,1) cumsum(extr(1:nx,2)) + extr(1:nx,3)]; 126 case 4,% Upcrossings + minima 127 lc=[extr(1:nx,1) cumsum(extr(1:nx,2))]; 128 lc(nx,2)=lc(nx-1,2); 129 end 130 131 %% Plots are made by lcplot 132 if plotflag 133 h=lcplot(lc,2,0,sa); 134 end 135 if nargout>1 136 cmax = max(lc(:,2)); 137 alpha = cmax/ncc;% approximately Tmaxima/Tz 138 end 139 return 140 141 142 143 144 % Old call: slow (kept just in case) 145 ii=1; 146 n=length(extremes); 147 extr=zeros(n,4); 148 extr(1,:)=extremes(1,:); 149 for i=2:n 150 if extremes(i,1)==extr(ii,1); 151 extr(ii,2:4) = extr(ii,2:4)+extremes(i,2:4); 152 else 153 ii=ii+1; 154 extr(ii,:) = extremes(i,:); 155 end 156 end 157 [xx nx]=max(extr(:,1));
Comments or corrections to the WAFO group