for_taup.m

著名的seismiclab的代码 是地震学研究人员必备的工具

function [d]=for_taup(m,dt,h,p,N,flow,fhigh);
%INV_TAUP: Forward  Radon transform. Given Radon panel (tau-p), 
%          this function computes the data in time-space
%
%  [d] = for_taup(m,dt,h,p,N,flow,fhigh);
%
%  IN   m:     the Radon panel (d(nt,np)
%       dt:    sampling in sec
%       h(nh): offset or position of traces in mts
%       p(np): ray parameter  to retrieve if N=1
%              curvature of the parabola if N=2
%       N:     N=1 linear tau-p
%              N=2 parabolic tau-p
%       flow:  min freq. in Hz
%       fhig:  max freq. in Hz
%
%  OUT  d:     data
%
%
%  SeismicLab
%  Version 1
%
%  written by M.D.Sacchi, last modified December 10,  1998.
%  [e-mail]
%
%  Copyright (C) 1998 Signal Analysis and Imaging Group
%                     Department of Physics
%                     The University of Alberta
%

[nt,nq] = size(m);
nh = length(h);
if N==2; h=h/max(abs(h));end;   
nfft = 2*(2^nextpow2(nt));

M = fft(m,nfft,1);
D = zeros(nfft,nh);
i = sqrt(-1);

ilow  = floor(flow*dt*nfft)+1; if ilow<1; ilow=1;end;
ihigh = floor(fhigh*dt*nfft)+1;
if ihigh>floor(nfft/2)+1; ihigh=floor(nfft/2)+1;end
 
for ifreq=ilow:ihigh
f = 2.*pi*(ifreq-1)/nfft/dt;
L = exp(i*f*(h.^N)'*p);
x = M(ifreq,:)';
y = L * x; 
D(ifreq,:) = y';
D(nfft+2-ifreq,:) = conj(y)';
end
D(nfft/2+1,:) = zeros(1,nh);
d = real(ifft(D,[],1));
d = d(1:nt,:);

return;