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

## CROSS-REFERENCE INFORMATION

This function calls:
 gravity returns the constant acceleration of gravity spec2spec Transforms between different types of spectra ttspec Toggle Transform between angular frequency and frequency spectrum w2k Translates from frequency to wave number ctranspose ' Complex conjugate transpose. diff Difference and approximate derivative. eigs Find a few eigenvalues and eigenvectors of a matrix using ARPACK. error Display message and abort function. lower Convert string to lowercase. meshgrid X and Y arrays for 3-D plots. trapz Trapezoidal numerical integration. warning Display warning message; disable or enable warning messages.
This function is called by:
 Chapter2 % CHAPTER2 Modelling random loads and stochastic waves Chapter3 % CHAPTER3 Demonstrates distributions of wave characteristics Chapter4 % CHAPTER4 contains the commands used in Chapter 4 of the tutorial hermitetr Calculate transformation, g, proposed by Winterstein

## 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 %
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.
068 %  -updated the help header accordingly
069 % revised pab 12.11.2001
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