bp_filter.m

GeoPhysics at the University of Alberta 提供的一个基于Matlab的地震处理软件

function  [o] =  bp_filter(d,dt,f1,f2,f3,f4);
% BP_Filter: apply a BP filter with corner freq. given
%            ny f1,f2,f3,f4. FFT construction (not good
%            for short time series)
% 
%  [o] = bp_filter(d,dt,f1,f2,f3,f4);
%
%  IN     d:    data (columns are traces)
%         dt:   sampling interval in sec
%         f1:   freq. in Hz
%         f2:   freq. in Hz
%         f3:   freq. in Hz
%         f4:   freq. in Hz
%
%   ^
%   |     ___________
%   |    /           \   Amplitude spectrum
%   |   /             \
%   |  /               \
%   |------------------------>
%      f1 f2        f3 f4
%
%  OUT    o:    output  (columns are traces)
%
%
%  Author(s): M.D.Sacchi ([e-mail])
%  Copyright 1998-2003 SeismicLab
%  Revision: 1.2  Date: Dec/2002
%
%  Signal Analysis and Imaging Group (SAIG)
%  Department of Physics, UofA
%


 [nt,nx] = size(d);
 k = nextpow2(nt);
 nf = 4*(2^k);

 i1 = floor(nf*f1*dt)+1;
 i2 = floor(nf*f2*dt)+1;
 i3 = floor(nf*f3*dt)+1;
 i4 = floor(nf*f4*dt)+1;

 up =  (1:1:(i2-i1))/(i2-i1);
 down = (i4-i3:-1:1)/(i4-i3);
 aux = [zeros(1,i1), up, ones(1,i3-i2), down, zeros(1,nf/2+1-i4) ];
 aux2 = fliplr(aux(1,2:nf/2));

 c = 0; % zero phase (could apply rotations as well)
 F = ([aux,aux2]');
 Phase = (pi/180.)*[0.,-c*ones(1,nf/2-1),0.,c*ones(1,nf/2-1)];
 Transfer = F.*exp(-i*Phase');


 D = fft(d,nf,1);

 for k = 1:nx
  Do(:,k) = Transfer.*D(:,k);
 end

 o = ifft(Do,nf,1);

 o = real(o(1:nt,:));