%% Ejemplo de filtro pasobajo de Chebyshev-I analógico

clear all;

fp=20e3;wp=2*pi*fp;
Ap=0.5;
fs=90e3;ws=2*pi*fs;
As=50;

% Calculo de e2
e2 = 10^(0.1*Ap)-1;

% Orden del filtro
n = acosh(sqrt((10^(0.1*As)-1)/e2))/acosh(fs/fp);
n = ceil(n);

% Parámetros  para el cálculo de los polos y ceros del filtro 

theta = (1/n)*asinh(1/sqrt(e2));
k=1:n;
alpha = (2*k-1)*pi/(2*n);

% Determinacion de los polos en el prototipo
pk = -sinh(theta)*sin(alpha) +j*cosh(theta)*cos(alpha);

[B,A]=zp2tf([],pk,1);

% Ajustar ganancia para que |H(0)|=1
if rem(n,2)==0
  K=A(length(A))*(1/sqrt(1+e2));
else
  K=A(length(A));
end

[B1,A1] = lp2lp(K*B,A,wp);

N=500;
wf=2*pi*100e3;

W=0:wf/N:wf;
H=freqs(B1,A1,W);

subplot(1,2,1); plot(W/(2*pi),20*log10(abs(H)));grid;
ylabel(['|H(w)|']);
xlabel(['Frecuencia (Hz)']);
subplot(1,2,2); plot(W/(2*pi),20*log10(abs(H)));grid;
ylabel(['|H(w)|']);
xlabel(['Frecuencia (Hz)']);
axis([0 1*fp -Ap-0.1 0.1]);
print -depsc2 cheb1_ej
pause;close;

%% Calcular el retraso de grupo
dW=W(2)-W(1);
phase=unwrap(angle(H));
dphase = -diff(phase)/dW;

plot(W(1:length(W)-1)/(2*pi),dphase);grid;
ylabel(['Retraso de Grupo (s)']);
xlabel(['Frecuencia (Hz)']);
print -depsc2 cheb1_groupdelay_ej
pause;close;
plot(pk,'x');axis([-1.1 0 -1.1 1.1]);axis equal;