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

## CROSS-REFERENCE INFORMATION

This function calls:
 freqtype returns the frequency type of a Spectral density struct. spec2cov Computes covariance function and its derivatives specinterp Interpolation and zero-padding of spectrum diff Difference and approximate derivative. error Display message and abort function. getfield Get structure field contents. isstruct True for structures. sprintf Write formatted data to string. strcmpi Compare strings ignoring case. warning Display warning message; disable or enable warning messages.
This function is called by:
 spec2AcAt Evaluates survival function R(h1,h2)=P(Ac>h1,At>h2). spec2Acdf Evaluates cdf of crests P(Ac<=h) or troughs P(At<=h). spec2mmtpdf Calculates joint density of Maximum, minimum and period. spec2tccpdf Evaluates densities of wave period Tcc, wave lenght Lcc. spec2thpdf Joint density of amplitude and period/wave-length characteristics spec2tpdf Evaluates densities for crest-,trough-period, length. spec2tpdf2 Evaluates densities for various wave periods or wave lengths wminmax Calculates joint density of minimum and following maximum

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