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
returns the frequency type of a Spectral density struct. | |
Computes covariance function and its derivatives | |
Interpolation and zero-padding of spectrum | |
Difference and approximate derivative. | |
Display message and abort function. | |
Get structure field contents. | |
True for structures. | |
Write formatted data to string. | |
Compare strings ignoring case. | |
Display warning message; disable or enable warning messages. |
Evaluates survival function R(h1,h2)=P(Ac>h1,At>h2). | |
Evaluates cdf of crests P(Ac<=h) or troughs P(At<=h). | |
Calculates joint density of Maximum, minimum and period. | |
Evaluates densities of wave period Tcc, wave lenght Lcc. | |
Joint density of amplitude and period/wave-length characteristics | |
Evaluates densities for crest-,trough-period, length. | |
Evaluates densities for various wave periods or wave lengths | |
Calculates joint density of minimum and following maximum |
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
Comments or corrections to the WAFO group