Home > wafo > wsim > cov2sdat.m

# cov2sdat

## PURPOSE

Simulates a Gaussian process and its derivative

## SYNOPSIS

[x, xder]=cov2sdat(R,np,iseed)

## DESCRIPTION

``` 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)

## CROSS-REFERENCE INFORMATION

This function calls:
 lagtype Returns the lag type of a Covariance struct. tranproc Transforms process X and up to four derivatives error Display message and abort function. fft Discrete Fourier transform. getfield Get structure field contents. isfield True if field is in structure array. nextpow2 Next higher power of 2.
This function is called by:
 cov2csdat generates conditionally simulated values mctrtest Test if a stochastic process is Gaussian. reconstruct reconstruct the spurious/missing points of timeseries spec2sdat Simulates a Gaussian process and its derivative from spectrum

## SOURCE CODE

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

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