Code covered by the BSD License  

Highlights from
Seis_Pick

image thumbnail

Seis_Pick

by

 

06 Feb 2012 (Updated )

Seis_Pick provides an interactive picking environment for processing seismic data.

varargout=predictive_filter(inSeis, NoisemaxWin, varargin)
function varargout=predictive_filter(inSeis, NoisemaxWin, varargin)
%
% function varargout=predictive_filter(inSeis, NoisemaxWin, varargin)
%
% e.g., filteredtrace=predictive_filter(trace,1000,0.001, [2 2 1 .5])
%
% Uses pre-event noise to create a notch filter, which is then used to 
% filter the data. This is an effective way of cleaning up signals
%
% Written by A. Wustefeld, University of Bristol, 2009 - 2011
%


if nargout==2 && nargin ~=3
    error ('Need sampling rate as third input argument!')
end

weight=varargin{2};


Amp = inSeis;
len = length(Amp);
o   = nextpow2(len);
warning off all
L   = 2^o;  % L extends the trace to the next power of 2 

Amp(end+1:L)   = 0;
Noise          = Amp(1:NoisemaxWin); % Take the first part of the trace as noise
Noise(end+1:L) = 0;  % Make the rest of the noise trace to be zeros


len2  = round(len*.1); %taper length is 1% of total seismogram length
nn   = 1:len2;
nn2  = (len-len2+1):len;
x    = linspace(pi, 2*pi, len2);
taper  = 0.5 * (cos(x')+1);
taper2 = flipud(taper);

%mPA   = 0;
loops = 0;
%oldstd=1;
FF=nan((2^(o-1))+1,3);
for k = 1:length(weight)
    % taper at begin             taper at end of seismogram
    Amp(nn) = Amp(nn).*taper;    Amp(nn2) = Amp(nn2).*taper2;


    Afft = fft(Amp,2^o);   % FFT the whole trace
%     PA   = Afft.* conj(Afft);

    x=Noise;
    %power spectrum of the noise cross-correlation
    Yfft = fft(x,2^o);
    Px   = Yfft.* conj(Yfft) / 2^o;

    loops=loops+1;
    if max(Px)==0
        mm = 1;
    else
        mm = max(Px);  %Normalise Px by the maximum
    end
    Px=Px/mm;


    Fil = 1 - (Px)*weight(k); % The notch is the inverse of the power

    Fil(Fil<0.001)=0.001;  %Limit Fil to between 0-1
    Fil(Fil>=1)=1;

    % Filtered signal
    Amp   = ifft(Afft.*Fil(:) );
    Noise = Amp(1:NoisemaxWin);
    Noise(end+1:L)=0;
    if nargout~=1
        FF(:,loops)=Fil(1:((2^(o-1))+1));
    end

end

warning on all

outsignal = Amp(1:len);
nargout;

switch nargout
    case 0
        t=1:len;
        subplot(2,1,1)
        plot(t,inSeis+max(abs(inSeis))-mean(inSeis),'r',  t,outsignal-max(abs(outsignal))-mean(outsignal),'k')
        %        xlim([200 600])
        set(gca,'ytick',[])
        xlabel('Samples')
        legend({'Original','Filtered'})

        PI = Afft.* conj(Afft);


        subplot(2,1,2)
        Afft = fft(Amp,2^o);
        PA = Afft.* conj(Afft);


        Amp = inSeis;
        Amp(end+1:L) = 0;
        Y   = fft(Amp,2^o);
        PS  = Y.* conj(Y);


        if nargin==3
            smpl=varargin{1};
        else
            smpl=6000;
        end
        f = smpl*(0:2^(o-1))/2^o;
        Power = .9*[PS(1:((2^(o-1))+1))/max(PS(1:((2^(o-1))+1))), PI(1:((2^(o-1))+1))/max(PI(1:((2^(o-1))+1))), PA(1:((2^(o-1))+1))/max(PA(1:((2^(o-1))+1)))];
        plot(  f, Power(:,1), 'r-',f, Power(:,2) , 'b-',f, Power(:,3) , 'k-')

        hold on
        plot( f, 1+FF(:,1),'r-',  f, 1+FF(:,end),'b-' )
        line([0 max(xlim)] ,[1 1],'color','k')
        hold off
        set(gca,'Ytick', [0 .9 1 2],'Yticklabel' ,'0|1|0|1')
        legend({'Initial Spectrum','First run','Second Run'})
        xlabel('Frequency [Hz]')


    case 1
        varargout{1} = outsignal;
    case 2

        smpl=varargin{1};
        f = smpl*(0:2^(o-1))/2^o;
        varargout{1} = outsignal;
        varargout{2} = [f', sum(FF,2)-loops+1];
end

Contact us