COV2SDAT Simulates a Gaussian process and its derivative from covariance func. CALL: [xs, xsder ] = cov2sdat(R,[np cases],iseed); xs = a cases+1 column matrix ( t,X1(t) X2(t) ...). xsder = a cases+1 column matrix ( t,X1'(t) X2'(t) ...). R = a covariance structure np = giving np load points. (default length(S/R)-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) iseed = starting seed number for the random number generator (default none is set) COV2SDAT performs a Fast and exact simulation of stationary zero mean Gaussian process through circulant embedding of the Covariance matrix. If the covariance structure 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. Example: np = 100; [x1, x2] = cov2sdat(spec2cov(jonswap),np); waveplot(x1,'r',x2,'g',1,1) See also spec2sdat, tranproc
Returns the lag type of a Covariance struct. | |
Transforms process X and up to four derivatives | |
Display message and abort function. | |
Discrete Fourier transform. | |
Get structure field contents. | |
True if field is in structure array. | |
Next higher power of 2. |
generates conditionally simulated values | |
Test if a stochastic process is Gaussian. | |
reconstruct the spurious/missing points of timeseries | |
Simulates a Gaussian process and its derivative from spectrum |
001 function [x, xder]=cov2sdat(R,np,iseed) 002 %COV2SDAT Simulates a Gaussian process and its derivative 003 % from covariance func. 004 % 005 % CALL: [xs, xsder ] = cov2sdat(R,[np cases],iseed); 006 % 007 % xs = a cases+1 column matrix ( t,X1(t) X2(t) ...). 008 % xsder = a cases+1 column matrix ( t,X1'(t) X2'(t) ...). 009 % R = a covariance structure 010 % np = giving np load points. (default length(S/R)-1=n-1). 011 % If np>n-1 it is assummed that R(k)=0 for all k>n-1 012 % cases = number of cases, i.e. number of replicates (default=1) 013 % iseed = starting seed number for the random number generator 014 % (default none is set) 015 % 016 % 017 % COV2SDAT performs a Fast and exact simulation of stationary zero mean 018 % Gaussian process through circulant embedding of the Covariance matrix. 019 % 020 % If the covariance structure has a non-empty field .tr, then the 021 % transformation is applied to the simulated data, the result is a 022 % simulation of a transformed Gaussian process. 023 % 024 % Example: 025 % np = 100; 026 % [x1, x2] = cov2sdat(spec2cov(jonswap),np); 027 % waveplot(x1,'r',x2,'g',1,1) 028 % 029 % See also spec2sdat, tranproc 030 031 % Reference 032 % C.R Dietrich and G. N. Newsam (1997) 033 % "Fast and exact simulation of stationary 034 % Gaussian process through circulant embedding 035 % of the Covariance matrix" 036 % SIAM J. SCI. COMPUT. Vol 18, No 4, pp. 1088-1107 037 038 039 % tested on: Matlab 5.3 040 % History: 041 % Revised pab Feb2004 042 % - changed seed to state 043 % revised pab 12.10.2001 044 % - added example 045 % - the derivative, xder, is now calculated correctly using fft. 046 % derivate.m is no longer needed! 047 % revised pab 11.10.2001 048 % - added call to lagtype.m 049 % - added a new solution on removing high frequency noise due to the fact 050 % that the fix of 16.01.2000 was not sufficient. 051 % revised by es 24.05.00 Some changes in help 052 % revised pab 13.03.2000 053 % -commented out warning messages 054 % revised pab 24.01.2000 055 % -added iseed 056 % revised pab 16.01.2000 057 % -added a fix in line 134: truncating the values of the spectral density to 058 % zero if S/max(S) > trunc. This is to ensure that high frequency 059 % noise is not added to the simulated timeseries. 060 % -fixed transformation for the derivatives 061 % revised pab 12.11.1999 062 % -fixed transformation 063 % revised pab 12.10.1999 064 % last modified by Per A. Brodtkorb 19.08.98 065 066 % add a nugget effect to ensure that round off errors 067 % do not result in negative spectral estimates 068 nugget=0;%10^-12; 069 070 ACF=R.R; 071 [n,m]=size(ACF); 072 073 if n<m 074 b=m;m=n;n=b; 075 ACF=ACF'; 076 end 077 078 if n<2, 079 error('The input vector must have more than 2 elements!') 080 end 081 082 istime=1; 083 switch m 084 case 1, % OK 085 otherwise, error('Only capable of 1D simulation') 086 end 087 088 [y I]=max(ACF); 089 switch I, 090 case 1,% ACF starts with zero lag 091 otherwise, error('ACF does not have a maximum at zero lag') 092 end 093 % new call pab 11.10.2001 094 ltype = lagtype(R); 095 %names=fieldnames(R); 096 %ind=find(strcmp(names,'x')+strcmp(names,'y')+strcmp(names,'t')); 097 %options are 'f' and 'w' and 'k' 098 %ltype=lower(names{ind}); 099 100 ti=getfield(R,ltype); 101 if isempty(ti) 102 dT=1; 103 else 104 dT=ti(2)-ti(1); % ti 105 end 106 107 if nargin<2|isempty(np) 108 np=n-1; 109 cases=1; 110 else 111 switch length(np) 112 case 1, cases=1; 113 case 2, cases=np(2); np=np(1); 114 otherwise, error('Wrong input. Too many arguments') 115 end 116 end 117 if nargin<3|isempty(iseed) 118 else 119 try 120 randn('state',iseed); 121 catch 122 randn('seed',iseed); % set the the seed 123 end 124 end 125 x=zeros(np,cases+1); 126 127 128 if nargout==2 129 xder=x; 130 end 131 132 % add a nugget effect to ensure that round off errors 133 % do not result in negative spectral estimates 134 ACF(1)=ACF(1)+nugget; 135 136 % Fast and exact simulation of simulation of stationary 137 % Gaussian process throug circulant embedding of the 138 % Covariance matrix 139 if (abs(ACF(n))>eps),% assuming ACF(n+1)==0 140 m2=2*n-1; 141 nfft=2^nextpow2(max(m2,2*np)); 142 ACF=[ACF(1:n);zeros(nfft-m2,1);ACF(n:-1:2)]; 143 if 0, %(n<np), 144 disp('Warning: I am now assuming that ACF(k)=0 ') 145 disp('for k>MAXLAG.') 146 end 147 else % ACF(n)==0 148 m2=2*n-2; 149 nfft=2^nextpow2(max(m2,2*np)); 150 ACF=[ACF(1:n);zeros(nfft-m2,1);ACF(n-1:-1:2)]; 151 end 152 m2=2*n-2; 153 154 S=real(fft(ACF,nfft));% periodogram 155 156 157 [maxS I]=max(S(:)); 158 k=find(S<0); 159 if any(k),% 160 % disp('Warning: Not able to construct a nonnegative circulant ') 161 % disp('vector from the ACF. Apply the parzen windowfunction ') 162 % disp('to the ACF in order to avoid this.') 163 % disp('The returned result is now only an approximation.') 164 165 % truncating negative values to zero to ensure that 166 % that this noise is not added to the simulated timeseries 167 168 S(k)=0; 169 % New call pab 11.10.2001 170 ix = min(k(find(k>2*I))); 171 if any(ix), 172 % truncating all oscillating values above 2 times the peak 173 % frequency to zero to ensure that 174 % that high frequency noise is not added to 175 % the simulated timeseries. 176 S(ix:end-ix) =0; 177 end 178 end 179 180 181 trunc=1e-5; 182 k=find(S(I:end-I)/maxS<trunc); 183 if any(k),% 184 S(k+I-1)=0; 185 % truncating small values to zero to ensure that 186 % that high frequency noise is not added to 187 % the simulated timeseries 188 end 189 190 cases2=ceil(cases/2); 191 % Generate standard normal random numbers for the simulations 192 % ----------------------------------------------------------- 193 if 1, 194 epsi=randn(nfft,cases2)+sqrt(-1)*randn(nfft,cases2); 195 Ssqr=sqrt(S/(nfft)); %sqrt(S(wn)*dw ) 196 ephat=epsi.*Ssqr(:,ones(1,cases2)); 197 y=fft(ephat,nfft); 198 x(:,2:cases+1)=[real(y(3:np+2,1:cases2)) imag(y(3:np+2,1:floor(cases/2)))]; 199 else 200 epsi=randn(nfft/2+1,cases2)+sqrt(-1)*randn(nfft/2+1,cases2); 201 epsi(1,:)=real(epsi(1,:)); 202 epsi(end,:)=real(epsi(end,:)); 203 epsi=[epsi; conj(epsi((end-1):-1:2,:))]; 204 end 205 206 207 x(:,1)=(0:dT:(dT*(np-1)))'; 208 209 if nargout==2, 210 211 % xder(:,1:(cases+1))=derivate((-2*dT:dT:(dT*(np+1)))', ... 212 % [real(y(1:np+4,1:cases2)) imag(y(1:np+4,1:floor(cases/2)))]); 213 214 % New call: pab 12.10.2001 215 Ssqr = Ssqr.*[(0:(nfft/2))';-((nfft/2-1):-1:1)']*2*pi/nfft/dT; 216 ephat = epsi.*Ssqr(:,ones(1,cases2)); 217 y = fft(ephat,nfft); 218 xder(:,2:(cases+1))=[imag(y(3:np+2,1:cases2)) -real(y(3:np+2,1:floor(cases/2)))]; 219 xder(:,1)=x(:,1); 220 221 end 222 223 if isfield(R,'tr') & ~isempty(R.tr) 224 disp(' Transforming data.') 225 g=R.tr; 226 G=fliplr(g); % the invers of g 227 if nargout==2, % gaus2dat 228 for ix=1:cases 229 tmp=tranproc([x(:,ix+1) xder(:,ix+1)],G); 230 x(:,ix+1)=tmp(:,1); 231 xder(:,ix+1)=tmp(:,2); 232 end 233 else 234 for ix=1:cases % gaus2dat 235 x(:,ix+1)=tranproc(x(:,ix+1),G); 236 end 237 end 238 end 239
Comments or corrections to the WAFO group