Home > wafo > papers > tutorcom > Chapter2.m

Chapter2

PURPOSE ^

% CHAPTER2 Modelling random loads and stochastic waves

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

 % 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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 %% CHAPTER2 Modelling random loads and stochastic waves
002 %
003 % Chapter2 contains the commands used in Chapter 2 of the tutorial and
004 % present some tools for analysis of randodm functions with
005 % respect to their correlation, spectral and distributional properties.
006 % The code is divided into three examples: 
007 %
008 % Example1 is devoted to estimation of different parameters in the model.
009 % Example2 deals with spectral densities and
010 % Example3 presents the use of WAFO to simulate samples of a Gaussian
011 % process.
012 %
013 % Some of the commands are edited for fast computation. 
014 % Each set of commands is followed by a 'pause' command.
015 % 
016 
017 % Tested on Matlab 5.3, 7.01
018 % History
019 % Revised pab sept2005
020 %  Added sections -> easier to evaluate using cellmode evaluation.
021 % Revised pab Dec2004
022 % Created by GL July 13, 2000
023 % from commands used in Chapter 2
024 %
025 
026 pstate =  'off';
027 
028 %% Section 2.1 Introduction and preliminary analysis
029 %% Example 1: Sea data
030 %% Observed crossings compared to the expected for Gaussian signals
031 xx = load('sea.dat');
032 me = mean(xx(:,2))
033 sa = std(xx(:,2))
034 xx(:,2) = xx(:,2) - me;
035 lc = dat2lc(xx);
036 pause(pstate)
037 
038 
039 %% Turningpoints and irregularity factor
040 T = max(xx(:,1))-min(xx(:,1))
041 f0 = interp1(lc(:,1),lc(:,2),0)/T
042 pause(pstate)
043 
044 tp = dat2tp(xx);
045 alfa = f0/(length(tp)/(2*T))
046 pause(pstate)
047 
048 % A part of sea data is visulized with the following commands
049 clf
050 waveplot(xx,tp,'k-','*',1,1)
051 axis([0 2 -inf inf])
052 wafostamp([],'(ER)')
053 pause(pstate)
054 
055 %% Finding possible spurious points
056 %However, if the amount of data is too large for visual examinations one
057 %could use the following criteria to find possible spurious points. One
058 %must be careful using the criteria for extremevalue analysis, because
059 %these criteria might remove the highest and steepest waves.
060 dt = diff(xx(1:2,1));
061 dcrit = 5*dt;
062 ddcrit = 9.81/2*dt*dt;
063 zcrit = 0;
064 [inds, indg] = findoutliers(xx,zcrit,dcrit,ddcrit);
065 pause(pstate)
066 
067 %% Section 2.2 Frequency Modeling of Load Histories
068 %% Periodogram: Raw spectrum
069 clf
070 S = dat2spec2(xx,9500);
071 wspecplot(S)
072 wafostamp([],'(ER)')
073 pause(pstate)
074 
075 %% Calculate moments  
076 mom = spec2mom(S,4)
077 [sa sqrt(mom(1))]
078 pause(pstate)
079 
080 %% Section 2.2.1 Random functions in Spectral Domain - Gaussian processes
081 %% Smoothing of spectral estimate 
082 clf
083 S1 = dat2spec2(xx,200);
084 S2 = dat2spec2(xx,50);
085 wspecplot(S1,[],'-.')
086 hold on
087 wspecplot(S2)
088 hold off
089 wafostamp([],'(ER)')
090 pause(pstate)
091 
092 %% Estimated autocovariance
093 clf
094 R2 = spec2cov(S1,1);
095 Rest = dat2cov(xx,80,[],'- -');
096 covplot(R2,80,[],'.')
097 hold on
098 covplot(Rest)
099 wafostamp([],'(ER)')
100 hold off
101 pause(pstate)
102 
103 %% Section 2.2.2 Transformed Gaussian models
104 rho3 = wskewness(xx(:,2))
105 rho4 = wkurtosis(xx(:,2))
106 
107 [sk, ku]=spec2skew(S1)
108 
109 %% Comparisons of 3 transformations
110 clf
111 gh = hermitetr([],[sa sk ku me]);
112 g  = gh; g(:,2)=g(:,1)/sa;
113 trplot(g)
114 
115 [glc, test0] = dat2tr(xx);
116 hold on
117 plot(glc(:,1),glc(:,2),'b-')
118 plot(gh(:,1),gh(:,2),'b-.')
119 hold off
120 wafostamp([],'(ER)')
121 pause(pstate)
122 
123 %%  Test Gaussianity of a stochastic process.
124 %MCTRTEST simulates  e(g(u)-u) = int (g(u)-u)^2 du  for Gaussian processes 
125 %  given the spectral density, S. The result is plotted if test0 is given.
126 %  This is useful for testing if the process X(t) is Gaussian.
127 %  If 95% of TEST1 is less than TEST0 then X(t) is not Gaussian at a 5% level.
128       
129 %the following test takes time
130 N = length(xx);
131 test1 = mctrtest(S1,[N,50],test0);
132 wafostamp([],'(CR)')
133 pause(pstate)
134 
135 %% Normalplot of data xx
136 % indicates that the underlying distribution has a "heavy" upper tail and a
137 % "light" lower tail. 
138 clf
139 wnormplot(xx(:,2))
140 wafostamp([],'(ER)')
141 pause(pstate)
142 
143 %% Section 2.2.3 Spectral densities of sea data
144 %% Example 2: Different forms of spectra
145 clf
146 Hm0 = 7; Tp = 11;
147 spec = jonswap([],[Hm0 Tp]);
148 spec.note
149 pause(pstate)
150 
151 
152 %% Directional spectrum and Encountered directional spectrum
153 clf
154 D = spreading(101,'cos2s',0,[],spec.w,1)
155 Sd = mkdspec(spec,D)
156 pause(pstate)
157 
158 clf
159 Se = spec2spec(Sd,'encdir',0,10);
160 wspecplot(Se), hold on
161 wspecplot(Sd,1,'--'), hold off
162 wafostamp([],'(ER)')
163 pause(pstate)
164 
165 %% Frequency spectra
166 clf
167 S1 =spec2spec(Sd,'freq');
168 S2 = spec2spec(Se,'enc');
169 wspecplot(spec), hold on
170 wspecplot(S1,1,'.'),
171 wspecplot(S2),
172 wafostamp([],'(ER)')
173 hold off
174 pause(pstate)
175 
176 %% Wave number spectrum
177 clf
178 Sk = spec2spec(spec,'k1d')
179 Skd = spec2spec(Sd,'k1d')
180 wspecplot(Sk), hold on
181 wspecplot(Skd,1,'--'), hold off
182 wafostamp([],'(ER)')
183 pause(pstate)
184 
185 %% Effect of waterdepth on spectrum
186 clf
187 wspecplot(spec,1,'--'), hold on
188 S20 = spec;
189 S20.S = S20.S.*phi1(S20.w,20);
190 S20.h = 20;
191 wspecplot(S20),  hold off
192 wafostamp([],'(ER)')
193 pause(pstate)
194 
195 %% Section 2.3 Simulation of transformed Gaussian process
196 %% Example 3: Simulation of random sea    
197 % Reconstruct replaces the spurious points of seasurface by simulated
198 % data on the basis of the remaining data and a transformed Gaussian
199 % process. As noted previously one must be careful using the criteria 
200 % for finding spurious points when reconstructing a dataset, because
201 % these criteria might remove the highest and steepest waves as we can see
202 % in this plot where the spurious points is indicated with a '+' sign:
203 %
204 
205 clf
206 [y, grec] = reconstruct(xx,inds);
207 waveplot(y,'-',xx(inds,:),'+',1,1)
208 axis([0 inf -inf inf])
209 wafostamp([],'(ER)')
210 pause(pstate)
211 
212 %%
213 clf
214 L = 200;
215 x = dat2gaus(y,grec);
216 Sx = dat2spec(x,L);
217 pause(pstate)
218       
219 clf
220 dt = spec2dt(Sx)
221 Sx.tr = grec;
222 ysim = spec2sdat(Sx,480);
223 waveplot(ysim,'-')
224 wafostamp([],'(CR)')
225 pause(pstate)
226  
227 %% Estimated spectrum compared to Torsethaugen spectrum
228 clf
229 Tp = 1.1;
230 H0 = 4*sqrt(spec2mom(S1,1))
231 St = torsethaugen([0:0.01:5],[H0  2*pi/Tp]);
232 wspecplot(S1)
233 hold on
234 wspecplot(St,[],'-.')
235 wafostamp([],'(ER)')
236 pause(pstate)
237 
238 %%
239 clf
240 Snorm = St;
241 Snorm.S = Snorm.S/sa^2;
242 dt = spec2dt(Snorm)
243 pause(pstate)
244 
245 clf
246 [Sk Su] = spec2skew(St);
247 sa = sqrt(spec2mom(St,1));
248 gh = hermitetr([],[sa sk ku me]);
249 Snorm.tr = gh;
250 pause(pstate)
251 
252 %% Transformed Gaussian model compared to Gaussian model
253 clf
254 dt = 0.5;
255 ysim_t = spec2sdat(Snorm,240,dt);
256 xsim_t = dat2gaus(ysim_t,Snorm.tr);
257 pause(pstate)
258 
259 clf
260 xsim_t(:,2) = sa*xsim_t(:,2);
261 waveplot(xsim_t,ysim_t,5,1,sa,4.5,'r.','b')
262 wafostamp([],'(CR)')
263 pause(pstate)
264 
265

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