
clear
close all

% f-domain representation of Butterworth filter


f=0:100;
f0=20;
fs=15;



for iplot=1:4

subplot(2,2,iplot )
% n=1
n=2*iplot;
df=5
Fb = 1./(( 1+((f-f0)/df).^(2*n)));
plot(f,Fb)
set(gca,'FontSize',fs)
title(sprintf(' Butterworth n=%i, f0=%i Hz ',n,f0))
xlabel(' Frequency (Hz) ')
ylabel(' Filter amplitude ')



end

% effect on a spike

dt=.005;
t=(1:1000)*dt;
s=t*0;s(100)=1;

figure

subplot(2,2,1)
plot(t,s)
set(gca,'FontSize',fs)
title(sprintf(' Original function ',n,f0))
xlabel(' Time (s) ')
ylabel(' Amplitude ')


for iplot=2:4

n=(iplot-1)^2;
subplot(2,2,iplot )
Fs=bandpass(s,dt,f0-df,f0+df,n);
plot(t,Fs)
set(gca,'FontSize',fs)
title(sprintf(' Filtered with n=%i, f0=%i Hz ',n,f0))
xlabel(' Time (s) ')
ylabel(' Amplitude ')
axis([ 0 1 min(Fs) max(Fs) ])

end


% effect on a seismogram


df=.5;

load germany

t=translation_time;
dt=.05;
f0=1;
s=translation_Z;

figure

subplot(2,2,1)
plot(t,s)
set(gca,'FontSize',fs)
title(sprintf(' Original function ',n,f0))
xlabel(' Time (s) ')
ylabel(' Amplitude ')
t1=40; t2=65;
rmax=max( abs(s(floor(t1/dt):floor(t2/dt))));
axis([ 40 65  -rmax rmax ])

for iplot=2:4

n=(iplot);
subplot(2,2,iplot )
Fs=bandpass(s,dt,f0-df,f0+df,n);
plot(t,Fs)
set(gca,'FontSize',fs)
title(sprintf(' Filtered with n=%i, f0=%i Hz ',n,f0))
xlabel(' Time (s) ')
ylabel(' Amplitude ')
t1=40; t2=65;
rmax=max( abs(Fs(floor(t1/dt):floor(t2/dt))))
axis([ 40 65  -rmax rmax ])

end


figure

subplot(2,2,1)
plot(t,s)
set(gca,'FontSize',fs)
title(sprintf(' Original function ',n,f0))
xlabel(' Time (s) ')
ylabel(' Amplitude ')
t1=50; t2=60;
rmax=max( abs(s(floor(t1/dt):floor(t2/dt))));
axis([ 50 60  -rmax rmax ])

for iplot=2:4

n=4;
f0=iplot
subplot(2,2,iplot )
Fs=bandpass(s,dt,f0-df,f0+df,n);
plot(t,Fs)
set(gca,'FontSize',fs)
title(sprintf(' Filtered with n=%i, f0=%g Hz ',n,f0))
xlabel(' Time (s) ')
ylabel(' Amplitude ')
t1=50; t2=60;
rmax=max( abs(Fs(floor(t1/dt):floor(t2/dt))))
axis([ t1 t2  -rmax rmax ])

end