Home > wafo > spec > private > dspec2dcov.m

# dspec2dcov

## PURPOSE

Return covariance function given a directional spectrum

## SYNOPSIS

R=dspec2dcov(S,nr,Nt,rate,Nx,Ny,dx,dy)

## DESCRIPTION

``` DSPEC2DCOV Return covariance function given a directional spectrum
(used in spec2cov.)

CALL: R=dspec2dcov(S,nr,Nt,rate,Nx,Ny,dx,dy)

DCOV = covariance structure (see datastructures)
S = spectrum struct, type 'dir' or 'freq' (see datastructures)
nr = 0,1,2, number of derivatives
Nt = number of time lags
rate = 1,2,4,8...2^r, interpolation rate for R
Nx,Ny = number of space lags, (default 10 resp. 0)
dx,dy = space lag step  (default such that a wave length is covered)```

## CROSS-REFERENCE INFORMATION

This function calls:
 createcov Covariance class constructor gravity returns the constant acceleration of gravity simpson Numerical integration with the Simpson method spec2mom Calculates spectral moments from spectrum spec2spec Transforms between different types of spectra w2k Translates from frequency to wave number clear Clear variables and functions from memory. ifft Inverse discrete Fourier transform. isfield True if field is in structure array. meshgrid X and Y arrays for 3-D plots. nextpow2 Next higher power of 2. strcat Concatenate strings. strcmpi Compare strings ignoring case. warning Display warning message; disable or enable warning messages.
This function is called by:
 spec2cov Computes covariance function and its derivatives

## SOURCE CODE

