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
Simulates a Gaussian process and its derivative | |
returns the frequency type of a Spectral density struct. | |
Returns the lag type of a Covariance struct. | |
Computes covariance function and its derivatives | |
Interpolation and zero-padding of spectrum | |
Transforms process X and up to four derivatives | |
Normalize a spectral density such that m0=m2=1 | |
Clear variables and functions from memory. | |
Display message and abort function. | |
Discrete Fourier transform. | |
Get structure field contents. | |
1-D interpolation (table lookup) | |
Quick 1-D linear interpolation. | |
True if field is in structure array. | |
Linearly spaced vector. | |
Compare strings ignoring case. |
% CHAPTER1 demonstrates some applications of WAFO | |
% CHAPTER2 Modelling random loads and stochastic waves | |
% CHAPTER4 contains the commands used in Chapter 4 of the tutorial | |
Script to computer exercises 3 | |
Quick test of the routines in module 'cycles' | |
setup all global variables of the WAFODEMO |
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).';
Comments or corrections to the WAFO group