% Cebisevljev filtar PO 
% PO:  0.4*pi <= w <= 0.6*pi
%      Rp = 3dB
% NPO: w <= 0.3*pi
%      w >= 0.7*pi
%      Rs = 30dB
%
% Projektovanje korak po korak (bilinearna transformacija)
%
wp1 = 0.4*pi;
wp2 = 0.6*pi;
ws1 = 0.3*pi;
ws2 = 0.7*pi;
Rp = 3;
Rs = 30;

%% Preslikavanje frekvencija u analogne
dt = 1;

Wp1 = 2/dt*tan(wp1/2)
Wp2 = 2/dt*tan(wp2/2)
Ws1 = 2/dt*tan(ws1/2)
Ws2 = 2/dt*tan(ws2/2)

%% centralna ucestanost i propusni opseg
W0 = sqrt(Wp1 * Wp2)
B = Wp2 - Wp1
Q = W0 / B

%% normalizacija
Wp1n = Wp1/W0;
Wp2n = Wp2/W0;
Ws1n = Ws1/W0;
Ws2n = Ws2/W0;

%% frekvencijska transformacija NP->PO
Wpp1 = Q*(Wp1n^2-1)/Wp1n
Wpp2 = Q*(Wp2n^2-1)/Wp2n
Wss1 = Q*(Ws1n^2-1)/Ws1n
Wss2 = Q*(Ws2n^2-1)/Ws2n

%% projektovanje filtra
% izabrali smo strozije uslove
[n, wn] = cheb1ord(Wpp2, Wss2, Rp, Rs, 's')
[bnp, anp] = cheby1(n, Rp, wn, 's')
[bpo, apo] = lp2bp(bnp, anp, W0, B)
[b, a] = bilinear(bpo, apo, 1/dt)
figure, freqz(b, a)
20*log10(abs(freqz(b, a, [0.3 0.4 0.6 0.7]*pi)))

%% projektovanje filtra direktno
[n1, wn1] = cheb1ord([wp1 wp2]/pi, [ws1 ws2]/pi, Rp, Rs)
[b1, a1] = cheby1(n, Rp, wn1)
figure, freqz(b1, a1)
20*log10(abs(freqz(b1, a1, [0.3 0.4 0.6 0.7]*pi)))

%% impulsna invarijansa

wp1 = 0.4*pi;
wp2 = 0.6*pi;
ws1 = 0.3*pi;
ws2 = 0.7*pi;
Rp = 3;
Rs = 30;

%% Preslikavanje frekvencija u analogne
dt = 1;

Wp1 = wp1 / dt
Wp2 = wp2 / dt
Ws1 = ws1 / dt
Ws2 = ws2 / dt

%% centralna ucestanost i propusni opseg
W0 = sqrt(Wp1 * Wp2)
B = Wp2 - Wp1
Q = W0 / B;

%% normalizacija
Wp1n = Wp1/W0;
Wp2n = Wp2/W0;
Ws1n = Ws1/W0;
Ws2n = Ws2/W0;

%% frekvencijska transformacija NP->PO
Wpp1 = Q*(Wp1n^2-1)/Wp1n
Wpp2 = Q*(Wp2n^2-1)/Wp2n
Wss1 = Q*(Ws1n^2-1)/Ws1n
Wss2 = Q*(Ws2n^2-1)/Ws2n

%% projektovanje filtra
% izabrali smo strozije uslove
[n, wn] = cheb1ord(Wpp2, Wss2, Rp, Rs, 's')
[bnp, anp] = cheby1(n, Rp, wn, 's')
[bpo, apo] = lp2bp(bnp, anp, W0, B)
[b, a] = impinvar(bpo, apo, 1/dt)
figure, freqz(b, a)
20*log10(abs(freqz(b, a, [0.3 0.4 0.6 0.7]*pi)))

%% impulsna invarijansa (jos jedan primjer u kojem se ukazuje na probleme)

wp1 = 0.4*pi;
wp2 = 0.6*pi;
ws1 = 0.3*pi;
ws2 = 0.8*pi;
Rp = 3;
Rs = 20;

%% Preslikavanje frekvencija u analogne
dt = 0.001;

Wp1 = wp1 / dt
Wp2 = wp2 / dt
Ws1 = ws1 / dt
Ws2 = ws2 / dt

%% centralna ucestanost i propusni opseg
W0 = sqrt(Wp1 * Wp2)
B = Wp2 - Wp1
Q = W0 / B;

%% normalizacija
Wp1n = Wp1/W0;
Wp2n = Wp2/W0;
Ws1n = Ws1/W0;
Ws2n = Ws2/W0;

%% frekvencijska transformacija NP->PO
Wpp1 = Q*(Wp1n^2-1)/Wp1n
Wpp2 = Q*(Wp2n^2-1)/Wp2n
Wss1 = Q*(Ws1n^2-1)/Ws1n
Wss2 = Q*(Ws2n^2-1)/Ws2n

%% projektovanje filtra
% izabrali smo strozije uslove
[n, wn] = cheb1ord(Wpp2, Wss2, Rp, Rs, 's')
[bnp, anp] = cheby1(n, Rp, wn, 's')
[bpo, apo] = lp2bp(bnp, anp, W0, B)
[b, a] = impinvar(bpo, apo, 1/dt)
figure, freqz(b, a)
20*log10(abs(freqz(b, a, [0.3 0.4 0.6 0.8]*pi)))

%% niskopropusni Batervortov filtar
% PO:  Fp = 1kHz
%      Rp = 1dB
% NPO: Fs = 2kHz
%      Rs = 30dB
%

Fp = 1000;
Fs = 2000;
Rp = 1;
Rs = 30;

%% analogni filtar
[n, wn] = buttord(2*pi*Fp, 2*pi*Fs, Rp, Rs, 's')
e = sqrt(10^0.1*Rp - 1)
wc = e^(-1/n);
[ba, aa] = butter(n, 2*pi*Fp*wc, 's')
20*log10(abs(freqs(ba, aa, 2*pi*[0 Fp Fp*wc Fs])))

%% impulsna invarijansa (pogresno!)
Fo = 4000;

[b, a] = impinvar(ba, aa, Fo);

20*log10(abs(freqz(b, a, [1000 2000], Fo)))

%% impulsna invarijansa (bolje!)
Fo = 6000;

[b, a] = impinvar(ba, aa, Fo);

20*log10(abs(freqz(b, a, [1000 2000], Fo)))

