% Problema 1
% fp=2 KHz Ap=0.5 dB
% fs=10 KHz As=50 dB
clear;
fp=2000;Ap=0.5;
fs=10000;As=50;
At=-[Ap As];
fm=44100;

% Filtro de Chebychev I
% Calculo de e2 y N
e2 = 10^(0.1*Ap)-1;
N = ceil(acosh(sqrt((10^(0.1*As))/e2)-1)/acosh(fs/fp));

% Prototipo de Filtro pasobajo

disp('******** Prototipo de Filtro Pasobajo **********');
[Z,P,K]=cheb1ap(N,Ap);
[B,A]=zp2tf(Z,P,K);	
[H,W]=freqs(B,A,500);
plot(W,20*log10(abs(H)));hold on;zoom;grid;
We=[1 fs/fp];
plot(We,-[Ap As],'*r');hold off;
title('Prototipo de Filtro Pasobajo');
disp(['N = ' int2str(N)]);
fprintf(1,'B = \n');
fprintf(1,'    %.4f\n',B);
fprintf(1,'A = \n');
fprintf(1,'    %.4f\n',A);
disp('***************************************************');
disp('Pulse una tecla para seguir');pause;

% Pasar el Filtro pasobajo

disp('********* Filtro pasobajo ************');
[B,A]=lp2lp(B,A,2*pi*fp);
[H,W]=freqs(B,A,500);
plot(W/(2*pi),20*log10(abs(H)));hold on;grid;
We=[fp fs]*2*pi;
plot(We/(2*pi),At,'*r');hold off;
title('Filtro Pasobajo');
fprintf(1,'B = \n');
fprintf(1,'    %.4f\n',B);
fprintf(1,'A = \n');
fprintf(1,'    %.4f\n',A);
disp('Pulse una tecla para seguir');pause;

% Filtro Digital a partir de la Transformación Bilineal
% w=2tan(2*pi*F/2);

Fp=fp/fm;Fs=fs/fm;

% Predistorsionado
disp('********** Filtro Digital-Transformación Bilineal ********');
Fe=[Fp Fs];fe=Fe*fm;
We=2*tan(2*pi*Fe/2);
[N,Wn]=cheb1ord(We(1),We(2),Ap,As,'s');
[Bs,As]=cheby1(N,Ap,Wn,'s');
[Bz,Az]=bilinear(Bs,As,1);
[Hz,Fz]=freqz(Bz,Az,500,fm);
plot(Fz,20*log10(abs(Hz)));hold on;grid;
plot(fe,At,'*r');
title('Filtro Digital, (-y) T. Bilineal, (:g) T. Inv. Imp.');
fprintf(1,'B = \n');
fprintf(1,'    %.4f\n',Bz);
fprintf(1,'A = \n');
fprintf(1,'    %.4f\n',Az);
[Gz,Fz]=grpdelay(Bz,Az,500,fm);

disp('Pulse una tecla para seguir');pause;

disp('********* Filtro Digital-Transformación Invariante al Impulso *********');
% Los coeficientes del filtro pasobajo eran B y A
[BzII,AzII]=impinvar(B,A,fm);
% Hay que normalizar para que H(0)=1
k=sum(real(AzII))/sum(real(BzII));
BzII=k*real(BzII);AzII=real(AzII);
[HzII,FzII]=freqz(BzII,AzII,500,fm);
plot(FzII,20*log10(abs(HzII)),':g');hold off;
fprintf(1,'B = \n');
fprintf(1,'    %.4f\n',BzII);
fprintf(1,'A = \n');
fprintf(1,'    %.4f\n',AzII);
disp('Pulse una tecla para seguir');pause;
[GzII,FzII]=grpdelay(BzII,AzII,500,fm);
plot(Fz,Gz,'-r',FzII,GzII,':g');grid;
title('Retraso de Fase - Filtro Digital, (-r) T. Bilineal, (:g) T. Inv. Impulso');

disp('Pulse una tecla para seguir');pause;

d=zeros(1,200);
n=[0:199];
d(1)=1;
yB=filter(Bz,Az,d);
yII=filter(BzII,AzII,d);
plot(n,yB,'-r',n,yII,':g');
title('Respuesta a Impulso, (-r) T. Bilineal, (:g) T. Inv. Imp.');zoom;grid;