Home > wafo > trgauss > spec2skew.m

spec2skew

PURPOSE ^

Estimates the moments of 2'nd order non-linear waves

SYNOPSIS ^

[skew, kurt, ma, sa, Hs ,Hd]=spec2skew(S,h,method)

DESCRIPTION ^

 SPEC2SKEW Estimates the moments of 2'nd order non-linear waves 
 
  CALL: [skew, kurt, mean, sigma] = spec2skew(S,h,method);
 
    skew, kurt, 
    mean, sigma = skewness, kurtosis, mean and standard deviation,
                  respectively, of 2'nd order waves to the leading
                  order. (skew=kurt=0 for a Gaussian process)
              S = spectral density structure
              h = water depth (default S.h) 
    method      = 'approximate' method due to Marthinsen & Winterstein (default)
                  'eigenvalue'  method due to Kac and Siegert
 
   Skewness = kurtosis-3 = 0 for a Gaussian process.
   The mean, sigma, skewness and kurtosis are determined as follows:
 
   method == 'approximate':  due to Marthinsen and Winterstein
     mean  = 2 * int Hd(w1,w1)*S(w1) dw1
     sigma = sqrt(int S(w1) dw1)
     skew  = 6 * int int [Hs(w1,w2)+Hd(w1,w2)]*S(w1)*S(w2) dw1*dw2/m0^(3/2)
     kurt  = (4*skew/3)^2
 
   where Hs = sum frequency effects  and Hd = difference frequency effects
 
  method == 'eigenvalue'
   
    mean  = sum(E);                    
    sigma = sqrt(sum(C^2)+2*sum(E^2));   
    skew  = sum((6*C^2+8*E^2).*E)/sigma^3;
    kurt  = 3+48*sum((C^2+E^2).*E^2)/sigma^4;
 
   where
    h1 = sqrt(S*dw/2);
    C  = (ctranspose(V)*[h1;h1]);
    and E and V is the eigenvalues and eigenvectors, respectively, of the 2'order 
    transfer matrix. S is the spectrum and dw is the frequency spacing of S.
 
  Example:  Simulate a Transformed Gaussian process:
   Hm0=7;Tp=11;
   S = jonswap([],[Hm0 Tp]); [sk, ku, me]=spec2skew(S);
   g=hermitetr([],[Hm0/4 sk ku me]);  g2=[g(:,1), g(:,2)*Hm0/4];
   ys = spec2sdat(S,15000);   % Simulated in the Gaussian world
   xs = gaus2dat(ys,g2);      % Transformed to the real world
 
  See also  hermitetr, ochitr, lc2tr, dat2tr

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

