Home > wafo > wsim > spec2nlsdat.m

spec2nlsdat

PURPOSE

Simulates a Randomized 2nd order non-linear wave X(t)

SYNOPSIS

[x2,x,svec,dvec,A]=spec2nlsdat(S,np,dt,iseed,method,truncationLimit)

DESCRIPTION

``` 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)

CROSS-REFERENCE INFORMATION

This function calls:
 disufq Return difference- and sum-frequency effects. freqtype returns the frequency type of a Spectral density struct. gravity returns the constant acceleration of gravity spec2char Evaluates spectral characteristics and their covariance specinterp Interpolation and zero-padding of spectrum w2k Translates from frequency to wave number clear Clear variables and functions from memory. error Display message and abort function. fft Discrete Fourier transform. getfield Get structure field contents. interp1q Quick 1-D linear interpolation. linspace Linearly spaced vector. sprintf Write formatted data to string. strncmpi Compare first N characters of strings ignoring case.
This function is called by:
 spec2linspec Separates the linear component of the Spectrum

SOURCE CODE

```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 %
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
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)));;
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```

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