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)   
    
  See also  cov2spec, datastructures

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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 %   
043 % See also  cov2spec, datastructures 
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