inv_taup.old.m

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

function [m,M,Q] = inv_taup(d,dt,h,q,N,flow,fhigh,mu,type);
%INV_TAUP: Inverse Radon transform
%
%  Given a cmp gather this subroutine is used
%  to compute the  Radon panel. By Radon Transform
%  we mean the Linear or Parabolic Radon Transform.
%
%  [m] = inv_taup(d,dt,h,q,N,flow,fhigh,mu,type)
% 
%  IN   d: seismic traces ( d(nt,nh) 
%       dt: sampling in sec
%       h(nh) offset or position of traces in meters
%       q(nq) ray parameters  if N=1
%             residual moveout at far offset if N=2
%       N:1 Linear tau-p  
%        :2 Parabolic tau-p
%       flow:  freq.  where the inversion starts in HZ (> 0Hz)
%       fhigh: freq.  where the inversion ends in HZ   (< Nyquist) 
%       mu: regularization parameter 
%       type: 'ls' Least-squares solution
%             'hr' non-iterative high resolution 
%
%  OUT  m: the linear or parabolic tau-p panel [m(nt,nq)]
%
%  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,nh]= size(d);
nq = max(size(q));

if N==2; h=h/max(abs(h));end;  

nfft = 2^nextpow2(nt);

D = fft(d,nfft,1);
M = zeros(nfft,nq);
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

Q = eye(nq);
Qtr = nq;
Mtr = nq*nh;
beta = mu*Mtr/Qtr;


for ifreq=ilow:ihigh

 f = 2.*pi*(ifreq-1)/nfft/dt;
 L = exp(i*f*(h.^N)'*q); 
 y = D(ifreq,:)';
 x = L'*y;

 MATRIX = L'*L; 


 A = MATRIX + beta*Q;
 x = inv(A) *L'* y; 

 M(ifreq,:) = x';
 M(nfft+2-ifreq,:) = conj(x)';

if type=='hr'; 
   xx = x /max(abs(x)); Q = diag(1./(abs(xx).^2+0.0001)); 
end

end
 
M(nfft/2+1,:) = zeros(1,nq);
m = real(ifft(M,[],1));
m = m(1:nt,:);

return