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';
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
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
%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
clf
S = dat2spec2(xx,9500);
wspecplot(S)
wafostamp([],'(ER)')
pause(pstate)
mom = spec2mom(S,4) [sa sqrt(mom(1))] pause(pstate)
mom = 0.2237 0.5266 8.1113 ans = 0.4730 0.4730
clf S1 = dat2spec2(xx,200); S2 = dat2spec2(xx,50); wspecplot(S1,[],'-.') hold on wspecplot(S2) hold off wafostamp([],'(ER)') pause(pstate)
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
rho3 = wskewness(xx(:,2)) rho4 = wkurtosis(xx(:,2)) [sk, ku]=spec2skew(S1)
rho3 = 0.2546 rho4 = 3.1739 sk = 0.1449 ku = 3.0373
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)
%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
indicates that the underlying distribution has a "heavy" upper tail and a "light" lower tail.
clf
wnormplot(xx(:,2))
wafostamp([],'(ER)')
pause(pstate)
clf Hm0 = 7; Tp = 11; spec = jonswap([],[Hm0 Tp]); spec.note pause(pstate)
ans = JONSWAP, Hm0 = 7, Tp = 11, gamma = 2.3853
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'
clf S1 =spec2spec(Sd,'freq'); S2 = spec2spec(Se,'enc'); wspecplot(spec), hold on wspecplot(S1,1,'.'), wspecplot(S2), wafostamp([],'(ER)') hold off pause(pstate)
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]
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)
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.
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
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.