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);
 
  See also  lcplot

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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

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