%% ===================================================================
%% FILTRO PASABAJOS CAJÓN
%%
%% Generación de un filtro pasabajos tomando muestras en frecuencia
%% entre [0,2pi]. El filtro obtenido es reordenado en tiempo para
%% hacerlo simétrico. Observar los problemas de ripple.
%% ===================================================================

clear all;

% --------------------------------------------------------------------
% Respuesta en frecuencia discretizada H(w_k)
% --------------------------------------------------------------------

N  = 31;      % longitud del filtro (numero impar para operador simétrico)
M  = (N-1)/2; % mitad de las muestras 
dw = 2*pi/N;  % intervalo de frecuencias
wc = pi/2;    % frecuencia de corte

H(1:N) = 0.;
for k = 1:N       % discretización de la respuesta en frecuencia ideal
  w = dw*(k-1);   % frecuencia angular digital
  if(w <= wc)
    H(k) = 1;
  elseif(w >= pi+wc)
    H(k) = 1;
  endif
endfor

% --------------------------------------------------------------------
% Operador en tiempo 
% --------------------------------------------------------------------

% Reacomodo coeficientes para que sea simetrico en tiempo

h = real(ifft(H));            % operador circular
h = [ h(M+1:-1:2) h(1:M+1)];  % operador simétrico

% --------------------------------------------------------------------
% GRAFICOS
% --------------------------------------------------------------------

figure(1); clf;

% Vector de la respuesta en frecuencia muestreada H_k = H(w_k)
ww = 0:2*pi/N:2*pi-2*pi/N; % [0,2*pi-dw]

% Ticks para los gráfico en frecuencia
xt      = 0:pi/2:2*pi;
xtlabel = ["0";"{/Symbol p}/2";"{/Symbol p}";"3/2{/Symbol p}";"{/Symbol 2p}"];
yt      = [0 1];

% Espectro del filtro
subplot(3,1,1)
stem(ww,abs(H),'r')
title('Filtro Pasabajos','fontsize',14)
ylabel('|H_k|')
axis([0 2*pi -0.2 1.4])
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt);

% Operador en tiempo
subplot(3,1,2)
stem((-M:M),h,'b')
axis([-M-1 M+1 4*min(h) 1.3*max(h)])
ylabel('h_k')

% Respuesta en frecuencia con zero padding 
HH = fft(h);

% Zero padding
hpad = [h,zeros(1,1000)];   % zero padding para observar la respuesta en
Hpad = fft(hpad); 	    % frecuencia del operador.
Npad = length(hpad);
wpad = 0:2*pi/Npad:2*pi-2*pi/Npad; % [0,2*pi-dw]

% Espectro con zero padding 
subplot(3,1,3)
stem(ww,abs(fft(h)),'r'); hold on
plot(wpad,abs(Hpad),'g'); hold off
ylabel('|H_k|')
axis([0 2*pi -0.2 1.4])
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt);

% --------------------------------------------------------------------
% Salida 
% --------------------------------------------------------------------

print('pasabajos.eps','-depsc');
