001 function R=dspec2dcov(S,nr,Nt,rate,Nx,Ny,dx,dy)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
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
046 m=spec2mom(S,4,'t',1);
047 m=[m(1) m(3)/g^2 m(3)/g^2];
048 end
049 if nargin<7|isempty(dx)
050 dx=4*pi*sqrt(m(1)/m(2))/Nx;
051 end
052
053 if nargin<8|isempty(dy),
054 if Ny>0
055 dy=dx;
056 if isinf(dy)
057 dy=4*pi*sqrt(m(1)/m(3))/Ny;
058 dx=dy;
059 end
060 if isinf(dy)
061 dx=1;dy=1;
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);
105 np= length(S.theta);
106
107 if Nx==Ny & dx==dy,
108 symmetry=1;
109
110 else
111 symmetry=0;
112 end
113 if symmetry,
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;
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);
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);
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
156 [kCtheta]=w2k(S.w,S.theta,S.h,S.g);
157 kCtheta=kCtheta(:)';
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,
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
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;
201
202 R.R=rshape(cov,Nt,Nx,Ny);
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);
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);
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);
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);
225 cov=real(ifft(i*w.*Sxyx,nfft,1)).'*rat;
226 R.Rxt=rshape(cov,Nt,Nx,Ny);
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);
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);
236 cov=real(ifft(i*w.*Sxyxx,nfft,1)).'*rat;
237 R.Rxxt=rshape(cov,Nt,Nx,Ny);
238 cov=real(ifft(-w.^2.*Sxyx,nfft,1)).'*rat;
239 R.Rxtt=rshape(cov,Nt,Nx,Ny);
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);
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);
249 cov=real(ifft(-i*w.^3.*Sxyx,nfft,1)).'*rat;
250 R.Rxttt=rshape(cov,Nt,Nx,Ny);
251 cov=real(ifft(-w.^2.*Sxyxx,nfft,1)).'*rat;
252 R.Rxxtt=rshape(cov,Nt,Nx,Ny);
253 cov=real(ifft(i*w.*Sxyxxx,nfft,1)).'*rat;
254 R.Rxxxt=rshape(cov,Nt,Nx,Ny);
255 end
256 end
257 end
258 end
259 end
260 return
261
262
263
264 function cov=rshape(cov,Nt,Nx,Ny)
265
266
267
268
269
270
271
272
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
286