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
Covariance class constructor | |
Return covariance function given a directional spectrum | |
returns the frequency type of a Spectral density struct. | |
Transforms between different types of spectra | |
Difference and approximate derivative. | |
Display message and abort function. | |
Discrete Fourier transform. | |
Get structure field contents. | |
True for structures. | |
Convert string to lowercase. | |
Next higher power of 2. | |
Set structure field contents. | |
Compare strings ignoring case. |
% CHAPTER2 Modelling random loads and stochastic waves | |
Test if a stochastic process is Gaussian. | |
Computes covariance function and its derivatives, alternative version | |
Simulates a Gaussian process and its derivative from spectrum | |
Calculates joint density of minimum and following maximum |
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
Comments or corrections to the WAFO group