fx.m

地震勘探的必备工具。大家相互交流

function [DATA_f] = fx(DATA,lf,mu,flow,fhigh,dt,type);
%FX: SNR enhancement using FX-AR modelling
%
%  FX processing using Canales' method. 
%
%  [DATA_f] = fx(DATA,lf,mu,flow,fhigh,dt,itype)
% 
%  IN    DATA: the data, columns are traces
%        lf: lenght of the ar process (lenght of the peo)
%        flow, fhigh: min and max freq. in the data in Hz
%        dt: sampling interval in sec
%        mu  : pre-whitening in 
%        type: 1 Forward, -1 Backward, 0 Transient free matrix
%              is used to design the PEO.
% 
%  OUT   DATA_f: filtered data
%       
%  NOTE: Columns of D are traces 
% 
%  Example:
%
%  nt = 200; nx=60; 
%  D = flat(nt,nx,40,0.004,0.1); D_f = fx(D,10,0.01,2,120,0.004,1;
%  subplot(211); wigb(D);   title('Before FX filtering')
%  subplot(212); wigb(D_f); title('After FX filtering') 
%
%
%  SeismicLab
%  Version 1
%
%  written by M.D.Sacchi, last modified December 10,  1998.
%  [e-mail]
%
%  Copyright (C) 1998 Seismic Processing and Imaging Group
%                     Department of Physics
%                     The University of Alberta
%

[nt,ntraces] = size(DATA);

nf = 2^nextpow2(nt);
 
DATA_FX_f = zeros(nf,ntraces);

% Lower and Upper samples of the DFT.

ilow  = floor(flow*dt*nf)+1; if ilow<1; ilow=1;end;
ihigh = floor(fhigh*dt*nf)+1;if ihigh>floor(nf/2)+1; ihigh=floor(nf/2)+1;end

% Go into FX

DATA_FX = fft(DATA,nf,1);

for if=ilow:ihigh;
  aux_in  = DATA_FX(if,:)';
  aux_out = ar_f_b(aux_in,lf,mu,type);
  DATA_FX_f(if,:) = aux_out';
end;
  for if=nf/2+2:nf
  DATA_FX_f(if,:) = conj(DATA_FX_f(nf-if+2,:));
end

% Back to TX (the output) 

DATA_f = real(ifft(DATA_FX_f,[],1));
DATA_f = DATA_f(1:nt,:);

return