%% Diseño de un filtro fir pasabanda por el metodo de las series de
%Fourier

clear all; close all;

fs1 = 2e3;
fp1 = 10e3;
fp2 = 12e3;
fs2 = 16e3;
fm = 44e3;
%% Interesa que las especificaciones del filtro sean simétricas
%% por lo que calculo cual de las dos zonas de transición es menor
%% y modifico la frecuencia de parabanda

[mm,jj] = min([fp1-fs1 fs2-fp2]);
if jj==1
  fs2 = fp2+mm;
else
  fs1 = fp1-mm;
end

%% Calculo de los parámetros del filtro pasabanda Fc y F0

fc1 = (fp1+fs1)/2;
fc2 = (fp2+fs2)/2;

Fc1 = fc1/fm;
Fc2 = fc2/fm;

F0 = (Fc1+Fc2)/2;
Fc = (Fc2-Fc1)/2;

windows=['boxcar  ';'bartlett';'hamming ';'hanning ';'blackman'];
colors=['b  ';'r  ';'m  ';'--b';'--r'];

%% Numero de puntos del filtro
%% N puede ser par o impar ya que se puede construir un filtro pasobajo
%% siempre que h[n] sea par 
N = 45;
n=0:N-1;

%% Filtro pasobajo ideal
hlpi = 2*Fc*sinc(2*Fc*(n-(N-1)/2));

%% Aplicamos ventanas espectrales a este filtro ideal

M=500; %% Nº de puntos para calcular la Resp Freq
for i=1:size(windows)
  eval(['w=' windows(i,:) '(N);']);
  hlp = hlpi.*w';
  %% Normalizar de forma que |H(0)|=1
  hlp = hlp/sum(hlp);
  
  hbp = 2*hlp.*cos(2*pi*F0*(n-(N-1)/2));
  
  [H,f]=freqz(hbp,1,M,fm);
  
  plot(f,abs(H),colors(i,:));
  hold on;
  if (i==3) %% ventana de hamming
    figure;stem(n,hbp,'o');grid;pause;
    eval(['print -depsc2 fir_bp_coef_N' num2str(N)]);
    close;
    figure;plot(f,20*log10(abs(H)));grid;
    axis([-inf inf -60 1]);
    xlabel(['Frecuencia (Hz)']);
    ylabel(['|H(F|']);pause;
    eval(['print -depsc2 fir_bp_rf_N' num2str(N)]);
    close;
  end
end
grid;hold off;
legend('rect','bartlett','hamming','hanning','blackman');
xlabel(['Frecuencia (Hz)']);
ylabel(['|H(F|']);
title(['Filtro Pasobanda, fc=' num2str(Fc*fm/1000) 'KHz, f0=' ...
       num2str(F0*fm/1000) 'KHz, fm=' num2str(fm/1000) 'KHz,' ...
		    ' N=' num2str(N)]);

eval(['print -depsc2 fir_bp_wind_N' num2str(N)]);

pause;close