%program13- inverzní Èebyšev postup dle uèebnice Pøíklad 4.5
clear
format long
ap=1
as=38
Os=2.35
epsilon=sqrt(10^(ap/10)-1)
k1=sqrt(10^(ap/10)-1)/sqrt(10^(as/10)-1)
n=ceil(acosh(1/k1)/acosh(Os))
k=1/cosh(acosh(1/k1)/n)
B=(sqrt(1+(epsilon/k1)^2)+(epsilon/k1))^(1/n)
a=(B-1/B)/2
b=(B+1/B)/2
%format long
mu=(1:2:2*n-1)'
p=1./(-a*sin(mu*pi/(2*n))+i*b*cos(mu*pi/(2*n)))/k
Hjme=real(poly(p))             %jmenovatel
if rem(n,2)                    %n - liche
    mu=mu([1:end]~=(n+1)/2)
    z=i./cos(mu*pi/(2*n))/k    %nuly pren. funkce
    Hcit=real(poly(z))*n*k1/epsilon/k %citatel pren. funkce
    Pcit=eye(1,n)*k/n/k1       %citatel char. funkce
else                           %n - sude
    z=i./cos(mu*pi/(2*n))/k    %nuly pren. funkce
    Hcit=real(poly(z))/sqrt(epsilon^2/k1^2+1) %citatel pren. funkce
    Pcit=eye(1,n)/k1           %citatel char. funkce
end
Pjme=real(poly(z))             %jmenovatel char. funkce

A=Hjme
B1=Hcit
wk=2*Os
w=linspace(0,wk,200);

H=freqs(B1,A,w);
figure(1)
plot(w,abs(H));
grid
xlabel w
ylabel '|H|'
title 'Modulová charakteristika'
figure(2)
plot(w,phase(H));
grid
xlabel w
ylabel 'PH'
title 'Fázová charakteristika'
t=diff(-phase(H))./diff(w);
figure(3)
plot(w(1:199),t);
t1=grdelays(B,A,w);
figure(4)
plot(w,t1)
grid
xlabel w
ylabel 't'
title 'Skupinové zpoždìní'
