% 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.
Plots the auto covariance function (ACF) 1D or 2D. | |
Estimate auto covariance function from data. | |
Transforms x using the transformation g. | |
Extracts level-crossing spectrum from data, | |
Estimate one-sided spectral density from data. | |
Estimate one-sided spectral density, version 2. | |
Extracts turning points from data, | |
Estimate transformation, g, from data. | |
Finds the indices to spurious points in a timeseries | |
Calculate transformation, g, proposed by Winterstein | |
Calculates (and plots) a JONSWAP spectral density | |
Test if a stochastic process is Gaussian. | |
Make a directional spectrum | |
factor for transforming spectra to finite water depth spectra | |
reconstruct the spurious/missing points of timeseries | |
Computes covariance function and its derivatives | |
Computes sampling interval from Nyquist frequency in spectrum | |
Calculates spectral moments from spectrum | |
Simulates a Gaussian process and its derivative from spectrum | |
Estimates the moments of 2'nd order non-linear waves | |
Transforms between different types of spectra | |
Directional spreading functions | |
Calculates a double peaked (swell + wind) spectrum | |
Plots transformation, g, eg. estimated with dat2tr. | |
Prints a caption "made by WAFO" in current figure. | |
Plots the surface elevation of timeseries. | |
Computes sample kurtosis | |
Plots data on a Normal distribution paper | |
Computes sample skewness | |
Plot a spectral density |
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
Comments or corrections to the WAFO group