001 function [skew, kurt, ma, sa, Hs ,Hd]=spec2skew(S,h,method)
002 %SPEC2SKEW Estimates the moments of 2'nd order non-linear waves 
003 %
004 % CALL: [skew, kurt, mean, sigma] = spec2skew(S,h,method);
005 %
006 %   skew, kurt, 
007 %   mean, sigma = skewness, kurtosis, mean and standard deviation,
008 %                 respectively, of 2'nd order waves to the leading
009 %                 order. (skew=kurt=0 for a Gaussian process)
010 %             S = spectral density structure
011 %             h = water depth (default S.h) 
012 %   method      = 'approximate' method due to Marthinsen & Winterstein (default)
013 %                 'eigenvalue'  method due to Kac and Siegert
014 %
015 %  Skewness = kurtosis-3 = 0 for a Gaussian process.
016 %  The mean, sigma, skewness and kurtosis are determined as follows:
017 %
018 %  method == 'approximate':  due to Marthinsen and Winterstein
019 %    mean  = 2 * int Hd(w1,w1)*S(w1) dw1
020 %    sigma = sqrt(int S(w1) dw1)
021 %    skew  = 6 * int int [Hs(w1,w2)+Hd(w1,w2)]*S(w1)*S(w2) dw1*dw2/m0^(3/2)
022 %    kurt  = (4*skew/3)^2
023 %
024 %  where Hs = sum frequency effects  and Hd = difference frequency effects
025 %
026 % method == 'eigenvalue'
027 %  
028 %   mean  = sum(E);                    
029 %   sigma = sqrt(sum(C^2)+2*sum(E^2));   
030 %   skew  = sum((6*C^2+8*E^2).*E)/sigma^3;
031 %   kurt  = 3+48*sum((C^2+E^2).*E^2)/sigma^4;
032 %
033 %  where
034 %   h1 = sqrt(S*dw/2);
035 %   C  = (ctranspose(V)*[h1;h1]);
036 %   and E and V is the eigenvalues and eigenvectors, respectively, of the 2'order 
037 %   transfer matrix. S is the spectrum and dw is the frequency spacing of S.
038 %
039 % Example:  Simulate a Transformed Gaussian process:
040 %  Hm0=7;Tp=11;
041 %  S = jonswap([],[Hm0 Tp]); [sk, ku, me]=spec2skew(S);
042 %  g=hermitetr([],[Hm0/4 sk ku me]);  g2=[g(:,1), g(:,2)*Hm0/4];
043 %  ys = spec2sdat(S,15000);   % Simulated in the Gaussian world
044 %  xs = gaus2dat(ys,g2);      % Transformed to the real world
045 %
046 % See also  hermitetr, ochitr, lc2tr, dat2tr
047 
048 % References:
049 % Langley, RS (1987)
050 % 'A statistical analysis of nonlinear random waves'
051 % Ocean Engineering, Vol 14, No 5, pp 389-407
052 %
053 % Marthinsen, T. and Winterstein, S.R (1992)
054 % 'On the skewness of random surface waves'
055 % In proceedings of the 2nd ISOPE Conference, San Francisco, 14-19 june.
056 %
057 % Winterstein, S.R, Ude, T.C. and Kleiven, G. (1994)
058 % "Springing and slow drift responses:
059 %  predicted extremes and fatigue vs. simulation"
060 % In Proc. 7th International behaviour of Offshore structures, (BOSS)
061 % Vol. 3, pp.1-15
062 
063 % tested on: matlab 5.2 
064 % History
065 % revised pab 22.07.2002
066 %  -fixed a bug in the calculaton of the mean for the 2'nd order process.
067 %  - added method
068 %  -updated the help header accordingly
069 % revised pab 12.11.2001
070 % - added Langley's version of the quadratic transfer functions.
071 % revised pab 09.01.2001
072 % -simplified calculation of quadratic transfer functions when h=inf
073 % revised pab 09.01.2001
074 % - changed kurtosis so that kurtosis correspond to what wkurtosis
075 % measures, i.e., kurtosis-3=0 for a Gaussian process
076 % by pab 01.03.2000
077 
078 
079 error(nargchk(1,3,nargin))
080 
081 % default options
082 opts.disp = 0;
083 
084 if nargin<2|isempty(h),  h = S.h; end
085 if nargin<3|isempty(method), method = 'approximate'; end
086 
087 S = spec2spec(S,'freq');
088 S = ttspec(S,'w');
089 
090 
091 S1 = S.S(:);
092 m0 = trapz(S.w,S1);
093 Nw = length(S.w);
094 
095 [Hs, Hd,Hdii] = qtf(S.w(:),h);
096 
097 %return
098 %skew=6/sqrt(m0)^3*simpson(S.w,simpson(S.w,(Hs+Hd).*S1(:,ones(1,Nw))).*S1.');
099 
100 Hspd = trapz(S.w,trapz(S.w,(Hs+Hd).*S1(:,ones(1,Nw))).*S1.');
101 switch lower(method(1))
102   case 'a', %approx : Marthinsen, T. and Winterstein, S.R (1992) method
103     if nargout>2
104       ma = 2*trapz(S.w,Hdii.*S1);
105     end
106     sa   = sqrt(m0);
107     skew = 6/sa^3*Hspd;
108     kurt = (4*skew/3).^2+3;
109   
110   case 'q', % quasi method
111     dw = diff(S.w(1:2));
112     tmp1 =sqrt(S1(:,ones(1,Nw)).*(S1(:,ones(1,Nw)).'))*dw; 
113     Hd = Hd.*tmp1;
114     Hs = Hs.*tmp1;
115     k = 6;
116     stop = 0;
117     while (~stop)
118       E = eigs([Hd,Hs;Hs,Hd],[],k);
119       stop = (length(find(abs(E)<1e-4))>0 | k>1200);
120       k = 2*k;
121     end
122   
123   
124     m02=2*sum(E.^2); % variance of 2'nd order contribution 
125   
126     %Hstd = 16*trapz(S.w,(Hdii.*S1).^2);
127     %Hstd = trapz(S.w,trapz(S.w,((Hs+Hd)+ 2*Hs.*Hd).*S1(:,ones(1,Nw))).*S1.');
128     ma   = 2*trapz(S.w,Hdii.*S1);
129     %m02  = Hstd-ma^2% variance of second order part
130     sa   = sqrt(m0+m02);
131     skew = 6/sa^3*Hspd;
132     kurt = (4*skew/3).^2+3;
133   case 'e', % Kac and Siegert eigenvalue analysis
134     dw = diff(S.w(1:2));
135     tmp1 =sqrt(S1(:,ones(1,Nw)).*(S1(:,ones(1,Nw)).'))*dw; 
136     Hd = Hd.*tmp1;
137     Hs = Hs.*tmp1;
138     k = 6;
139     stop = 0;
140     E2 = 1;
141  
142     while (~stop)
143       [V,D] = eigs([Hd,Hs;Hs,Hd],[],k);
144       E = diag(D);
145       stop = (length(find(abs(E)<1e-4))>0 | k>=min(2*Nw,1200));
146       k = min(2*k,2*Nw);
147     end
148   
149     
150     h1 = sqrt(S1*dw/2);
151     C  = (ctranspose(V)*[h1;h1])
152     
153     E2 = E.^2;
154     C2 = C.^2;
155   
156     ma   = sum(E);                     % mean 
157     sa   = sqrt(sum(C2)+2*sum(E2));    % standard deviation
158     skew = sum((6*C2+8*E2).*E)/sa^3;   % skewness
159     kurt = 3+48*sum((C2+E2).*E2)/sa^4; % kurtosis
160     otherwise error('Method is not available')
161 end
162 
163 
164 
165 return
166 
167 function [Hs, Hd,Hdii]=qtf(w,h)
168 % QTF Quadratic Transfer Function
169 %
170 % CALL: [Hs, Hd, Hdii]=qtf(w,h)
171 %
172 %  Hs   = sum frequency effects
173 %  Hd   = difference frequency effects
174 %  Hdii = diagonal of Hd
175 %  w    = angular frequencies
176 %  h    = water depth
177 
178 
179 Nw = length(w);
180 g  = gravity;
181 kw = w2k(w,0,h,g);
182 [k1, k2] = meshgrid(kw);
183 
184 
185 
186 Hd = zeros(size(k1));
187 if h==inf,% go here for faster calculations
188   Hs   = 0.25*(abs(k1)+abs(k2));
189   Hd   = -0.25*abs(abs(k1)-abs(k2));
190   Hdii = zeros(size(w));
191   return
192 end
193 
194 [w1, w2]=meshgrid(w);
195 
196 warning off %off % suppress warnings on division by zero
197 
198 w12  = (w1.*w2);
199 w1p2 = (w1+w2);
200 w1m2 = (w1-w2);
201 k12  = (k1.*k2);
202 k1p2 = (k1+k2);
203 k1m2 = abs(k1-k2);
204 if 0, % Langley
205   p1 = (-2*w1p2.*(k12*g^2-w12.^2)+...
206       w1.*(w2.^4-g^2*k2.^2)+w2.*(w1.^4-g^2*k1.^2))./(4.*w12);
207   p2= w1p2.^2.*cosh((k1p2).*h)-g*(k1p2).*sinh((k1p2).*h);
208   
209   Hs = -p1./p2.*w1p2.*cosh((k1p2).*h)/g-...
210       (k12*g^2-w12.^2)./(4*g*w12)+(w1.^2+w2.^2)/(4*g);
211   
212   p3 = (-2*w1m2.*(k12*g^2+w12.^2)-...
213       w1.*(w2.^4-g^2*k2.^2)+w2.*(w1.^4-g^2*k1.^2))./(4.*w12);
214   p4= w1m2.^2.*cosh(k1m2.*h)-g*(k1m2).*sinh((k1m2).*h);
215   
216    
217   Hd = -p3./p4.*(w1m2).*cosh((k1m2).*h)/g-...
218       (k12*g^2+w12.^2)./(4*g*w12)+(w1.^2+w2.^2)/(4*g);  
219 
220 else  % Marthinsen & Winterstein
221   tmp1 = 0.5*g*k12./w12;
222   tmp2 = 0.25/g*(w1.^2+w2.^2+w12);
223   Hs   = (tmp1-tmp2+0.25*g*(w1.*k2.^2+w2.*k1.^2)./(w12.*(w1p2))).....
224       ./(1-g*(k1p2)./(w1p2).^2.*tanh((k1p2).*h))+tmp2-0.5*tmp1; % OK
225   
226   tmp2 = 0.25/g*(w1.^2+w2.^2-w12); %OK
227   Hd   = (tmp1-tmp2-0.25*g*(w1.*k2.^2-w2.*k1.^2)./(w12.*(w1m2))).....
228       ./(1-g*(k1m2)./(w1m2).^2.*tanh((k1m2).*h))+tmp2-0.5*tmp1; % OK
229 end  
230 
231 tmp1 = 0.5*g*kw./(w.*sqrt(g*h));
232 tmp2 = 0.25*w.^2/g;
233 
234 
235 Cg   = 0.5*g*(tanh(kw*h) +kw*h.*(1- tanh(kw*h).^2))./w; %Wave group velocity
236 Hdii = 0.5*(0.5*g*(kw./w).^2-0.5*w.^2/g+g*kw./(w.*Cg))....
237       ./(1-g*h./Cg.^2)-0.5*kw./sinh(2*kw*h); % OK
238 Hd(1:Nw+1:end) = Hdii;
239 
240 %k    = find(w1==w2);
241 %Hd(k) = Hdii;
242 
243 % The NaN's occur due to division by zero. => Set the isnans to zero
244 Hdii(isnan(Hdii))=0;
245 Hd(isnan(Hd))=0;
246 Hs(isnan(Hs))=0;
247 
248 warning on
249 return
250 
251 
252

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