Code covered by the BSD License

# Seis_Pick

### James Verdon (view profile)

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
```