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)  
 
  See also  spec2sdat, tranproc

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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

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