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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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 %  
035 % See also  cov2sdat, gaus2dat
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
058 % - added example
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)
074 % - added iseed
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
078 %    last modified by Per A. Brodtkorb 19.08.98
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