Home > wafo > spec > spec2cov2.m

spec2cov2

PURPOSE ^

Computes covariance function and its derivatives, alternative version

SYNOPSIS ^

R = spec2cov2(S,nr,Nt,dt)

DESCRIPTION ^

 SPEC2COV2 Computes covariance function and its derivatives, alternative version 
   
     CALL:  R = spec2cov2(S,nr,Nt,dt); 
   
   R    = [R0, R1,...Rnr] matrix with autocovariance and its 
           derivatives, i.e., Ri (i=1:nr) are column vectors with 
           the 1'st to nr'th derivatives of R0.  size Nt+1 x Nr+1 
   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*(n-1)) where rate = round(1/(2*f(end)*dt)) or  
           rate = round(pi/(w(n)*dt)) depending on S. 
   dt   = time spacing for R 
    
    NB! This routine requires that the spectrum grid is equidistant 
        starting from zero frequency. 
    Example:     
       S = jonswap; 
       dt = 0.1;   
       R = spec2cov2(S,3,256,dt); 
  
  See also  spec2cov, specinterp, datastructures

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function R = spec2cov2(S,nr,Nt,dt) 
002 %SPEC2COV2 Computes covariance function and its derivatives, alternative version 
003 %  
004 %    CALL:  R = spec2cov2(S,nr,Nt,dt); 
005 %  
006 %  R    = [R0, R1,...Rnr] matrix with autocovariance and its 
007 %          derivatives, i.e., Ri (i=1:nr) are column vectors with 
008 %          the 1'st to nr'th derivatives of R0.  size Nt+1 x Nr+1 
009 %  S    = a spectral density structure (See datastructures) 
010 %  Nr   = number of derivatives in output, Nr<=4          (default 0) 
011 %  Nt   = number in time grid, i.e., number of time-lags.  
012 %         (default rate*(n-1)) where rate = round(1/(2*f(end)*dt)) or  
013 %          rate = round(pi/(w(n)*dt)) depending on S. 
014 %  dt   = time spacing for R 
015 %   
016 %   NB! This routine requires that the spectrum grid is equidistant 
017 %       starting from zero frequency. 
018 %   Example:     
019 %      S = jonswap; 
020 %      dt = 0.1;   
021 %      R = spec2cov2(S,3,256,dt); 
022 % 
023 % See also  spec2cov, specinterp, datastructures   
024    
025 % History 
026 % revised pab Sept 2005 
027 %  - rate is sometimes zero -> spec2cov2 crash: made an ad-hoc solution. 
028 % by pab 24.11.2003 
029 % based on code from spec2XXpdf programmes 
030    
031 error(nargchk(1,4,nargin)) 
032 if ~isstruct(S) 
033   error('Incorrect input spectrum, see help datastructures') 
034 end 
035    
036 if nargin<2|isempty(nr) 
037   nr=0; % number of derivatives 
038 end 
039 ftype = freqtype(S); %options are 'f' and 'w' and 'k' 
040 freq  = getfield(S,ftype); 
041 n     = length(freq); 
042  
043 if nargin<4|isempty(dt), 
044   if strcmpi(ftype,'f') 
045     dt=1/(2*freq(n)); % sampling interval=1/Fs 
046   else 
047     dt=pi/(freq(n)); 
048   end 
049 end 
050 if strcmpi(ftype,'f') 
051   rate = round(1/(2*freq(n)*dt)); % sampling interval=1/Fs 
052 else 
053   rate = round(pi/(freq(n)*dt)); 
054 end 
055 % TODO % rate sometimes is zero -> spec2cov2 sometimes crash 
056 rate = max(rate,1); 
057  
058 if nargin<3|isempty(Nt), 
059   Nt = rate*(n-1); 
060 else %check if Nt is ok 
061   Nt = min(Nt,rate*(n-1)); 
062 end 
063  
064 checkdt = 1.2*min(diff(freq))/2/pi; 
065 switch ftype 
066  case {'w'}, 
067   vari = 't';  
068  case {'f'},  
069   vari    = 't';  
070   checkdt = checkdt*2*pi; 
071  case {'k'}, 
072   vari = 'x'; 
073 end 
074 msg1 = sprintf(['The step dt = %g in computation of the density is too' ... 
075         ' small.'],dt); 
076 msg2 = sprintf(['The step dt = %g step is small, may cause numerical' ... 
077         ' inaccuracies.'],dt); 
078  
079 if (checkdt < 2^-16/dt) 
080   disp(msg1) 
081   disp('The computed covariance (by FFT(2^K)) may differ from the theoretical.') 
082   disp('Solution:') 
083   error('use larger dt or sparser grid for spectrum.') 
084 end 
085  
086 % Calculating covariances 
087 %~~~~~~~~~~~~~~~~~~~~~~~~ 
088 S2 = specinterp(S,dt); 
089 R2 = spec2cov(S2,nr,Nt,1,0,0); 
090 R  = zeros(Nt+1,nr+1); 
091 fieldname = ['R'  vari(ones(1,nr))]; 
092 idx = {1:Nt+1}; 
093 for ix = 1:nr+1 
094   R(:,ix) = getfield(R2,fieldname(1:ix),idx); 
095 end 
096  
097 EPS0 = 0.0001; 
098 cc  = R(1,1)-R(2,1)*(R(2,1)/R(1,1)); 
099 if Nt+1>=5 
100   %cc1=R(1,1)-R(3,1)*(R(3,1)/R(1,1))+R(3,2)*(R(3,2)/R(1,3)); 
101   %cc3=R(1,1)-R(5,1)*(R(5,1)/R(1,1))+R(5,2)*(R(5,2)/R(1,3)); 
102    
103   cc2 = R(1,1)-R(5,1)*(R(5,1)/R(1,1)); 
104   if (cc2<EPS0) 
105     warning(msg1) 
106   end 
107 end 
108 if (cc<EPS0) 
109   disp(msg2) 
110 end 
111 return

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