clear;
close all;

% Especificaciones del filtro digital
wp = [3000 5000];Rp = 1; 
ws = [2000 7500];Rs = 60;
fm = 5000;

% 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 = wpx(2)-wpx(1);
% Fijamos las frecuencias de pasabanda
wx2 = wpx(1)*wpx(2);wx=sqrt(wx2);
if wsx(1)*wsx(2) < wx2,
	wsx(2) = wx2/wsx(1);
else
	wsx(1) = wx2/wsx(2);
end
Bw2 = wsx(2)-wsx(1);

% Criterio de MATLAB
% Bw1 = 1;
% Bw2 = (wsx.^2-wx2)/(wsx*(wpx(2)-wpx(1)));
% 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 Butterworth
e2 = k1;
n = ceil((log10(sqrt(k2/k1)))/log10(Ws/Wp));
v3 = (1/k1)^(1/(2*n));

[Z,P,K] = buttap(n);
[Bp,Ap] = zp2tf(Z,P,K);
[Bs,As] = lp2bp(Bp,Ap,wx,Wp*v3);
[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,'r');hold;plot([wp ws]/(2*pi),-[Rp Rp Rs Rs],'*k');grid;zoom;
figure;plot(F*fm,Gd);zoom;grid;

% Hacer el mismo filtro pero utilizando funciones de Matlab
% wpx = 2*tan(wp/(2*fm));
% wsx = 2*tan(ws/(2*fm));

% [N,Wn] = buttord(wpx,wsx,Rp,Rs,'s');
% [Bs,As] = butter(N,Wn,'s');
% [B,A] = bilinear(Bs,As,1);

% H = freqz(B,A,F*fm,fm);
% H = 20*log10(abs(H));

% plot(F*fm,H,'b');hold off;