Contents

CHAPTER2 Modelling random loads and stochastic waves

Chapter2 contains the commands used in Chapter 2 of the tutorial and present some tools for analysis of randodm functions with respect to their correlation, spectral and distributional properties. The code is divided into three examples:

Example1 is devoted to estimation of different parameters in the model. Example2 deals with spectral densities and Example3 presents the use of WAFO to simulate samples of a Gaussian process.

Some of the commands are edited for fast computation. Each set of commands is followed by a 'pause' command.

% Tested on Matlab 5.3, 7.01
% History
% Revised pab sept2005
%  Added sections -> easier to evaluate using cellmode evaluation.
% Revised pab Dec2004
% Created by GL July 13, 2000
% from commands used in Chapter 2
%

pstate =  'off';

Section 2.1 Introduction and preliminary analysis

Example 1: Sea data

Observed crossings compared to the expected for Gaussian signals

xx = load('sea.dat');
me = mean(xx(:,2))
sa = std(xx(:,2))
xx(:,2) = xx(:,2) - me;
lc = dat2lc(xx);
wafostamp([],'(ER)')
pause(pstate)
me =

  1.5441e-009


sa =

    0.4730

Turningpoints and irregularity factor

T = max(xx(:,1))-min(xx(:,1))
f0 = interp1(lc(:,1),lc(:,2),0)/T
pause(pstate)

tp = dat2tp(xx);
alfa = f0/(length(tp)/(2*T))
pause(pstate)

% A part of sea data is visulized with the following commands
clf
waveplot(xx,tp,'k-','*',1,1)
axis([0 2 -inf inf])
wafostamp([],'(ER)')
pause(pstate)
T =

  2.3808e+003


f0 =

    0.2263


alfa =

    0.4961

Finding possible spurious points

%However, if the amount of data is too large for visual examinations one
%could use the following criteria to find possible spurious points. One
%must be careful using the criteria for extremevalue analysis, because
%these criteria might remove the highest and steepest waves.
dt = diff(xx(1:2,1));
dcrit = 5*dt;
ddcrit = 9.81/2*dt*dt;
zcrit = 0;
[inds, indg] = findoutliers(xx,zcrit,dcrit,ddcrit);
pause(pstate)
Found 0 missing points
Found 0 spurious positive jumps of Dx
Found 0 spurious negative jumps of Dx
Found 37 spurious positive jumps of D^2x
Found 200 spurious negative jumps of D^2x
Found 244 consecutive equal values
Found the total of 1152 spurious points

Section 2.2 Frequency Modeling of Load Histories

Periodogram: Raw spectrum

clf
S = dat2spec2(xx,9500);
wspecplot(S)
wafostamp([],'(ER)')
pause(pstate)

Calculate moments

mom = spec2mom(S,4)
[sa sqrt(mom(1))]
pause(pstate)
mom =

    0.2237    0.5266    8.1113


ans =

    0.4730    0.4730

Section 2.2.1 Random functions in Spectral Domain - Gaussian processes

Smoothing of spectral estimate

clf
S1 = dat2spec2(xx,200);
S2 = dat2spec2(xx,50);
wspecplot(S1,[],'-.')
hold on
wspecplot(S2)
hold off
wafostamp([],'(ER)')
pause(pstate)

Estimated autocovariance

clf
R2 = spec2cov(S1,1);
Rest = dat2cov(xx,80,[],'- -');
covplot(R2,80,[],'.')
hold on
covplot(Rest)
wafostamp([],'(ER)')
hold off
pause(pstate)
ans =

    81     1

Section 2.2.2 Transformed Gaussian models

rho3 = wskewness(xx(:,2))
rho4 = wkurtosis(xx(:,2))

[sk, ku]=spec2skew(S1)
rho3 =

    0.2546


rho4 =

    3.1739


sk =

    0.1449


ku =

    3.0373

Comparisons of 3 transformations

clf
gh = hermitetr([],[sa sk ku me]);
g  = gh; g(:,2)=g(:,1)/sa;
trplot(g)

[glc, test0] = dat2tr(xx);
hold on
plot(glc(:,1),glc(:,2),'b-')
plot(gh(:,1),gh(:,2),'b-.')
hold off
wafostamp([],'(ER)')
pause(pstate)

Test Gaussianity of a stochastic process.

%MCTRTEST simulates  e(g(u)-u) = int (g(u)-u)^2 du  for Gaussian processes
%  given the spectral density, S. The result is plotted if test0 is given.
%  This is useful for testing if the process X(t) is Gaussian.
%  If 95% of TEST1 is less than TEST0 then X(t) is not Gaussian at a 5% level.

%the following test takes time
N = length(xx);
test1 = mctrtest(S1,[N,50],test0);
wafostamp([],'(CR)')
pause(pstate)
finished 1 of 3
finished 2 of 3
finished 3 of 3

Normalplot of data xx

indicates that the underlying distribution has a "heavy" upper tail and a "light" lower tail.

clf
wnormplot(xx(:,2))
wafostamp([],'(ER)')
pause(pstate)

Section 2.2.3 Spectral densities of sea data

Example 2: Different forms of spectra

clf
Hm0 = 7; Tp = 11;
spec = jonswap([],[Hm0 Tp]);
spec.note
pause(pstate)
ans =

JONSWAP, Hm0 = 7, Tp = 11, gamma = 2.3853

Directional spectrum and Encountered directional spectrum

clf
D = spreading(101,'cos2s',0,[],spec.w,1)
Sd = mkdspec(spec,D)
pause(pstate)

