clear;

% Especificaciones del filtro digital
wp = 4000;Rp = 2; 
ws = 7000;Rs = 40;
fm = 22000;
%Transformación bilineal
% Predistorsionar las especificaciones
wpx = 2*tan(wp/(2*fm));
wsx = 2*tan(ws/(2*fm));

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

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] = lp2lp(Bp,Ap,Ws);


%Transformación bilineal
[B,A] = bilinear(Bs,As,1);
[H F]= freqz(B,A,1000,1);
H = 20*log10(abs(H));

plot(F*fm,H,'r');hold on;plot([wp ws]/(2*pi),-[Rp Rs],'*k');grid;zoom;
title('Transformación bilineal');
%Retardo de grupo
F = [1:1000]/3000;
Gd = grpdelay(B,A,F*fm,fm);
figure;plot(F*fm,Gd,'b');zoom;
title('retardo de grupo');