% Problema 1
% fp=2 KHz Ap=0.5 dB
% fs=10 KHz As=50 dB

fp=2000;Ap=0.5;
fs=10000;As=50;
At=-[Ap As];
fm=44100;

% Filtro de Butterworth

% Calculo de e2, N y v3
e2=10^(0.1*Ap)-1;
n=(log10(sqrt((10^(0.1*As)-1)/e2)))/log10(fs/fp);
N=ceil(n);
v3=(1/e2)^(1/(2*N));

% Prototipo de Filtro pasobajo

disp('******** Prototipo de Filtro Pasobajo **********');
[Z,P,K]=buttap(N);
[B,A]=zp2tf(Z,P,K);
[H,W]=freqs(B,A,500);
plot(W,20*log10(abs(H)));hold on;zoom;grid;
We=[1 v3 fs/fp]/v3;
plot(We,-[Ap 3 As],'*r');hold off;
title('Prototipo de Filtro Pasobajo');
disp(['N = ' int2str(N)]);
disp(['e2 = ' num2str(e2)]);
disp(['v3 = ' num2str(v3)]);
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,v3*2*pi*fp);
[H,W]=freqs(B,A,500);
plot(W/(2*pi),20*log10(abs(H)));hold on;zoom;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]=buttord(We(1),We(2),Ap,As,'s');
[Bs,As]=butter(N,Wn,'s');
[Bz,Az]=bilinear(Bs,As,1);
[Hz,Fz]=freqz(Bz,Az,500,fm);
plot(Fz,20*log10(abs(Hz)));hold on;zoom;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*BzII;
[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;zoom;
title('Retraso de Fase - Filtro Digital, (-r) T. Bilineal, (:g) T. Inv. Impulso');

% Respuesta a impulso
delta=zeros(500,1);delta(1)=1;
hbl=filter(Bz,Az,delta);
hii=filter(BzII,AzII,delta);
t=0:499;
disp('Pulse una tecla para seguir');pause;
plot(t,hbl,'r',t,hii,'g');zoom;grid;
title('Respuesta a Impulso del filtro digital, (-r) T. Bilineal, (g) T. Inv. Impulso');