% Define
PI=3.14159165;
dt=0.004;
dx=30;
dz=14;
% Get Parameters:
%---------------
nt=1024;
nz=200;
nx=257;
%---------------
ntfft =1040;
nxfft =260;
fw=0.0;
fz=0;
%---------------
nw = ntfft/2+1;
dw = 2.0*PI/(ntfft*dt);
fw = 0; % w ranges from 0 to Pi
nk = nxfft;
dk = 2.0*PI/(nxfft*dx);
fk = -PI/dx; % k rangs from -Pi to Pi, coz the Fourier Kernel is 1��
% clarification matrix
cresult=zeros(nx,nz);
CQ1=zeros(nk,nw);
% Read Data
% Load traces into the zero-offset array
% [D,H] = readsegy('spike.su');
load data.mat
% Load velicoty file
% old mathod
load v.mat;
% FID = fopen('vfile');
% status = fseek(FID,0,'bof');
% v=fread(FID,[200,257],'float');
% Pad data with zeros and apply FFT
T=zeros(ntfft-nt,nx);
P=[D;T];
CP=ntfft*ifft(P);
%----test----- % not correct! Still cannot figure out '+1'&'-1'
%CP=fft(P); % means what to the Fourier Kernel !!!
CP=CP(1:nw,:);
clear T;
% clear D;
% loops over depth
profile on -detail 'builtin'
iz=1;
tz=fz;
while iz<=nz
for ix=1:nx
cresult(ix,iz)=0.0;
for iw=1:nw
cresult(ix,iz)=cresult(ix,iz)+real(CP(iw,ix))/ntfft;
end
end
vmin=0;
for ix=1:nx
vmin=vmin+1.0/v(iz,ix)/nx;
end
vmin=1.0/vmin*0.5;
% ???----------------------------------
% data modified to fit SU-style FFT
CQ=zeros(nw,nk);
for ik=1:nx
if mod(ik,2)==1
CQ(:,ik)=CP(:,ik);
else
CQ(:,ik)=-CP(:,ik);
end
end
% CQ(:,1:nx)=CP;
% ???----------------------------------
% FFT to W-K domain
for iw=1:nw
CQ(iw,:)=fft(CQ(iw,:));
% CQ(iw,:)=fftshift(fft(CQ(iw,:))); % No difference will be shown ?!
end
v1=vmin;
%_________Phase Shift_________
k=fk;
ik=1;
while ik<=nk
w=fw;
iw=1;
while iw<=nw
%-------------------------------
if w==0.0
w=1.0e-10/dt;
end
kz1=1.0-(v1*k/w)^2;
if kz1>0
phase1=-w*sqrt(kz1)*dz/v1;
cshift1=cos(phase1)+sqrt(-1)*sin(phase1);
CQ1(ik,iw)=CQ(iw,ik)*cshift1;
else
cshift1=0.0;
CQ1(ik,iw)=CQ(iw,ik)*cshift1;
end
%--------------------------------
iw=iw+1;
w=w+dw;
end
ik=ik+1;
k=k+dk;
end
%______________________________
%-----FFT to W-X domain------
CQ1=ifft(CQ1);
%_________Time shift___________
for ix=1:nx
w=fw;
for iw=1:nw
if mod(ix,2)==0
CQ1(ix,iw)=-CQ1(ix,iw);
else
CQ1(ix,iw)=CQ1(ix,iw);
end
kz2=(1.0/v1-2.0/v(iz,ix))*w*dz;
cshift2=cos(kz2)+sqrt(-1)*sin(kz2);
CP(iw,ix)=CQ1(ix,iw)*cshift2;
w=w+dw;
end
end
%__________________________________
tz=tz+dz;
iz=iz+1
end
profile viewer
cresult=cresult';
%---------------------------PLOT-------------------------------
x=(0:(nx-1))*dx;
z=linspace(0,(nt-1)*dt,nz);
[xx zz]=meshgrid(x,z);
clf;
pcolor(xx,zz,cresult);
set(gca, 'YDir', 'reverse','XAxisLocation', 'top')
shading interp
colormap winter
评论0