Band pass Filter in seismic data

Hello everyone,

Can i ask this question to all of you,

here i attach seismic file in .txt that contain the amplitude of signal earthquake, i use it to analyze with band pass filter, but i found it wrong with the matlab code, please help me:(

the matlab code:

% Input Sinyal
fid = fopen('LEM_BHN.TXT');% Open data
data = textscan(fid,'%f'); % Read data file
fclose(fid); % Close data
x = data{1}; % Amplitudo signal in time domain
N = length(x); % Length sample
dt = 0.01 % Delta time
t = [0:dt:(N-1)*dt]; % Time series

% Spectral analysis using FFT
y = data(:,N);
m = sum(y,2);
Fs = 250;
LD = length(m);
YD = fft(m);
P2D = abs(YD/LD);
P1D = P2D(1:LD/2+1);
fD = Fs*(0:(LD/2))/LD;
figure(1)
plot(fD,P1D)
hold on
title ('Analysis Spectrum')
ylabel ('Amplitude')
xlabel ('Frequency (Hz)')
axis([min(fD) max(fD) min(P1D) max(P1D)])
 
% m and n top and bottom of range limit bandpass
% Syarat m < n 
m = input('Input the bottom freq value = ');
n = input('Input the top freq value = ');
 
% Plot Raw Data
figure(2)
for i = 1:N
    subplot(1,2,1)
    plot(data(:,i)+300*i,t,'k');
    set(gca,'Ydir','reverse');
    title ('Before Bandpass Filter')
    legend('Raw')
    ylabel('Time (s)')
    xlabel('Offset (m)')
    hold on
    axis([min(data(:,1)+300*1) max(data(:,i)+300*i) min(t) max(t)])
end
 
% Bandpass filter and plotting filter
for i = 1:N
    y = data(:,i);
    Fs = 250; %Sampling Freq
    [b,a] = butter(2,[m n]/(Fs/2)); %Bandpass butter function
    yabp = filter(b,a,y);%the result of Bandpass 
    subplot(1,2,2)
    plot(yabp+300*i,t,'k')%Plot Bandpass Filter
    set(gca,'Ydir','reverse');
    legend('band')
    title ('After Bandpass Filter')
    legend('Bandpass')
    ylabel ('Time (s)')
    xlabel ('Offset (m)')
    hold on
    axis([min(yabp+300*1) max(yabp+300*i) min(t) max(t)])
end

Answer 1

[Surya Talluri]

I understand that the data provided is for a single sample and code is for multiple samples. As you are plotting the signal and filter in loop on same figure, plot is overriding the existing plot.

[b,a] = butter(2,[m n]/125, 'bandpass');
xfilt = filter(b,a,x)
 
subplot(2,1,1)
plot(t, x)
title ('Before Bandpass Filter')
  
subplot(2,1,2)
plot(t,xfilt)
title ('After Bandpass Filter')

You can observe the response of filter with “fvtool” function

fvtool(b,a)

We can use “designfilt” function to design a filter by mentioning even more parameters of filter. You can go for higher order filters, if you need sharper response.

d = designfilt('bandpassiir', 'FilterOrder', 2,...
     'HalfPowerFrequency1',m,'HalfPowerFrequency2',n,...
     'DesignMethod','butter','SampleRate',250);
  
fvtool(d)

xfilt = filter(d, x);

You can refer to following documentation for further understanding:

• butter - https://www.mathworks.com/help/signal/ref/butter.html
• subplot - https://www.mathworks.com/help/matlab/ref/subplot.html