Home > wafo > wsim > spec2sdat.m

# spec2sdat

## PURPOSE

Simulates a Gaussian process and its derivative from spectrum

## SYNOPSIS

[x,xder]=spec2sdat(S,np,dt,iseed,method)

## DESCRIPTION

``` SPEC2SDAT Simulates a Gaussian process and its derivative from spectrum

CALL: [xs, xsder] = spec2sdat(S,[np cases],dt,iseed,method);

xs    = a cases+1 column matrix  ( t,X1(t) X2(t) ...).
xsder = a cases+1 column matrix  ( t,X1'(t) X2'(t) ...).
S     = a spectral density structure
np    = giving np load points.  (default length(S)-1=n-1).
If np>n-1 it is assummed that R(k)=0 for all k>n-1
cases = number of cases, i.e. number of replicates (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 = 'exact'  : simulation using cov2sdat
'random' : random phase and amplitude simulation (default)

SPEC2SDAT performs a fast and exact simulation of stationary zero mean
Gaussian process through circulant embedding of the covariance matrix
or by summation of sinus functions with random amplitudes and random
phase angle.

If the spectrum has a non-empty field .tr, then the transformation is
applied to the simulated data, the result is a simulation of a transformed
Gaussian process.

NB! The method 'exact' simulation may give high frequency ripple when
used with a small dt. In this case the method 'random' works better.

Example:
np =100; dt = .2;
[x1 x2] = spec2sdat(jonswap,np,dt);
waveplot(x1,'r',x2,'g',1,1)

## CROSS-REFERENCE INFORMATION

This function calls:
 cov2sdat Simulates a Gaussian process and its derivative freqtype returns the frequency type of a Spectral density struct. lagtype Returns the lag type of a Covariance struct. spec2cov Computes covariance function and its derivatives specinterp Interpolation and zero-padding of spectrum tranproc Transforms process X and up to four derivatives wnormspec Normalize a spectral density such that m0=m2=1 clear Clear variables and functions from memory. error Display message and abort function. fft Discrete Fourier transform. getfield Get structure field contents. interp1 1-D interpolation (table lookup) interp1q Quick 1-D linear interpolation. isfield True if field is in structure array. linspace Linearly spaced vector. strcmpi Compare strings ignoring case.
This function is called by:
 Chapter1 % CHAPTER1 demonstrates some applications of WAFO Chapter2 % CHAPTER2 Modelling random loads and stochastic waves Chapter4 % CHAPTER4 contains the commands used in Chapter 4 of the tutorial itmkurs_lab3 Script to computer exercises 3 test_cycles Quick test of the routines in module 'cycles' wafoinit setup all global variables of the WAFODEMO

## SOURCE CODE

```001 function [x,xder]=spec2sdat(S,np,dt,iseed,method)
002 %SPEC2SDAT Simulates a Gaussian process and its derivative from spectrum
003 %
004 %   CALL: [xs, xsder] = spec2sdat(S,[np cases],dt,iseed,method);
005 %
006 %         xs    = a cases+1 column matrix  ( t,X1(t) X2(t) ...).
007 %         xsder = a cases+1 column matrix  ( t,X1'(t) X2'(t) ...).
008 %         S     = a spectral density structure
009 %         np    = giving np load points.  (default length(S)-1=n-1).
010 %                 If np>n-1 it is assummed that R(k)=0 for all k>n-1
011 %         cases = number of cases, i.e. number of replicates (default=1)
012 %            dt = step in grid (default dt is defined by the Nyquist freq)
013 %         iseed = starting seed number for the random number generator
014 %                 (default none is set)
015 %        method = 'exact'  : simulation using cov2sdat
016 %                 'random' : random phase and amplitude simulation (default)
017 %
018 %  SPEC2SDAT performs a fast and exact simulation of stationary zero mean
019 %  Gaussian process through circulant embedding of the covariance matrix
020 %  or by summation of sinus functions with random amplitudes and random
021 %  phase angle.
022 %
023 %  If the spectrum has a non-empty field .tr, then the transformation is
024 %  applied to the simulated data, the result is a simulation of a transformed
025 %  Gaussian process.
026 %
027 %  NB! The method 'exact' simulation may give high frequency ripple when
028 %  used with a small dt. In this case the method 'random' works better.
029 %
030 % Example:
031 %  np =100; dt = .2;
032 % [x1 x2] = spec2sdat(jonswap,np,dt);
033 % waveplot(x1,'r',x2,'g',1,1)
034 %
036
037 % Reference
038 % C.R Dietrich and G. N. Newsam (1997)
039 % "Fast and exact simulation of stationary
040 %  Gaussian process through circulant embedding
041 %  of the Covariance matrix"
042 %  SIAM J. SCI. COMPT. Vol 18, No 4, pp. 1088-1107
043 %
044 % Hudspeth, R.T. and Borgman, L.E. (1979)
045 % "Efficient FFT simulation of Digital Time sequences"
046 % Journal of the Engineering Mechanics Division, ASCE, Vol. 105, No. EM2,
047 %
048
049 % Tested on: Matlab 5.3
050 % History:
051 % Revised pab Mar2005
052 %
053 % Revised pab Feb2004
054 % - changed seed to state
055 % revised pab 21.07.2002
056 % -
057 % revised pab 12.10.2001
059 % - the derivative, xder, is now calculated
060 %   correctly using fft for method = 'random'.
061 %   derivate.m is no longer needed!
062 % revised pab 11.10.2001
063 % - added a trick to avoid adding high frequency noise to spectrum when
064 %   method = 'exact'
065 % revised pab 21.01.2001
066 % - random is now available for 1D wavenumber spectra as well
067 % revised ir 11.06.00 - introducing dt
068 % revised es 24.05.00 - fixed default value for np, and small changes to help
069 % revised pab 13.03.2000
070 % - fixed default value for np
071 % revised pab 24.01.2000
072 % - added method random from L. Borgman fwavsim (slightly faster and needs
073 %   less memory)
075 % revised es 12.01.2000:  enable dir. spectrum input (change of spec2cov call)
076 % revised pab 12.10.1999
077 %  simplified call by calling cov2sdat
079
080 if nargin<5|isempty(method)
081   method='random';% 'exact' is slightly slower than method=='random'
082 end
083 if nargin<4|isempty(iseed)
084  iseed=[];
085 else
086   try
087     randn('state',iseed)
088   catch
089     randn('seed',iseed); % set the the seed
090   end
091 end
092
093 if nargin<3|isempty(dt)
094  S1=S;
095 else
096   S1=specinterp(S,dt); % interpolate spectrum
097 end
098
099 ftype = freqtype(S);
100 freq  = getfield(S,ftype);
101 Nt = length(freq);
102
103 S = S1;
104
105 if nargin<2|isempty(np)
106   np=Nt-1;
107   cases=1;
108 end
109 if strcmpi(method,'exact')
110   R = spec2cov(S,0,[],1,0,0);
111
112   if strcmpi(ftype,'f')
113     T = Nt/freq(end)/2;
114   else
115     T = Nt*pi/freq(end);
116   end
117   ltype = lagtype(R);
118   ix = min(find(getfield(R,ltype)>T));
119   % Trick to avoid adding high frequency noise to the spectrum
120   if ~isempty(ix),R.R(ix:end)=0; end
121
122   if nargout>1
123     [x, xder]=cov2sdat(R,np,iseed);
124   else
125     x=cov2sdat(R,np,iseed);
126   end
127   return
128 end
129
130 switch  length(np)
131   case 1, cases=1;
132   case 2, cases=np(2); np=np(1);
133   otherwise, error('Wrong input. Too many arguments')
134 end
135 if mod(np,2),np=np+1;end % make sure it is even
136
137 %ftype = freqtype(S); %options are 'f' and 'w' and 'k'
138
139
140
141 if 0
142   % Do the simulation with normalized spectrum
143   % Normalize the spectrum
144   [Sn,mn4,m0,m2] = wnormspec(S);
145   % Normalization constants
146   tnorm = 2*pi*sqrt(m0/m2);
147   xnorm = sqrt( tnorm*m0/(2*pi));
148   fs    = getfield(Sn,ftype);
149   Si    = Sn.S(2:end-1);
150 else
151   tnorm = 1;
152   xnorm = 1;
153
154   fs    = getfield(S,ftype);
155   Si    = S.S(2:end-1);
156
157   switch ftype
158     case 'f'
159     case {'w','k'}
160       Si=Si*2*pi;
161       fs=fs/2/pi;
162     otherwise
163       error('Not implemented for wavenumber spectra')
164   end
165 end
166
167 x=zeros(np,cases+1);
168 if nargout==2
169   xder=x;
170 end
171
172
173 dT = 1/(2*fs(end)); % dT
174
175 df = 1/(np*dT);
176
177
178 % interpolate for freq.  [1:(N/2)-1]*df and create 2-sided, uncentered spectra
179 % ----------------------------------------------------------------------------
180 f = [1:(np/2)-1]'*df;
181
182 fs(1)   = [];
183 fs(end) = [];
184 Fs      = [0; fs(:); (np/2)*df];
185 Su      = [0; abs(Si(:))/2; 0];
186
187 Smax = max(Su);
188 %Si = interp1(Fs,Su,f,'linear');
189 Si = interp1q(Fs,Su,f);
190
191 tmp  = find(Si>Smax*1e-3);
192 nmax = max(tmp)+1;
193 nmin = min(tmp)+1;
194 if isempty(nmax),nmax = np/2;end
195 if isempty(nmin),nmin = 2;end % Must always be greater than 1
196
197 Su=[0; Si; 0; Si((np/2)-1:-1:1)];
198
199 clear Si Fs
200
201 % Generate standard normal random numbers for the simulations
202 % -----------------------------------------------------------
203 Zr = randn((np/2)+1,cases);
204 Zi = [zeros(1,cases); randn((np/2)-1,cases); zeros(1,cases)];
205
206 A                = zeros(np,cases);
207 A(1:(np/2+1),:)  = Zr - sqrt(-1)*Zi; clear Zr Zi
208 A((np/2+2):np,:) = conj(A(np/2:-1:2,:));
209 A(1,:)           = A(1,:)*sqrt(2);
210 A((np/2)+1,:)    = A((np/2)+1,:)*sqrt(2);
211
212
213 % Make simulated time series
214 % --------------------------
215
216 T    = (np-1)*dT;
217 Ssqr = sqrt(Su*df/2);
218 %max(Ssqr)
219 % stochastic amplitude
220 A    = A.*Ssqr(:,ones(1,cases));
221 %max(abs(A))
222
223 % Deterministic amplitude
224 %A     =  sqrt(2)*Ssqr(:,ones(1,cases)).*exp(sqrt(-1)*atan2(imag(A),real(A)));
225 clear Su Ssqr
226
227
228 x(:,2:end) = xnorm*real(fft(A));
229
230 x(:,1)     = linspace(0,T*tnorm,np)'; %(0:dT:(np-1)*dT).';
231
232
233
234
235 if nargout==2,
236   %xder=derivate(x(:,1),x(:,2:end));
237
238   % new call pab 12.10.2001
239   w = 2*pi*[0; f;0;-f(end:-1:1)] ;
240   A = -i*A.*w(:,ones(1,cases));
241   xder(:,2:(cases+1)) = real(fft(A))*xnorm/tnorm;
242   xder(:,1)           = x(:,1);
243 end
244
245
246 if isfield(S,'tr') & ~isempty(S.tr)
247   disp('   Transforming data.')
248   g=S.tr;
249   G=fliplr(g); % the invers of g
250   if nargout==2, % gaus2dat
251     for ix=1:cases
252       tmp=tranproc([x(:,ix+1) xder(:,ix+1)],G);
253       x(:,ix+1)=tmp(:,1);
254       xder(:,ix+1)=tmp(:,2);
255     end
256   else
257     for ix=1:cases % gaus2dat
258       x(:,ix+1)=tranproc(x(:,ix+1),G);
259     end
260   end
261 end
262
263
264
265 return
266
267 % Old call kept just in case
268
269
270
271 df=1/(np*dT);
272
273
274 % Interpolate for freq.  [1:(N/2)-1]*df and create 2-sided, uncentered spectra
275 % ----------------------------------------------------------------------------
276 f=[1:(np/2)-1]'*df;
277
278 fs(1)=[];fs(end)=[];
279 Fs=[0; fs(:); (np/2)*df];
280 Su=[0; Si(:); 0];
281 Si = interp1(Fs,Su,f,'linear')/2;
282
283 Su=[0; Si; 0; Si((np/2)-1:-1:1)];
284 clear Si
285
286 % Generate standard normal random numbers for the simulations
287 % -----------------------------------------------------------
288 %Zr=randn((np/2)+1,cases);
289 %Zi=[zeros(1,cases); randn((np/2)-1,cases); zeros(1,cases)];
290 %AA=Zr-sqrt(-1)*Zi;
291
292 AA = randn((np/2)+1,cases)-sqrt(-1)*[zeros(1,cases); randn((np/2)-1,cases); zeros(1,cases)];
293 A=[AA;  conj(AA(np/2:-1:2,:))];
294 clear AA
295 A(1,:)=A(1,:)*sqrt(2);
296 A((np/2)+1,:)=A((np/2)+1,:)*sqrt(2);
297
298 % Make simulated time series
299 % --------------------------
300 T    = np*dT;
301 Ssqr = sqrt(Su*df/2);
302 A    = A.*Ssqr(:,ones(1,cases));
303 clear Su Ssqr
304
305 x(:,2:end)=real(fft(A));
306
307 x(:,1)=(0:dT:(np-1)*dT).';```

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