SPEC2NLSDAT Simulates a Randomized 2nd order non-linear wave X(t) given the spectral density S. CALL: [xs2 xs1] = spec2nlsdat(S,[np cases],dt,iseed,method,fnLimit); xs2 = a cases+1 column matrix ( t,X1(t) X2(t) ...). (1'st + 2'nd order components) xs1 = a cases+1 column matrix ( t,X1(t) X2(t) ...). (1'st order component) S = a spectral density structure np = giving np load points. (default length(S)-1=n-1). If np>n-1 it is assummed that S(k)=0 for all k>n-1 cases = number of cases (default=1) dt = step in grid (default dt is defined by the Nyquist freq) iseed = starting seed number for the random number generator (default none is set) method = 'apStochastic' : Random amplitude and phase (default) 'aDeterministic' : Deterministic amplitude and random phase 'apDeterministic' : Deterministic amplitude and phase fnLimit = normalized upper frequency limit of spectrum for 2'nd order components. The frequency is normalized with sqrt(gravity*tanh(kbar*waterDepth)/Amax)/(2*pi) (default sqrt(2), i.e., Convergence criterion). Other possible values are: sqrt(1/2) : No bump in trough criterion sqrt(pi/7) : Wave steepness criterion SPEC2NLSDAT performs a Fast simulation of Randomized 2nd order non-linear waves by summation of sinus functions with random amplitudes and phase angles. The extent to which the simulated result are applicable to real seastates are dependent on the validity of the assumptions: 1) Seastate is unidirectional 2) the surface elevation is adequately represented by 2nd order random wave theory 3) The first order component of the surface elevation is a Gaussian random process. NOTE : If the spectrum does not decay rapidly enough towards zero, the contribution from the 2nd order wave components at the upper tail can be very large and unphysical. To ensure convergence of the perturbation series, the upper tail of the spectrum is truncated at FNLIMIT in the calculation of the 2nd order wave components, i.e., in the calculation of sum and difference frequency effects. This may also be combined with the elimination of second order effects from the spectrum, i.e., extract the linear components from the spectrum. One way to do this is to use SPEC2LINSPEC. Example: np =100; dt = .2; [x1, x2] = spec2nlsdat(jonswap,np,dt); waveplot(x1,'r',x2,'g',1,1) See also spec2linspec, spec2sdat, cov2sdat
Return difference- and sum-frequency effects. | |
returns the frequency type of a Spectral density struct. | |
returns the constant acceleration of gravity | |
Evaluates spectral characteristics and their covariance | |
Interpolation and zero-padding of spectrum | |
Translates from frequency to wave number | |
Clear variables and functions from memory. | |
Display message and abort function. | |
Discrete Fourier transform. | |
Get structure field contents. | |
Quick 1-D linear interpolation. | |
Linearly spaced vector. | |
Write formatted data to string. | |
Compare first N characters of strings ignoring case. |
Separates the linear component of the Spectrum |
001 function [x2,x,svec,dvec,A]=spec2nlsdat(S,np,dt,iseed,method,truncationLimit) 002 %SPEC2NLSDAT Simulates a Randomized 2nd order non-linear wave X(t) 003 % given the spectral density S. 004 % 005 % CALL: [xs2 xs1] = spec2nlsdat(S,[np cases],dt,iseed,method,fnLimit); 006 % 007 % xs2 = a cases+1 column matrix ( t,X1(t) X2(t) ...). 008 % (1'st + 2'nd order components) 009 % xs1 = a cases+1 column matrix ( t,X1(t) X2(t) ...). 010 % (1'st order component) 011 % S = a spectral density structure 012 % np = giving np load points. (default length(S)-1=n-1). 013 % If np>n-1 it is assummed that S(k)=0 for all k>n-1 014 % cases = number of cases (default=1) 015 % dt = step in grid (default dt is defined by the Nyquist freq) 016 % iseed = starting seed number for the random number generator 017 % (default none is set) 018 % method = 'apStochastic' : Random amplitude and phase (default) 019 % 'aDeterministic' : Deterministic amplitude and random phase 020 % 'apDeterministic' : Deterministic amplitude and phase 021 % fnLimit = normalized upper frequency limit of spectrum for 2'nd order 022 % components. The frequency is normalized with 023 % sqrt(gravity*tanh(kbar*waterDepth)/Amax)/(2*pi) 024 % (default sqrt(2), i.e., Convergence criterion). 025 % Other possible values are: 026 % sqrt(1/2) : No bump in trough criterion 027 % sqrt(pi/7) : Wave steepness criterion 028 % 029 % SPEC2NLSDAT performs a Fast simulation of Randomized 2nd order non-linear 030 % waves by summation of sinus functions with random amplitudes and 031 % phase angles. The extent to which the simulated result are applicable 032 % to real seastates are dependent on the validity of the assumptions: 033 % 034 % 1) Seastate is unidirectional 035 % 2) the surface elevation is adequately represented by 2nd order random 036 % wave theory 037 % 3) The first order component of the surface elevation is a Gaussian 038 % random process. 039 % 040 % NOTE : If the spectrum does not decay rapidly enough towards zero, the 041 % contribution from the 2nd order wave components at the upper tail can 042 % be very large and unphysical. 043 % To ensure convergence of the perturbation series, the upper tail of the 044 % spectrum is truncated at FNLIMIT in the calculation of the 2nd order 045 % wave components, i.e., in the calculation of sum and difference 046 % frequency effects. This may also be combined with the elimination of 047 % second order effects from the spectrum, i.e., extract the linear 048 % components from the spectrum. One way to do this is to use SPEC2LINSPEC. 049 % 050 % Example: 051 % np =100; dt = .2; 052 % [x1, x2] = spec2nlsdat(jonswap,np,dt); 053 % waveplot(x1,'r',x2,'g',1,1) 054 % 055 % See also spec2linspec, spec2sdat, cov2sdat 056 057 % Reference 058 % Nestegaard, A and Stokka T (1995) 059 % A Third Order Random Wave model. 060 % In proc.ISOPE conf., Vol III, pp 136-142. 061 % 062 % R. S Langley (1987) 063 % A statistical analysis of non-linear random waves. 064 % Ocean Engng, Vol 14, pp 389-407 065 % 066 % Marthinsen, T. and Winterstein, S.R (1992) 067 % 'On the skewness of random surface waves' 068 % In proc. ISOPE Conf., San Francisco, 14-19 june. 069 070 % tested on: Matlab 5.3 071 % History: 072 % Revised pab Feb2004 073 % - changed seed to state 074 % revised pab 25Jan2004 075 % - changed the truncation at the upper tail of the spectrum 076 % - added truncationLimit to input 077 % revised pab 11Nov2003 078 % changed call from disufq1 to disufq 079 % revised pab 22.07.2002 080 % revised pab 15.03.2002 081 % -new call to disufq 082 % - added nargchk 083 % by pab 21.01.2001 084 085 % TODO % Check the methods: 'apdeterministic' and 'adeterministic' 086 087 % Variables controlling the truncation of the spectrum for sum and 088 % difference frequency effects 089 reltol2ndorder = 1e-3; % 090 %truncationLimit = 1.5; 091 092 error(nargchk(1,6,nargin)) 093 094 ftype = freqtype(S); %options are 'f' and 'w' and 'k' 095 n = length(getfield(S,ftype)); 096 097 numWaves = 1000; % Typical number of waves in 3 hour seastate 098 %C = pi/7; % Wave steepness criterion 099 %C = 1/2 ; % No bump in trough 100 C = 2; % Convergence criterion as given in Nestegaard and Stokka (1995) 101 if (nargin<6 | isempty(truncationLimit)), 102 truncationLimit = sqrt(C); 103 end 104 if (nargin<5 | isempty(method)), method = 'apstochastic'; end 105 if (nargin>3 & ~isempty(iseed)), 106 try 107 randn('state',iseed); 108 catch 109 randn('seed',iseed); 110 end 111 end % set the the seed 112 if (nargin<2 | isempty(np)), np = n-1; cases = 1;end 113 if (nargin>2 & ~isempty(dt)), S = specinterp(S,dt);end % interpolate spectrum 114 115 switch length(np) 116 case 1, cases=1; 117 case 2, cases=np(2); np=np(1); 118 otherwise, error('Wrong input. Too many arguments') 119 end 120 np = np + mod(np,2); % make sure np is even 121 122 fs = getfield(S,ftype); 123 Si = S.S(2:end-1); 124 h = S.h; 125 if isempty(h), h = inf;end 126 127 128 switch ftype 129 case 'f' 130 case {'w','k'} 131 Si = Si*2*pi; 132 fs = fs/2/pi; 133 otherwise 134 error('Not implemented for wavenumber spectra') 135 end 136 dT = 1/(2*fs(end)); % dT 137 138 139 140 df = 1/(np*dT); 141 142 143 % interpolate for freq. [1:(N/2)-1]*df and create 2-sided, uncentered spectra 144 % ---------------------------------------------------------------------------- 145 f = [1:(np/2)-1]'*df; 146 147 fs(1) = []; 148 fs(end) = []; 149 Fs = [0; fs(:); (np/2)*df]; 150 Su = [0; abs(Si(:))/2; 0]; 151 152 Smax = max(Su); 153 %Si = interp1(Fs,Su,f,'linear'); 154 Si = interp1q(Fs,Su,f); 155 156 % If the spectrum does not decay rapidly enough towards zero, the 157 % contribution from the wave components at the upper tail can be very 158 % large and unphysical. 159 % To ensure convergence of the perturbation series, the upper tail of the 160 % spectrum is truncated in the calculation of sum and difference 161 % frequency effects. 162 % Find the critical wave frequency to ensure convergence. 163 tmp = find(Si>Smax*reltol2ndorder); 164 switch 2 165 case 1 166 nmax = max(tmp)+1; 167 nmin = min(tmp)+1; 168 case 2, 169 Hm0 = spec2char(S,'Hm0'); 170 Tm02 = spec2char(S,'Tm02'); 171 waterDepth = abs(S.h); 172 kbar = w2k(2*pi/Tm02,0,waterDepth); 173 174 Amax = sqrt(2*log(numWaves))*Hm0/4; % Expected maximum amplitude for 1000 waves seastate 175 176 fLimitUp = truncationLimit*sqrt(gravity*tanh(kbar*waterDepth)/Amax)/(2*pi); 177 fLimitLo = sqrt(gravity*tanh(kbar*waterDepth)*Amax/waterDepth^3)/(2*pi); 178 179 nmax = min(max(find(f<=fLimitUp)),max(find(Si>0)))+1; 180 nmin = max(min(find(fLimitLo<=f)),min(tmp))+1; 181 %nmin = min(tmp)+1; 182 end 183 if isempty(nmax),nmax = np/2;end 184 if isempty(nmin),nmin = 2;end % Must always be greater than 1 185 fLimitUp = df*nmax; 186 fLimitLo = df*nmin; 187 188 disp(sprintf('2nd order frequency Limits = %g,%g',fLimitLo, fLimitUp)) 189 190 Su = [0; Si; 0; Si((np/2)-1:-1:1)]; 191 192 clear Si Fs 193 194 T = (np-1)*dT; 195 x = zeros(np,cases+1); 196 x(:,1) = linspace(0,T,np)'; %(0:dT:(np-1)*dT).'; 197 x2 = x; 198 199 w = 2*pi*[0; f; np/2*df]; 200 g = gravity; 201 kw = w2k(w ,[],h,g); 202 203 % Generate standard normal random numbers for the simulations 204 % ----------------------------------------------------------- 205 Zr = randn((np/2)+1,cases); 206 Zi = [zeros(1,cases); randn((np/2)-1,cases); zeros(1,cases)]; 207 208 A = zeros(np,cases); 209 A(1:(np/2+1),:) = Zr - sqrt(-1)*Zi; clear Zr Zi 210 A((np/2+2):np,:) = conj(A(np/2:-1:2,:)); 211 A(1,:) = A(1,:)*sqrt(2); 212 A((np/2)+1,:) = A((np/2)+1,:)*sqrt(2); 213 214 % Make simulated time series 215 % -------------------------- 216 217 %mean(abs(A(:))) 218 Ssqr = sqrt(Su*df/2); 219 %max(Ssqr) 220 if strncmpi(method,'apdeterministic',3) 221 % Deterministic amplitude and phase 222 A(2:(np/2),:) = A(2,1); 223 A((np/2+2):np,:) = conj(A(2,1)); 224 A = sqrt(2)*Ssqr(:,ones(1,cases)).*exp(sqrt(-1)*atan2(imag(A),real(A)));; 225 elseif strncmpi(method,'adeterministic',3) 226 % Deterministic amplitude and random phase 227 A = sqrt(2)*Ssqr(:,ones(1,cases)).*... 228 exp(sqrt(-1)*atan2(imag(A),real(A))); 229 else 230 % stochastic amplitude and phase 231 A = A.*Ssqr(:,ones(1,cases)); 232 end 233 %max(abs(A)) 234 clear Su Ssqr 235 236 237 x(:,2:end) = real(fft(A)); 238 239 if nargout>3, 240 %compute the sum and frequency effects separately 241 [svec, dvec] = disufq((A.'),w,kw,min(h,10^30),g,nmin,nmax); 242 svec = svec.'; 243 dvec = dvec.'; 244 245 x2s = fft(svec); % 2'nd order sum frequency component 246 x2d = fft(dvec); % 2'nd order difference frequency component 247 248 % 1'st order + 2'nd order component. 249 x2(:,2:end) =x(:,2:end)+ real(x2s(1:np,:))+real(x2d(1:np,:)); 250 else 251 svec = disufq((A.'),w,kw,min(h,10^30),g,nmin,nmax).'; 252 253 x2o = fft(svec); % 2'nd order component 254 255 256 % 1'st order + 2'nd order component. 257 x2(:,2:end)=x(:,2:end)+ real(x2o(1:np,:)); 258 end 259 return 260 261 262 263 264
Comments or corrections to the WAFO group