Band pass Filter in seismic data

11 views (last 30 days)
Mila Apriani
Mila Apriani on 29 May 2020
Answered: Surya Talluri on 14 Aug 2020
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

Answers (1)

Surya Talluri
Surya Talluri on 14 Aug 2020
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:

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!