www.pudn.com > fk_migration.rar > fkmigration.m, change:2016-03-08,size:3374b

function [migratedsection_zx]=fkmigration(section_tx,dt,dx,v)
%进行fk偏移
%section_tx二维数组
%常数dt，时间采样间隔
%向量t
%常数dx，空间采样间隔
%向量x
%向量v，
if(nargin<4)
error('函数输入参数不足');
end
%记时开始
timestart=clock;
%开始进行fk偏移
disp('fkmigration.m程序开始运行......');
[nt,nx]=size(section_tx);
%补零，补成原先的一倍，并判断道数是否为2^n.如果不是补成2^n道.
new_nx=2*nx;
n=0;
while(2^n<new_nx)
n=n+1;
end
new_nx=2^n;
new_nt=2*nt;
n=0;
while(2^n<new_nt)
n=n+1;
end
new_nt=2^n;
section_tx=[section_tx;zeros(new_nt-nt,nx)];
section_tx=[section_tx,zeros(new_nt,new_nx-nx)];
disp(['补零已完成，t方向补零' num2str(new_nt-nt) '个。x方向补零' num2str(new_nx-nx) '个。']);
%做傅氏变换和反傅氏变换
InterpoPoint=2;
CArray=1-ceil(InterpoPoint/2):ceil(InterpoPoint/2);
spectrum=fft(section_tx);
spectrum=fft(spectrum.').';
df=2*pi/dt/new_nt;
nf=new_nt;
nf_positive=floor(1+new_nt/2);
dkz=df/v;
f=[(0:nf_positive-1)*df];
kz=f/v;
spectrum=spectrum(1:nf_positive,:);
dkx=2*pi/dx/new_nx;
nkx_positive=floor(1+new_nx/2);
kx=[(0:nkx_positive-1)*dkx];
new_nxadd2=new_nx+2;
for(kxi=2:nkx_positive)
tempf1=zeros(1,nf_positive).';
tempf2=tempf1;
for(kzi=1:nf_positive)
tempf=v*sqrt(kz(kzi)^2+kx(kxi)^2);
factor=f(kzi)*v/tempf;
%%%插值
n=tempf/df+1;
n1=floor(n);
n2=ceil(n);
if(n1==n2)
tempf1(kzi)=factor*spectrum(n,kxi);
tempf2(kzi)=factor*spectrum(n,new_nxadd2-kxi);
else

distance1=n-n1;

IP=zeros(1,InterpoPoint);
IN=IP;
for(ii=1:InterpoPoint)
flag=n1+CArray(ii);
if(flag>=1&&flag<=nf_positive)
IP(ii)=spectrum(flag,kxi);
IN(ii)=spectrum(flag,new_nxadd2-kxi);
end
end
C1=distance1-CArray;
IP=IP./C1;
IN=IN./C1;
for(nn=1:InterpoPoint)

C2(nn)=sin(C1(nn)*pi)*exp(-j*C1(nn)*pi)/pi;
end
spectrum_i_m_ikxP=sum(C2.*IP);
spectrum_i_m_ikxN=sum(C2.*IN);
tempf1(kzi)=factor*spectrum_i_m_ikxP;
tempf2(kzi)=factor*spectrum_i_m_ikxN;

end
end
spectrum(:,kxi)=tempf1;
spectrum(:,new_nxadd2-kxi)=tempf2;
end
spectrum(:,1)=spectrum(:,1)*v;
spectrum=ifft(spectrum.').';

spectrum=[spectrum;conj(flipud(spectrum(2:(new_nt-nf_positive+1),:)))];
if(nf_positive==new_nt/2+1);
spectrum(nf_positive,:)=0;
end
migratedsection_zx=ifft(spectrum);
migratedsection_zx(:,(nx+1):end)=[];
migratedsection_zx(nt+1:end,:)=[];
tend=etime(clock,timestart);
disp(['程序结束，总共用时:' num2str(tend) '秒。']);%显示用时
function dout=attenuation(l1,l0)%用cos函数计算从1到0的一个衰减值
%l1...靠近1那段的长度
%l0...靠近0那段的长度
if(l1==0&&l0==0)
error('l1和l0不能同时为零');
else
dout=0.5*(1+cos(l1*pi/(l1+l0)));
end