bp_filter.rar_bp filter_bp_filter_地震_地震 滤波_滤波 地震

function [o] = bp_filter(d,dt,f1,f2,f3,f4); %BP_Filter: Apply a band-pass filter to a group of traces. % % [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) % % Example: % % d=linear_events; % dout = bp_filter(d,0.004,1,3,30,40); % wigb([d,dout]); % % Copyright (C) 2008, Signal Analysis and Imaging Group % For more information: http://www-geo.phys.ualberta.ca/saig/SeismicLab % Author: M.D.Sacchi % % This program is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published % by the Free Software Foundation, either version 3 of the License, or % any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details: http://www.gnu.org/licenses/ % [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,:));