clear;

% Especificaciones del filtro digital
wp = 3000;Rp = 1; 
ws = 2000;Rs = 60;
fm = 4000;

% Predistorsionar las especificaciones
wpx = 2*tan(wp/(2*fm));
wsx = 2*tan(ws/(2*fm));

% Se hace paso paso el filtro analógico
Wp = wsx;
Ws = wpx;

k1 = 10^(0.1*Rp)-1;
k2 = 10^(0.1*Rs)-1;

e2 = 1/k2;
n = ceil((acosh(sqrt(k2/k1)))/(acosh(Ws/Wp)));

[Z,P,K] = cheb2ap(n,Rs);
[Bp,Ap] = zp2tf(Z,P,K);
[Bs,As] = lp2hp(Bp,Ap,Wp);

[B,A] = bilinear(Bs,As,1);
F = [1:1000]/2000;
H = freqz(B,A,F*fm,fm);
H = 20*log10(abs(H));
Gd = grpdelay(B,A,F*fm,fm);
plot(F*fm,H,'b');hold;plot([wp ws]/(2*pi),-[Rp Rs],'*k');grid;zoom;
figure;plot(F*fm,Gd,'r');