Home > wafo > onedim > mm2lc.m

# mm2lc

## PURPOSE

Extracts level-crossing spectrum from min2Max cycles.

## SYNOPSIS

[lc,alpha] = mm2lc(mm,def,plotflag,sa)

## DESCRIPTION

``` 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);

## CROSS-REFERENCE INFORMATION

This function calls:
 lcplot Plots level-crossing spectrum (lc) clear Clear variables and functions from memory. diff Difference and approximate derivative. error Display message and abort function. unique Set unique.
This function is called by:
 dat2lc Extracts level-crossing spectrum from data, dat2tr Estimate transformation, g, from data.

## SOURCE CODE

```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 %
029
030
031 % Tested on: matlab 5.3
032 % History:
033 % revised pab Feb2004
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));```

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