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: This function is called by:

SUBFUNCTIONS ^

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