clear;

fp = [5000 12000];
fs = [7000 9000];
N = input('Introduzca el orden del filtro =');
if rem(N,2)~=0,
	fprintf(1,'El orden de los filtros pasabanda debe ser par!!');
	N=N+1;
	fprintf(1,'Tomamos N = %d\n',N);
end
	

% Elegir una frecuencia de muestreo
fm = 60000;fny = fm/2;

% Frecuencias Digitales
Fp = fp/fm;Fs = fs/fm;

% Pasar a especificaciones de Pasobajo

F0 = (Fs(1)+Fs(2))/2;

% Hacer que las especificaciones sean simétricas
[minimum,idx]=min([Fs(1)-Fp(1) Fp(2)-Fs(2)]);
if idx==1,
	Fp(2)=Fs(2)+minimum;
else
      Fp(1)=Fs(1)-minimum;
end
FpLP=(Fs(2)-Fs(1))/2;FsLP=(Fp(2)-Fp(1))/2;

% Tomamos una frecuencia de corte intermedia (k<1)
k = 0.2;
Fc = FpLP+k*(FsLP-FpLP);

% Creamos el vector de indices
% El número de puntos es 41
% El número de puntos de un filtro Parabanda debe ser impar (secuencias tipo 1)
n=-N/2:N/2;

% Coeficientes del prototipo de filtro pasobajo ideal
hLPi = 2*Fc*sinc(2*n*Fc);

% Pasar al filtro Parabanda
delta = zeros(1,N+1);delta((N+2)/2)=1;
hbsi = delta-2*cos(2*pi*n*F0).*hLPi;

% Aplicar la ventana espectral
window = input('Ventana spectral? :','s');
window = window(1:4);
if window=='hamm',
	w = hamming(N+1)';
elseif window=='hann',
	w = hanning(N+1)';
elseif window=='blac',
	w = blackman(N+1)';
elseif window=='boxc',
	w = boxcar(N+1)';
end
B = hbsi.*w;
[H,f]=freqz(B,1,1500,fm);
[Hsp]=-20*log10(abs(freqz(B,1,[fp fs],fm)));
fprintf(1,'Atenuaciones en Pasabanda = %f\t%f\n',Hsp(1),Hsp(2));
fprintf(1,'Atenuaciones en Parabanda = %f\t%f\n',Hsp(3),Hsp(4));
plot(f,20*log10(abs(H)));zoom;grid;hold;
plot([fp fs],-Hsp,'*r');hold off;