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