% Demonstracija zadataka iz poglavlja o aproksimacionim metodima sa
% racunskih vjezbi

%% 2. zadatak
% Batervort
% Prvo normalizovani filtar
[n, wn] = buttord(1, 1.6, 3, 23, 's');
[bn, an] = butter(n, 1, 's');
[hn, w] = freqs(bn, an);
figure, plot(w, abs(hn));
axis([0 5 -0.2 1.2])
title('Amplitudna karakteristika normalizovanog filtra');
xlabel('\omega')
ylabel('|H(j\omega)|');
pause
figure, semilogx(w, 20*log10(abs(hn)));
title('Funkcija pojacanja normalizovanog filtra')
xlabel('\omega (log)')
ylabel('20*log10(|H(j\omega)|) [dB]')
pause
figure, semilogx(w, angle(hn));
title('Fazna karakteristika normalizovanog filtra')
xlabel('\omega (log)')
ylabel('radijana')
pause

%% Sada prava stvar
w = 2*pi*(0:10:100e3);
[n, wn] = buttord(2*pi*1200, 2*pi*1920, 3, 23, 's');
[b, a] = butter(n, 2*pi*1200, 's');
h = freqs(b, a, w);
figure, plot(w/2/pi, abs(h));
axis([0 3000 -0.2 1.2])
title('Amplitudna karakteristika filtra');
xlabel('f [Hz]')
ylabel('|H(j2\pif)|');
pause
figure, semilogx(w/2/pi, 20*log10(abs(h)));
axis([10 100e3 -250 20]);
title('Funkcija pojacanja filtra')
xlabel('f (log)')
ylabel('20*log10(|H(j2\pif)|) [dB]')
pause
figure, semilogx(w/2/pi, angle(h));
axis([10 100e3 -4 4]);
title('Fazna karakteristika filtra')
xlabel('f (log)')
ylabel('radijana')
pause

%% 3. zadatak
% Maksimalno ravna magnituda
% Varijanta u skladu sa teorijom
Wp = 2*pi*1200;
Ws = 2*pi*1920;
Ap = 0.5;
As = 23;
[n, wn] = buttord(Wp, Ws, Ap, As, 's');
wc = Wp*(sqrt(10^(0.1*Ap) - 1))^(-1/n);
[b, a] = butter(n, wc, 's');
h = freqs(b, a, w);
figure, plot(w/2/pi, abs(h));
axis([0 3000 -0.2 1.2])
title('Amplitudna karakteristika filtra');
xlabel('f [Hz]')
ylabel('|H(j2\pif)|');
pause
figure, semilogx(w/2/pi, 20*log10(abs(h)));
axis([10 100e3 -350 20]);
title('Funkcija pojacanja filtra')
xlabel('f (log)')
ylabel('20*log10(|H(j2\pif)|) [dB]')
pause
figure, semilogx(w/2/pi, angle(h));
axis([10 100e3 -4 4]);
title('Fazna karakteristika filtra')
xlabel('f (log)')
ylabel('radijana')
pause

%% 3. zadatak
% Maksimalno ravna magnituda
% MATLAB varijanta
[n, wn] = buttord(Wp, Ws, Ap, As, 's');
[b, a] = butter(n, wn, 's');
h = freqs(b, a, w);
figure, plot(w/2/pi, abs(h));
axis([0 3000 -0.2 1.2])
title('Amplitudna karakteristika filtra');
xlabel('f [Hz]')
ylabel('|H(j2\pif)|');
pause
figure, semilogx(w/2/pi, 20*log10(abs(h)));
axis([10 100e3 -350 20]);
title('Funkcija pojacanja filtra')
xlabel('f (log)')
ylabel('20*log10(|H(j2\pif)|) [dB]')
pause
figure, semilogx(w/2/pi, angle(h));
axis([10 100e3 -4 4]);
title('Fazna karakteristika filtra')
xlabel('f (log)')
ylabel('radijana')
pause

%% 5. zadatak
% Cebisevljev tip 1
% zahtjevi isti kao za prethodni
[n, wn] = cheb1ord(Wp, Ws, Ap, As, 's');
[b, a] = cheby1(n, Ap, wn, 's');
h = freqs(b, a, w);
figure, plot(w/2/pi, abs(h));
axis([0 3000 -0.2 1.2])
title('Amplitudna karakteristika filtra');
xlabel('f [Hz]')
ylabel('|H(j2\pif)|');
pause
figure, semilogx(w/2/pi, 20*log10(abs(h)));
axis([10 100e3 -350 20]);
title('Funkcija pojacanja filtra')
xlabel('f (log)')
ylabel('20*log10(|H(j2\pif)|) [dB]')
pause
figure, semilogx(w/2/pi, angle(h));
axis([10 100e3 -4 4]);
title('Fazna karakteristika filtra')
xlabel('f (log)')
ylabel('radijana')
pause

%% 5. zadatak
% Cebisevljev tip 2
% zahtjevi isti kao za prethodni
[n, wn] = cheb2ord(Wp, Ws, Ap, As, 's');
[b, a] = cheby2(n, As, wn, 's');
h = freqs(b, a, w);
figure, plot(w/2/pi, abs(h));
axis([0 5000 -0.2 1.2])
title('Amplitudna karakteristika filtra');
xlabel('f [Hz]')
ylabel('|H(j2\pif)|');
pause
figure, semilogx(w/2/pi, 20*log10(abs(h)));
axis([10 100e3 -100 20]);
title('Funkcija pojacanja filtra')
xlabel('f (log)')
ylabel('20*log10(|H(j2\pif)|) [dB]')
pause
figure, semilogx(w/2/pi, angle(h));
axis([10 100e3 -4 4]);
title('Fazna karakteristika filtra')
xlabel('f (log)')
ylabel('radijana')
pause

