Home > wafo > spec > spec2cov.m

spec2cov

PURPOSE

Computes covariance function and its derivatives

SYNOPSIS

R = spec2cov(S,nr,Nt,rate,Nx,Ny,dx,dy)

DESCRIPTION

``` SPEC2COV Computes covariance function and its derivatives

CALL:  R = spec2cov(S,nr,Nt,rate,Nx,Ny,dx,dy);

R    = a covariance structure (See datastructures)
S    = a spectral density structure (See datastructures)
nr   = number of derivatives in output, nr<=4 (default = 0).
Nt   = number in time grid, i.e., number of time-lags
(default rate*(length(S.S)-1)).
rate = 1,2,4,8...2^r, interpolation rate for R
(default = 1, no interpolation)
Nx,Ny   = number in space grid (default = )
dx,dy   = space grid step (default: depending on S)

The input 'rate' gives together with the spectrum
the t-grid-spacing: dt=pi/(S.w(end)*rate), S.w(end) is the Nyquist freq.
This results in the t-grid: 0:dt:Nt*dt.

What output is achieved with different S and choices of Nt,Nx and Ny:
1) S.type='freq' or 'dir', Nt set, Nx,Ny not set: then result R(t) (one-dim)
2) S.type='k1d' or 'k2d', Nt set, Nx,Ny not set: then result R(x) (one-dim)
3) Any type, Nt and Nx set =>R(x,t); Nt and Ny set =>R(y,t)
4) Any type, Nt, Nx and Ny set => R(x,y,t)
5) Any type, Nt not set, Nx and/or Ny set => Nt set to default, goto 3) or 4)

NB! This routine requires that the spectrum grid is equidistant
starting from zero frequency.
NB! If you are using a model spectrum, S, with sharp edges
to calculate covariances then you should probably do like this:
Example:
S   = jonswap;
S.S([1:40 100:end]) = 0;
Nt  = length(S.S)-1;
R   = spec2cov(S,0,Nt);
win = parzen(2*Nt+1);
R.R = R.R.*win(Nt+1:end);
S1  = cov2spec(R);
R2  = spec2cov(S1);
plot(R2.R-R.R)
plot(S1.S-S.S)

CROSS-REFERENCE INFORMATION

This function calls:
 createcov Covariance class constructor dspec2dcov Return covariance function given a directional spectrum freqtype returns the frequency type of a Spectral density struct. spec2spec Transforms between different types of spectra diff Difference and approximate derivative. error Display message and abort function. fft Discrete Fourier transform. getfield Get structure field contents. isstruct True for structures. lower Convert string to lowercase. nextpow2 Next higher power of 2. setfield Set structure field contents. strcmpi Compare strings ignoring case.
This function is called by:
 Chapter2 % CHAPTER2 Modelling random loads and stochastic waves mctrtest Test if a stochastic process is Gaussian. spec2cov2 Computes covariance function and its derivatives, alternative version spec2sdat Simulates a Gaussian process and its derivative from spectrum wminmax Calculates joint density of minimum and following maximum

SOURCE CODE

```001 function  R = spec2cov(S,nr,Nt,rate,Nx,Ny,dx,dy)
002 %SPEC2COV Computes covariance function and its derivatives
003 %
004 %  CALL:  R = spec2cov(S,nr,Nt,rate,Nx,Ny,dx,dy);
005 %
006 %       R    = a covariance structure (See datastructures)
007 %       S    = a spectral density structure (See datastructures)
008 %       nr   = number of derivatives in output, nr<=4 (default = 0).
009 %       Nt   = number in time grid, i.e., number of time-lags
010 %              (default rate*(length(S.S)-1)).
011 %       rate = 1,2,4,8...2^r, interpolation rate for R
012 %               (default = 1, no interpolation)
013 %     Nx,Ny   = number in space grid (default = )
014 %     dx,dy   = space grid step (default: depending on S)
015 %
016 % The input 'rate' gives together with the spectrum
017 % the t-grid-spacing: dt=pi/(S.w(end)*rate), S.w(end) is the Nyquist freq.
018 % This results in the t-grid: 0:dt:Nt*dt.
019 %
020 % What output is achieved with different S and choices of Nt,Nx and Ny:
021 % 1) S.type='freq' or 'dir', Nt set, Nx,Ny not set: then result R(t) (one-dim)
022 % 2) S.type='k1d' or 'k2d', Nt set, Nx,Ny not set: then result R(x) (one-dim)
023 % 3) Any type, Nt and Nx set =>R(x,t); Nt and Ny set =>R(y,t)
024 % 4) Any type, Nt, Nx and Ny set => R(x,y,t)
025 % 5) Any type, Nt not set, Nx and/or Ny set => Nt set to default, goto 3) or 4)
026 %
027 % NB! This routine requires that the spectrum grid is equidistant
028 %     starting from zero frequency.
029 % NB! If you are using a model spectrum, S, with sharp edges
030 %     to calculate covariances then you should probably do like this:
031 % Example:
032 %    S   = jonswap;
033 %    S.S([1:40 100:end]) = 0;
034 %    Nt  = length(S.S)-1;
035 %    R   = spec2cov(S,0,Nt);
036 %    win = parzen(2*Nt+1);
037 %    R.R = R.R.*win(Nt+1:end);
038 %    S1  = cov2spec(R);
039 %    R2  = spec2cov(S1);
040 %    plot(R2.R-R.R)
041 %    plot(S1.S-S.S)
042 %
044
045 % NB! requires simpson
046
047 % tested on: matlab 5.3
048 % history:
049 % revised pab 21.11.2003
050 % - streamlined some code
051 % - updated help header
052 % - fixed bug in example
053 % revised by es 25.05.00, error if frequencies are not equidistant or  do not
054 %                         start from zero
055 % revised by es 23.05.00, call of freqtype, R.norm=S.norm
056 % revised by pab 18.11.1999
057 %   - fixed a bug when S=S(f)
058 % revised by es 13.10.1999
059 % revised by pab 23.08.1999
060
061 if ~isstruct(S)
062   error('Incorrect input spectrum, see help datastructures')
063 end
064
065 if nargin<2|isempty(nr)
066   nr=0; % number of derivatives
067 end
068
069 ftype = freqtype(S);
070 freq  = getfield(S,ftype);
071 n     = length(freq);
072
073 if all(freq>0) % for .type='k2d', negative frequencies are allowed
074   disp('Spectrum does not start at zero frequency/wave number.')
075   error('Correct it with specinterp, for example.')
076 %   disp('trying to fix it, but numerical problems may occur')
077 %   wfix=(0:dw:(ws-dw)).';
078 %   nfix=length(wfix);
079 %   specn=[zeros(nfix,1);specn];
080 %   w=[wfix;w];
081 %   n=n+nfix;
082 end
083 if any(abs(diff(diff(freq)))>1.0e-8)
084   disp('Not equidistant frequencies/wave numbers in spectrum.')
085   error('Correct it with specinterp, for example.')
086 end
087
088 if nargin<4|isempty(rate),
089   rate=1; %interpolation rate
090 else
091   if rate>16
092     rate=16;
093   end
094   rate=2^nextpow2(rate);%make sure rate is a power of 2
095 end
096 if nargin<3|isempty(Nt),
097   Nt = rate*(n-1);
098 else %check if Nt is ok
099   Nt = min(Nt,rate*(n-1));
100 end
101
102 if prod(size(S.S))~=length(S.S)&(nargin>4 & Nx==0)&(nargin>5 & Ny==0)
103   S = spec2spec(S,'freq');
104   % If Nx=Ny=0 then a twodimensional spectrum gives no extra information
105   % transform to type 'freq', and you do not have to call dspec2dcov
106 end
107 if prod(size(S.S))~=length(S.S)|(nargin>4 & Nx>0)|(nargin>5 & Ny>0)
108   if nargin < 5
109     Nx=[];Ny=[];dx=[];dy=[];
110   end
111   if nargin < 6
112     Ny=[];
113   end
114   if nargin < 7
115     dx=[];
116   end
117   if nargin < 8
118     dy=[];
119   end
120
121   R = dspec2dcov(S,nr,Nt,rate,Nx,Ny,dx,dy);
122   return % !!!!
123 end
124
125 if strcmpi(ftype,'k')
126   vari='x';
127 else
128   vari='t';
129 end
130
131 R      = createcov(nr,vari);
132 %R.R    = zeros(Nt+1,1);
133 R.tr   = S.tr;
134 R.h    = S.h;
135 R.norm = S.norm;
136 R.note = S.note;
137
138
139 %normalize spec so that sum(specn)/(n-1)=R(0)=var(X)
140 switch lower(ftype)
141  case {'w','k'},
142   w     = freq(:);
143   dT    = pi/w(n);
144   specn = S.S(:)*freq(n); %S.S(:)*pi/dT;
145  case 'f',
146   w     = 2*pi*freq(:);
147   dT    = 1/(2*freq(n));  % sampling interval=1/Fs
148   specn = S.S(:)*freq(n); %S.S(:)/(2*dT);
149  otherwise
150   error('unknown frequency type')
151 end
152
153 nfft = rate*2^nextpow2(2*n-2);
154
155 Rper = [specn; zeros(nfft-(2*n)+2,1) ; conj(specn(n-1:-1:2))]; % periodogram
156 t    = (0:Nt)'*dT*((2*n-2)/nfft);
157 %eval(['R.', vari,'=t;']);
158 R   = setfield(R,vari,t);
159 r   = real(fft(Rper,nfft))/(2*n-2);
160 %r   = real(fft(Rper/(2*n-2),nfft));
161 R.R = r(1:Nt+1);
162 if nr>0
163   w         = [w ; zeros(nfft-2*n+2,1) ;-w(n-1:-1:2) ];
164   fieldname = ['R' vari(ones(1,nr)) ];
165   for ix=1:nr,
166     Rper = (-i*w.*Rper);
167     r    = real(fft(Rper,nfft))/(2*n-2);
168     R    = setfield(R,fieldname(1:ix+1),r(1:Nt+1));
169     %eval(['R.R' vari(ones(1,ix)) '=r(1:(Nt+1));']);
170   end
171 end
172
173
174
175
176
177
178```

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