%% Ejemplo filtro pasabanda de butterworth discreto
clear all;

% Especif.
fp1=45e3;fp2=90e3;
Ap=0.5;
fs1=20e3;fs2=150e3;
As=50;
fm = 500e3;

wp1 = 2*pi*fp1;
wp2 = 2*pi*fp2;
ws1 = 2*pi*fs1;
ws2 = 2*pi*fs2;

%% Predistorsionar
%% Obtenemos especs de filtro analógico
wap1 = 2*tan(wp1/2/fm);
wap2 = 2*tan(wp2/2/fm);
was1 = 2*tan(ws1/2/fm);
was2 = 2*tan(ws2/2/fm);

%% Fijamos frecuencias de pasabanda
wax2 = wap1*wap2;

if (was1*was2 < wax2) 
  was1 = wax2/was2;
else
  was2 = wax2/was1;
end

wlpp = wap2-wap1;
wlps = was2-was1;


%% Prototipo Pasobajo

% Calculo de e2
e2 = 10^(0.1*Ap)-1;

% Orden del filtro
n = log(sqrt((10^(0.1*As)-1)/e2))/log(wlps/wlpp);
n = ceil(n);

% Frecuencia de corte
wac = wlpp/(e2^(1/(2*n)));

k=1:n;
alpha = (2*k-1)*pi/(2*n);

% Determinacion de los polos en el prototipo
pk = -sin(alpha) +j*cos(alpha);

[B,A]=zp2tf([],pk,1);

[B1,A1] = lp2bp(B,A,sqrt(wax2),wac); 

[Bz,Az] = bilinear(B1,A1,1);

N=500;
ff=160e3;

f=(5:N-1)*ff/N;
Hz=freqz(Bz,Az,f,fm);

subplot(3,1,1); plot(f,20*log10(abs(Hz)));grid;
ylabel(['|H(w)|']);
xlabel(['Frecuencia (Hz)']);
subplot(3,1,2); plot(f,unwrap(angle(Hz)));grid;
ylabel(['Fase H(w) (rad)']);
xlabel(['Frecuencia (Hz)']);

% Retraso de grupo
gd = grpdelay(Bz,Az,f,fm);
subplot(3,1,3);plot(f,gd);grid;
ylabel(['Retraso de Grupo (muestreos)']);
xlabel(['Frecuencia (Hz)']);
print -depsc2 bp_butter_discr_ej