%% 6. zadatak
% Elipticki
% zahtjevi isti kao za prethodni
[n, wn] = ellipord(Wp, Ws, Ap, As, 's');
[b, a] = ellip(n, Ap, As, wn, 's');
h = freqs(b, a, w);
figure, plot(w/2/pi, abs(h));
axis([0 5000 -0.2 1.2])
title('Amplitudna karakteristika filtra');
xlabel('f [Hz]')
ylabel('|H(j2\pif)|');
pause
figure, semilogx(w/2/pi, 20*log10(abs(h)));
axis([10 100e3 -100 20]);
title('Funkcija pojacanja filtra')
xlabel('f (log)')
ylabel('20*log10(|H(j2\pif)|) [dB]')
pause
figure, semilogx(w/2/pi, angle(h));
axis([10 100e3 -4 4]);
title('Fazna karakteristika filtra')
xlabel('f (log)')
ylabel('radijana')
pause

%% 7. zadatak
% Besel-Tompsonov filtar
w = 2*pi*(0:10:1000e3);
W2 = 2*pi*6800;
W1 = 2*pi*9000;
Ap = 1.5;
dt = 1;
tau0 = 45e-6;
w1 = W1*tau0;
w2 = W2*tau0;
n = besselord(w1, dt, w2, Ap);
b0 = factorial(2*n)/(2^n * factorial(n));
wn = b0^(1/n);
[b, a] = besself(n, wn);
[b, a] = lp2lp(b, a, 1/tau0);
h = freqs(b, a, w);
figure, plot(w/2/pi, abs(h));
axis([0 50000 0 1])
title('Amplitudna karakteristika filtra');
xlabel('f [Hz]')
ylabel('|H(j2\pif)|');
pause
figure, semilogx(w/2/pi, 20*log10(abs(h)));
axis([10 1000e3 -250 20]);
title('Funkcija pojacanja filtra')
xlabel('f (log)')
ylabel('20*log10(|H(j2\pif)|) [dB]')
pause
figure, semilogx(w/2/pi, angle(h));
axis([10 1000e3 -4 4]);
title('Fazna karakteristika filtra')
xlabel('f (log)')
ylabel('radijana')
pause
s = j*w;
ad = polyder(a);
tau = 0.5*(polyval(ad, s)./polyval(a, s) + polyval(ad, -s)./polyval(a, -s));
tau = real(tau);
figure, semilogx(w/2/pi, 1e6*tau);
title('Grupno kasnjenje')
xlabel('f (log)')
ylabel('\mus')
pause
figure, semilogx(w/2/pi, 100*(tau0-tau)/tau0)
title('Greska grupnog kasnjenja')
xlabel('f (log)')
ylabel('%')

%% Cebisevljev PO
[n, wn] = cheb1ord(2*pi*[6600 8400], 2*pi*[5200 11500], 1, 26, 's');
[b, a] = cheby1(n, 1, wn, 's');
w = 2*pi*(1000:100:15000);
h = freqs(b, a, w);
figure, plot(w/2/pi, abs(h)); 
title('Amplitudna karakteristika filtra');
xlabel('f [Hz]')
ylabel('|H(j2\pif)|');
pause

%% Frekvencijske transformacije
% NP->VP
% Normalizovani Cebisevljev treceg reda
Wp = 1;
Ws = 2.5;
Ap = 3;
As = 20;

[n, wn] = cheb1ord(Wp, Ws, Ap, As, 's');
[bnp, anp] = cheby1(n, Ap, wn, 's');
[h, w] = freqs(bnp, anp);
figure, semilogx(w, abs(h)); 
title('Amplitudna karakteristika NP prototipa');
xlabel('\omega (log)')
ylabel('|H(j\omega)|');
pause
figure, semilogx(w, 20*log10(abs(h)));
title('Funkcija pojacanja NP prototipa')
xlabel('\omega (log)')
ylabel('20*log10(|H(j\omega)|) [dB]')
pause
figure, semilogx(w, angle(h));
title('Fazna karakteristika NP prototipa')
xlabel('\omega (log)')
ylabel('radijana')
pause

w = 2*pi*(0:10:1000e3);
[bvp, avp] = lp2hp(bnp, anp, 2*pi*2500);
h = freqs(bvp, avp, w);
figure, semilogx(w/2/pi, abs(h)); 
title('Amplitudna karakteristika filtra');
xlabel('f (log)')
ylabel('|H(j2\pif)|');
pause
figure, semilogx(w/2/pi, 20*log10(abs(h)));
axis([10 100e3 -100 20]);
title('Funkcija pojacanja filtra')
xlabel('f (log)')
ylabel('20*log10(|H(j2\pif)|) [dB]')
pause
figure, semilogx(w/2/pi, angle(h));
axis([10 100e3 -4 4]);
title('Fazna karakteristika filtra')
xlabel('f (log)')
ylabel('radijana')
pause
