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)   
  
  See also  spec2linspec, spec2sdat, cov2sdat

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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

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