```001 function R=dspec2dcov(S,nr,Nt,rate,Nx,Ny,dx,dy)
002 %DSPEC2DCOV Return covariance function given a directional spectrum
003 %            (used in spec2cov.)
004 %
005 %  CALL: R=dspec2dcov(S,nr,Nt,rate,Nx,Ny,dx,dy)
006 %
007 %   DCOV = covariance structure (see datastructures)
008 %      S = spectrum struct, type 'dir' or 'freq' (see datastructures)
009 %     nr = 0,1,2, number of derivatives
010 %     Nt = number of time lags
011 %   rate = 1,2,4,8...2^r, interpolation rate for R
012 %  Nx,Ny = number of space lags, (default 10 resp. 0)
013 %  dx,dy = space lag step  (default such that a wave length is covered)
014 %
015
016 % NB! requires simpson
017
018 % Tested on: Matlab 5.3
019 % History:
020 % revised ir 05.04.01, theta modified by a rotation angle phi.
021 % revised by es 24.05.00, fixing shape of output from w2k, since w2k changed
022 % revised by es 10.11.1999, nr>2
023 % revised by es 13.10.1999
024 % by pab 16.04.1999
025
026 % TODO % Needs more work
027 %initialize
028 if strcmpi(S.type(end-2:end),'k2d')
029   S=spec2spec(S,'dir');
030 elseif strcmpi(S.type(end-2:end),'k1d')
031   S=spec2spec(S,'freq');
032 end
033 if nargin<5|isempty(Nx)
034   Nx=10;
035 end
036 if nargin<6 |isempty(Ny)
037   Ny=0;
038 end
039
040 g=gravity;
041 if isfield(S,'g')
042   g=S.g;
043 end
044 m=spec2mom(S,2,'xyt',1);
045 if length(m)==2 % result when S of type 'freq'
046   m=spec2mom(S,4,'t',1);  %m=[m0 m002 m004]
047   m=[m(1) m(3)/g^2 m(3)/g^2]; %m=[m0 m200 m020];
048 end
049 if nargin<7|isempty(dx)
050   dx=4*pi*sqrt(m(1)/m(2))/Nx; % divide the average wave into Nx/2
051 end
052
053 if nargin<8|isempty(dy),
054   if Ny>0
055     dy=dx;
056     if isinf(dy) % in case of Nx==0
057       dy=4*pi*sqrt(m(1)/m(3))/Ny;
058       dx=dy;
059     end
060     if isinf(dy) % in case of both Nx and Ny ==0
061       dx=1;dy=1; % just to get some value
062     end
063   else
064     dy=1;
065     if isinf(dx)
066       dx=dy;
067     end
068   end
069 end
070
071 vari='t';
072 if Nx>0
073   vari=strcat(vari,'x');
074 end
075 if Ny>0
076   vari=strcat(vari,'y');
077 end
078 R=createcov(nr,vari);
079 R.tr=S.tr;
080 R.h=S.h;
081 R.norm=S.norm;
082 if isfield(S,'note')
083   R.note=S.note;
084 end
085
086 warning off
087 comment=1;
088 if comment,
089   disp('   Calculating covariances.')
090 end
091 comment=0;
092 if isfield(S,'f')
093   S.w=S.f*2*pi;
094   S.S=S.S/2/pi;
095 end
096 if ~isfield(S,'g')
097   S.g=gravity;
098 end
099 if ~isfield(S,'theta')
100   S.theta=0;
101   S.S=S.S(:)';
102 end
103
104 nf= length(S.w); %number of frequencies
105 np= length(S.theta); % number of angles
106
107 if Nx==Ny & dx==dy,
108   symmetry=1;%only able to exploit the symmetry S(x,y,w)=conj(S(y,x,w))
109          %if Nx=Ny,dx==dy
110 else
111   symmetry=0;
112 end
113 if symmetry,  % exploit the symmetry S(x,y,w)=conj(S(y,x,w))
114   xi=((-Nx):0)'*dx;
115 else
116   xi=((-Nx):Nx)'*dx;
117 end
118 yi=((-Ny):Ny)'*dy;
119
120 [Xi Yi]=meshgrid(xi,yi);
121 NXYi=length(Xi(:));
122
123 if ~isfield(S,'phi') | isempty(S.phi), S.phi=0; end
124
125 theta=S.theta(:)-S.phi;     %make sure it is a column
126 S.theta=S.theta-S.phi;
127 Yi=Yi(:);Xi=Xi(:);
128
129 if comment,
130   disp('   ... initiate')
131 end
132
133 DS2=(S.S)*S.w(end); % Normalizing
134 Sxy=zeros(nf,NXYi);
135 Sxyx=Sxy;Sxyxx=Sxy;Sxyxxx=Sxy;Sxyxxxx=Sxy;
136 if np>1
137   [kCtheta,kStheta]=w2k(S.w,S.theta,S.h,S.g); %size np x nf
138   for ix=1:NXYi,
139     Stheta=DS2.*exp(i*(Xi(ix).*kCtheta +Yi(ix).*kStheta));
140     Sxy(:,ix)=simpson(theta,Stheta).';
141     if nr>0
142       Sxyx(:,ix)=simpson(theta,i*kCtheta.*Stheta).';
143       if nr>1
144     Sxyxx(:,ix)=simpson(theta,-kCtheta.^2.*Stheta).';
145     if nr>2
146       Sxyxxx(:,ix)=simpson(theta,(-i*kCtheta.^3).*Stheta).';
147       if nr>3
148         Sxyxxxx(:,ix)=simpson(theta,(kCtheta.^4).*Stheta).';
149       end
150     end
151       end
152     end
153   end
154   clear kStheta
155 else %theta=0
156   [kCtheta]=w2k(S.w,S.theta,S.h,S.g);
157   kCtheta=kCtheta(:)'; %size 1 x nf
158   Stheta=(ones(size(Xi))*DS2).*exp(i*(Xi*kCtheta));
159   Sxy=Stheta.';
160   if nr>0
161     Sxyx=(i*(ones(NXYi,1)*kCtheta).*Stheta).';
162     if nr>1
163       Sxyxx=(-(ones(NXYi,1)*kCtheta.^2).*Stheta).';
164       if nr>2
165     Sxyxxx=(-i*(ones(NXYi,1)*kCtheta.^3).*Stheta).';
166     if nr>3
167       Sxyxxxx=((ones(NXYi,1)*kCtheta.^4).*Stheta).';
168     end
169       end
170     end
171   end
172 end
173
174 clear DS2 kCtheta  Xi Yi theta Stheta
175 if symmetry,   % exploit the symmetry in calculating Sxy
176   Sxy=[Sxy conj(Sxy(1:nf,end-2*Ny-1:-1:1))];
177   if nr>0
178     Sxyx=[Sxyx -conj(Sxyx(1:nf,end-2*Ny-1:-1:1))];
179     if nr>1
180       Sxyxx=[Sxyxx conj(Sxyxx(1:nf,end-2*Ny-1:-1:1))];
181       if nr>2
182     Sxyxxx=[Sxyxxx -conj(Sxyxxx(1:nf,end-2*Ny-1:-1:1))];
183     if nr>3
184       Sxyxxxx=[Sxyxxxx conj(Sxyxxxx(1:nf,end-2*Ny-1:-1:1))];
185     end
186       end
187     end
188   end
189   xi=((-Nx):Nx)'*dx;
190   NXYi=length(xi)*length(yi);
191 end
192 % size(Sxy)=[nf NXYi]
193
194 if comment,
195   disp('   ... FFT')
196 end
197 nfft=rate*2^nextpow2(2*nf-2);
198 rat=nfft/(2*nf-2);
199 Sxy=[Sxy; zeros(nfft-(2*nf)+2,NXYi); conj(Sxy(nf-1:-1:2,:))];
200 cov=real(ifft(Sxy,nfft,1)).'*rat; % size(cov) is NXYi  x nfft
201 % size(cov) will now become 2*Ny+1 x 2*Nx+1 x Nt+1
202 R.R=rshape(cov,Nt,Nx,Ny);         % for rshape: see end of this file
203
204 if Nx>0, R.x=xi; end, if Ny>0, R.y=yi; end, R.t=(0:Nt)'*pi/(S.w(nf))/rat;
205
206 if nr>0
207   w=S.w(:);
208   w=[w ; zeros(nfft-2*nf+2,1) ;-w(nf-1:-1:2) ]*ones(1,NXYi);
209   Sxy=i*w.*Sxy;
210   cov=real(ifft(Sxy,nfft,1)).'*rat;
211   R.Rt=rshape(cov,Nt,Nx,Ny);                               %r_t
212   if Nx>0
213     Sxyx=[Sxyx; zeros(nfft-(2*nf)+2,NXYi); conj(Sxyx(nf-1:-1:2,:))];
214     cov=real(ifft(Sxyx,nfft,1)).'*rat;
215     R.Rx=rshape(cov,Nt,Nx,Ny);                             %r_x
216   end
217   if nr>1
218     Sxy=i*w.*Sxy;
219     cov=real(ifft(Sxy,nfft,1)).'*rat;
220     R.Rtt=rshape(cov,Nt,Nx,Ny);                            %r_tt
221     if Nx>0
222       Sxyxx=[Sxyxx; zeros(nfft-(2*nf)+2,NXYi); conj(Sxyxx(nf-1:-1:2,:))];
223       cov=real(ifft(Sxyxx,nfft,1)).'*rat;
224       R.Rxx=rshape(cov,Nt,Nx,Ny);                          %r_xx
225       cov=real(ifft(i*w.*Sxyx,nfft,1)).'*rat;
226       R.Rxt=rshape(cov,Nt,Nx,Ny);                          %r_xt
227     end
228     if nr>2
229       Sxy=i*w.*Sxy;
230       cov=real(ifft(Sxy,nfft,1)).'*rat;
231       R.Rttt=rshape(cov,Nt,Nx,Ny);                         %r_ttt
232       if Nx>0
233     Sxyxxx=[Sxyxxx; zeros(nfft-(2*nf)+2,NXYi); conj(Sxyxxx(nf-1:-1:2,:))];
234     cov=real(ifft(Sxyxxx,nfft,1)).'*rat;
235     R.Rxxx=rshape(cov,Nt,Nx,Ny);                       %r_xxx
236     cov=real(ifft(i*w.*Sxyxx,nfft,1)).'*rat;
237     R.Rxxt=rshape(cov,Nt,Nx,Ny);                       %r_xxt
238     cov=real(ifft(-w.^2.*Sxyx,nfft,1)).'*rat;
239     R.Rxtt=rshape(cov,Nt,Nx,Ny);                       %r_xtt
240       end
241        if nr>3
242     Sxy=i*w.*Sxy;
243     cov=real(ifft(Sxy,nfft,1)).'*rat;
244     R.Rtttt=rshape(cov,Nt,Nx,Ny);                      %r_tttt
245     if Nx>0
246       cov=real(ifft([Sxyxxxx; zeros(nfft-(2*nf)+2,NXYi);...
247            conj(Sxyxxxx(nf-1:-1:2,:))],nfft,1)).'*rat;
248       R.Rxxxx=rshape(cov,Nt,Nx,Ny);                    %r_xxxx
249       cov=real(ifft(-i*w.^3.*Sxyx,nfft,1)).'*rat;
250       R.Rxttt=rshape(cov,Nt,Nx,Ny);                    %r_xttt
251       cov=real(ifft(-w.^2.*Sxyxx,nfft,1)).'*rat;
252       R.Rxxtt=rshape(cov,Nt,Nx,Ny);                    %r_xxtt
253       cov=real(ifft(i*w.*Sxyxxx,nfft,1)).'*rat;
254       R.Rxxxt=rshape(cov,Nt,Nx,Ny);                    %r_xxxt
255     end
256       end
257     end
258   end
259 end
260 return %end dspec2dcov
261
262
263 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
264 function cov=rshape(cov,Nt,Nx,Ny)
265 % RSHAPE Change size of matrix, special for covariance function routine
266 %
267 %  CALL:   Cout=rshape(Cin,Nt,Nx,Ny)
268 %
269 %     Cout = matrix, size [2*Nx+1 2*Ny+1 Nt+1], singleton dimensions removed
270 %     Cin  = matrix, size=[m n] where m>=(2*Nx+1)*(2*Ny+1), n>=Nt+1
271 %           rows below (2*Nx+1)*(2*Ny+1) and columns beyond Nt+1 are neglected
272 %     Nt,Nx,Ny = integers defining  size of output matrix
273
274 if Nx>0
275   if Ny>0
276     cov=reshape(cov(1:(2*Nx+1)*(2*Ny+1),1:Nt+1),2*Ny+1,2*Nx+1,Nt+1);
277   else
278     cov=reshape(cov(1:2*Nx+1,1:Nt+1),2*Nx+1,Nt+1);
279   end
280 elseif Ny>0
281   cov=reshape(cov(1:2*Ny+1,1:Nt+1),2*Ny+1,Nt+1);
282 else
283   cov=reshape(cov(1:Nt+1),Nt+1);
284 end
285 return  %end rshape
286```

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