% Program Perencanaan dan Simulasi Filter BPF IIR

% Oleh : Rizal Dwi Prasetyo

clc; clear all; clf;

k1=-3; k2=-20;
f1=1400; fu=1700;
f1=480; f2=5000;
fs=16000; Ts=1/fs;
f=(0:255)/256*(fs/2);
t=0:Ts:0.03; % iterasi untuk sinyal input kontinyu
nn=0:Ts:0.03; % iterasi untuk sinyal input hasil sampling

% SINYAL INPUT xt
x=sin(2*pi*216*t)+sin(2*pi*1832*t)+sin(2*pi*7000*t);

% SINYAL INPUT xn
xnn=sin(2*pi*216*nn)+sin(2*pi*1832*nn)+sin(2*pi*7000*nn);
xn=xnn(1:200);

% LISTING PROGRAM UNTUK MENGHITUNG NILAI ORDE FILTER n

w1=(2*pi*f1)/fs;
w2=(2*pi*f2)/fs;
w1=(2*pi*f1)/fs;
wu=(2*pi*fu)/fs;
omega1=2/Ts*tan(w1/2);
omega2=2/Ts*tan(w2/2);
omega1=2/Ts*tan(w1/2);
omegau=2/Ts*tan(wu/2);
A=(-omega1^2+(omega1*omegau))/(omega1*(omegau-omega1));
B=(omega2^2-(omega1*omegau))/(omega2*(omegau-omega1));
omegar=min(abs(A),abs(B));
n=(log10((10^(-k1/10)-1)/(10^(-k2/10)-1)))/(2*log10(1/omegar));

% KOEFISIEN FILTER DARI HASIL PERHITUNGAN MANUAL

B=[0.0554 0 -0.554]; A=[1 -1.8044 1.0333];

% LISTING PROGRAM MENGHITUNG RESPON IMPULS FILTER hn TABELNYA

xx=zeros(50,1);
xx(1)=1;
for n=1:50;
if n==1;
y(n)=0.0554*xx(n);
elseif n==2;
y(n)=0.0554*xx(n)+1.8044*y(n-1);
elseif n>=3;
y(n)=0.0554*xx(n)-0.0554*xx(n-2)+1.8044*y(n-1)-1.0333*y(n-2);
end
h(n)=y(n);
end
fprintf(' Tabel Respon Impuls h(n)\n');
fprintf('---------------------------\n');
fprintf(' N h(n) \n');
fprintf('---------------------------\n');
for i=1:50;
fprintf('%5.0f \b',(i-1));
fprintf('%15.4f \n',h(i));
end
fprintf('---------------------------\n');

% LISTING PROGRAM MENGHITUNG RESPON FREKUENSI FILTER Hejw

[H,w]=freqz(B,A,100);

% LISTING PROGRAM MENGHITUNG MAGNITUDO SINYAL INPUT X(ejw)

xf=fft(xn,512); % Transformasi Fourier sinyal x
X=xf(1:256);

% LISTING PROGRAM MENGHITUNG RESPON IMPULS OUTPUT yn

for n=1:200;
if n==1;
y(n)=0.554*xn(n);
elseif n==2;
y(n)=0.0554*xn(n)+1.8044*y(n-1);
elseif n>=3;
y(n)=0.0554*xn(n)-0.0554*xn(n-2)+1.8044*y(n-1)-1.0333*y(n-2);
end
end

% LISTING PROGRAM MENGHITUNG MAGNITUDO SINYAL OUTPUT Yejw

yf=fft(y,512); % Transformasi Fourier sinyal output y
Y=yf(1:256);
hertz=w/(2*pi*Ts);

% LISTING PROGRAM UNTUK PLOT GRAFIK

figure(1);
plot(t,x);
title('Grafik Sinyal Input dalam domain waktu x(t)')
xlabel('Waktu t, (detik)')
ylabel('Amplitudo, X(f)')
grid on;
figure(2)
plot(f,abs(X));
title('Grafik Sinyal Input dalam domain frekuensi X(f)')
xlabel('Frekuensi, hertz')
ylabel('Amplitudo, X(f)')
grid on;
figure(3);
stem(1:200,xn);
title('Grafik Respon Impuls Sinyal Input xn')
xlabel('Sampling ke-n')
ylabel('Respon Impuls sinyal input xn')
axis([1 100 -2 2]);
grid on;
figure(4)
plot(f*2*Ts,abs(X));
title('Grafik Respon Magnitudo Sinyal Input dalam domain frekuensi, Xejw')
xlabel('Frekuensi Ternomalisasi, W (pi rad)')
ylabel('Magnitudo,[X]')
grid on;
figure(5);
stem(0:49,h);
title('Grafik Respon Impuls Filter h(n)')
xlabel('Sampling ke-n')
ylabel('Respon Impuls h(n)')
grid on;
figure(6)
plot(w/pi,abs(H));
title('Grafik Response Magnitudo, |Hejw|^2')
xlabel('Frekuensi Ternormalisasi, W (pi rad)')
ylabel('Magnitudo, |Hejw|^2 (pi rad)')
grid on;
figure(7)
plot(w/pi,20*log10(abs(H)));
title('Grafik Response Magnitudo, 20log|Hejw|')
xlabel('Frekuensi Ternomalisasi, W (pi rad)')
ylabel('Magnitudo, 20log|Hejw| (dB)')
grid on;
figure(8);
plot(w/pi,angle(H)/pi);
title('Respon fasa')
xlabel('Frekuensi Ternomalisasi, W (pi rad)')
ylabel('Phase (pi rad)')
grid on;
figure(9);
stem(1:200,y);
title('Grafik Respon Impuls Sinyal Output y(n)')
xlabel('Sampling ke-n')
ylabel('Respon Impuls y(n)')
axis([1 100 -1 1]);
grid on;
figure(10);
plot(f*2*Ts,abs(Y));
title('Grafik Respon Magnitudo Sinyal Output dalam domain frekuensi, Yejw')
xlabel('Frekuensi Ternomalisasi, W (pi rad)')
ylabel('Magnitudo, |Yejw|')
grid on;
figure(11);
plot(t(1:200),y);
title('Grafik Sinyal Output dalam domain waktu y(t)')
xlabel('Waktu t, (detik)');
ylabel('Magnitudo, y(t))')
axis([0 0.01 -1 1]);
grid on;
figure(12);
plot(f,abs(Y));
title('Grafik Sinyal Output dalam domain frekuensi Y(f)')
xlabel('Frekuensi, hertz');
ylabel('Magnitudo, Y(f))')
grid on;