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

clear all; close all;

fp1 = 2e3;
fs1 = 10e3;
fs2 = 12e3;
fp2 = 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([fs1-fp1 fp2-fs2]);
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 solo puede ser impar y la secuencia deber ser par ya que es la
%unica forma de diseñar un filtro FIR parabanda
N = 41;
n=0:N-1;

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

%% Calculo de delta
delta = zeros(1,N);
delta((N+1)/2) = 1;

%% 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);
  
  hbs = delta - 2*hlp.*cos(2*pi*F0*(n-(N-1)/2));
  
  [H,f]=freqz(hbs,1,M,fm);
  
  plot(f,abs(H),colors(i,:));
  hold on;
  if (i==5) %% ventana de blackman
    figure;stem(n,hbs,'o');grid;pause;
    eval(['print -depsc2 fir_bs_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_bs_rf_N' num2str(N)]);
    close;
  end
  
end
grid;hold off;
legend('rect','bartlett','hamming','hanning','blackman');
xlabel(['Frecuencia (Hz)']);
ylabel(['|H(F|']);
title(['Filtro Parabanda, fc=' num2str(Fc*fm/1000) 'KHz, f0=' ...
       num2str(F0*fm/1000) 'KHz, fm=' num2str(fm/1000) 'KHz,' ...
		    ' N=' num2str(N)]);

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

pause;close