% Especificaciones en Hz
fp = [10000 20000]; fs = [2000 35000]; Ap = 1; As = 40;

% Frecuencia de muestreo
fm = 100e3;

% Convertir a rad/s
wp = 2*pi*fp;
ws = 2*pi*fs;

% Predistorsion
wp = 2*tan(wp/(2*fm));
ws = 2*tan(ws/(2*fm));

% Frecuencia Central
wx2 = wp(1)*wp(2); wx = sqrt(wx2);

% Mantengo las espec. de pasabanda y recalculo las de parabanda
if (ws(1)*ws(2) < wx2) ws(1) = wx2/ws(2); else ws(2) = wx2/ws(1); end

% Conversion a espec de prototipo filtro pasobajo
wp_p = wp(2)-wp(1);
ws_p = ws(2)-ws(1);

% Calculo e2
e2 = 1/(10^(0.1*As)-1) 

% Calculo del orden del filtro de Chebyshev II
n = acosh(sqrt(1/((10^(0.1*Ap)-1)*e2)))/acosh(ws_p/wp_p);

n = ceil(n);

% Prototipo de filtro pasobajo Chebyshev II
[Z,P,K] = cheb2ap(n,As)

[B,A] = zp2tf(Z,P,K)


% Conversion a pasabanda
% En Bw se pone ws_p porque Cheb2 normaliza con respecto a parabanda
[Bbp, Abp] = lp2bp(B,A,wx,ws_p);

% Vector de frecuencias para calcular resp. freq. (en rad/s)
W = 0:999; W = W*2*pi*fs(2)*1.5/1000;

H = freqs(Bbp, Abp, W);

% En db
Hdb = 20*log10(abs(H));

figure; plot(W/(2*pi), Hdb, [fp fs], -[Ap Ap As As],'r*'); hold

% Transformacion Bilineal
[Bd, Ad] = bilinear(Bbp, Abp, 1);

% Resp freq del filtro digital
[Hd, f] = freqz(Bd,Ad,1000, fm);

Hd_db = 20*log10(abs(Hd));

plot(f, Hd_db,'r--');