clear all;
close all;

% Especificaciones del filtro digital
wp = [1000 3000];Rp = 3; 
ws = [2000 2500];Rs = 40;
fm = 10000;

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

% Paso a especificaciones de pasobajo
% Tablas 1 y 2
Bw1 = wsx(2)-wsx(1);
% Fijamos las frecuencias de parabanda
wx2 = wsx(1)*wsx(2);wx=sqrt(wx2);
if wpx(1)*wpx(2) < wx2,
	wpx(1) = wx2/wpx(2);
else
	wpx(2) = wx2/wpx(1);
end
Bw2 = wpx(2)-wpx(1);

% Criterio de MATLAB
% Bw1 = 1;
% Bw2 = (wsx*(wpx(2)-wpx(1)))./(wsx.^2-wx2);
% Bw2 = min(abs(Bw2)); 


% Especificiones de pasobajo
Wp = Bw1;
Ws = Bw2;

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

% Filtro de Cheby1
e2 = k1;
n = ceil((acosh(sqrt(k2/k1)))/acosh(Ws/Wp));
v3 = (1/k1)^(1/(2*n));

[Z,P,K] = cheb1ap(n,Rp);
[Bp,Ap] = zp2tf(Z,P,K);
[Bs,As] = lp2bs(Bp,Ap,wx,Ws);
[B,A] = bilinear(Bs,As,1);

F = [0:499];F= F*2*wp(2)/(2*pi*500);
H = freqz(B,A,F,fm);
H = 20*log10(abs(H));

figure;plot(F,H,'r',[wp ws]/(2*pi),-[Rp Rp Rs Rs],'*k');grid;zoom;

% Hacer el mismo filtro pero utilizando funciones de Matlab


[N,Wn] = cheb1ord(2*wp/(2*pi*fm),2*ws/(2*pi*fm),Rp,Rs);
[B,A] = cheby1(N,Rp,Wn,'stop');

H = freqz(B,A,F,fm);
Hdb = 20*log10(abs(H));
hold;
plot(F,Hdb,'b',[wp ws]/(2*pi),-[Rp Rp Rs Rs],'*k');