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

clear all;

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

% Cálculo de e2
e2 = 1/(10^(0.1*As)-1);

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

% Parámetros para el cálculo de polo y ceros 
theta = (1/n)*asinh(1/sqrt(e2));
k=1:n;
alpha = (2*k-1)*pi/(2*n);

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

k = 1:n/2;
phi = (2*k-1)*pi/(2*n);
zk = [j./cos(phi) -j./cos(phi)];

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

K=A(length(A))/B(length(B));

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

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.1*fp -Ap-0.1 0.1]);
print -depsc2 cheb2_ej
pause;
close;

% Retraso de grupo
dW=W(2)-W(1);
phase=unwrap(angle(H));
grpdelay = -diff(phase)/dW;
%% Eliminar valores negativos (discontinuadades)
idx = find(grpdelay<0);
grpdelay(idx) = grpdelay(idx-1);

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