clf
Se = spec2spec(Sd,'encdir',0,10);
wspecplot(Se), hold on
wspecplot(Sd,1,'--'), hold off
wafostamp([],'(ER)')
pause(pstate)
D = 

        S: [101x257 double]
        w: [257x1 double]
    theta: [101x1 double]
     type: 'dir'
      phi: 0
     note: 'Spreading: cos2s'
      th0: 0
     data: [15 15 0.5200 5 -2.5000 0 1 Inf]


Sd = 

        S: [101x257 double]
        w: [257x1 double]
    theta: [101x1 double]
       tr: []
        h: Inf
     type: 'dir'
      phi: 0
     norm: 0
     note: 'JONSWAP, Hm0 = 7, Tp = 11, gamma = 2.3853; Spreading: cos2s'
     date: '03-Sep-2005 04:11:54'

Frequency spectra

clf
S1 =spec2spec(Sd,'freq');
S2 = spec2spec(Se,'enc');
wspecplot(spec), hold on
wspecplot(S1,1,'.'),
wspecplot(S2),
wafostamp([],'(ER)')
hold off
pause(pstate)

Wave number spectrum

clf
Sk = spec2spec(spec,'k1d')
Skd = spec2spec(Sd,'k1d')
wspecplot(Sk), hold on
wspecplot(Skd,1,'--'), hold off
wafostamp([],'(ER)')
pause(pstate)
Sk = 

       S: [1x514 double]
      tr: []
       h: Inf
    type: 'k1D'
     phi: 0
    norm: 0
    note: 'JONSWAP, Hm0 = 7, Tp = 11, gamma = 2.3853'
    date: '03-Sep-2005 04:11:57'
       g: 9.8063
       k: [1x514 double]


Skd = 

       S: [1x257 double]
      tr: []
       h: Inf
    type: 'k1d'
     phi: 0
    norm: 0
    note: 'JONSWAP, Hm0 = 7, Tp = 11, gamma = 2.3853; Spreading: cos2s'
    date: '03-Sep-2005 04:11:58'
       g: 9.8063
       k: [1x257 double]

Effect of waterdepth on spectrum

clf
wspecplot(spec,1,'--'), hold on
S20 = spec;
S20.S = S20.S.*phi1(S20.w,20);
S20.h = 20;
wspecplot(S20),  hold off
wafostamp([],'(ER)')
pause(pstate)

Section 2.3 Simulation of transformed Gaussian process

Example 3: Simulation of random sea

Reconstruct replaces the spurious points of seasurface by simulated data on the basis of the remaining data and a transformed Gaussian process. As noted previously one must be careful using the criteria for finding spurious points when reconstructing a dataset, because these criteria might remove the highest and steepest waves as we can see in this plot where the spurious points is indicated with a '+' sign:

clf
[y, grec] = reconstruct(xx,inds);
waveplot(y,'-',xx(inds,:),'+',1,1)
axis([0 inf -inf inf])
wafostamp([],'(ER)')
pause(pstate)
First reconstruction attempt,    e(g-u)=0.80807
Simulation nr: 1 of 20   e(g-g_old)=0.1611,  e(g-u)=0.80043
Simulation nr: 2 of 20   e(g-g_old)=0.16878,  e(g-u)=0.92456
Simulation nr: 3 of 20   e(g-g_old)=0.06987,  e(g-u)=0.96482
Simulation nr: 4 of 20   e(g-g_old)=0.036037,  e(g-u)=0.98469
Simulation nr: 5 of 20   e(g-g_old)=0.024314,  e(g-u)=0.9973
Simulation nr: 6 of 20   e(g-g_old)=0.020549,  e(g-u)=1.0087
Simulation nr: 7 of 20   e(g-g_old)=0.0068146,  e(g-u)=1.0124
Simulation nr: 8 of 20   e(g-g_old)=0.0012065,  e(g-u)=1.013
Simulation nr: 9 of 20   e(g-g_old)=0.00019174,  e(g-u)=1.013
Elapsed time is 37.714000 seconds.
clf
L = 200;
x = dat2gaus(y,grec);
Sx = dat2spec(x,L);
pause(pstate)

clf
dt = spec2dt(Sx)
Sx.tr = grec;
ysim = spec2sdat(Sx,480);
waveplot(ysim,'-')
wafostamp([],'(CR)')
pause(pstate)
dt =

    0.2500

   Transforming data.

Estimated spectrum compared to Torsethaugen spectrum

clf
Tp = 1.1;
H0 = 4*sqrt(spec2mom(S1,1))
St = torsethaugen([0:0.01:5],[H0  2*pi/Tp]);
wspecplot(S1)
hold on
wspecplot(St,[],'-.')
wafostamp([],'(ER)')
pause(pstate)
H0 =

    7.0000

Warning: Hm0 is outside the valid range
The validity of the spectral density is questionable
     Spectrum for Wind dominated sea
clf
Snorm = St;
Snorm.S = Snorm.S/sa^2;
dt = spec2dt(Snorm)
pause(pstate)

clf
[Sk Su] = spec2skew(St);
sa = sqrt(spec2mom(St,1));
gh = hermitetr([],[sa sk ku me]);
Snorm.tr = gh;
pause(pstate)
dt =

    0.6283

Transformed Gaussian model compared to Gaussian model

clf
dt = 0.5;
ysim_t = spec2sdat(Snorm,240,dt);
xsim_t = dat2gaus(ysim_t,Snorm.tr);
pause(pstate)

clf
xsim_t(:,2) = sa*xsim_t(:,2);
waveplot(xsim_t,ysim_t,5,1,sa,4.5,'r.','b')
wafostamp([],'(CR)')
pause(pstate)
   Transforming data.