%% ===================================================================
%% FILTRO PASABAJOS CAJÓN USANDO LA TRANSFORMADA DE HILBERT
%%
%% Generación de un filtro pasabajos a partir de la respuesta en 
%% frecuencia de su función analítica GH(w). Si H(w) es la respuesta
%% en frecuencia del pasabajon cajón, entonces su función analítica
%% en frecuencia está dada por:
%%          _
%%         |    0,     w < 0  (w > w_Nyquist)
%% GH(w) = |    H(0),  w = 0
%%         |_ 2*H(w),  w > 0  (0 < w < w_Nyquist)
%%
%% El filtro deseado se obtiene antitransformando al dominio del 
%% tiempo,  reordenandodo el operador para que sea simétrico, y como  
%% 
%% gh = h - i*Hilb{h}
%% 
%% se toma finalmente la parte real para obtener el filtro en tiempo.
%% Con esta forma de generar el filtro, no es necesario pedir la paridad
%% del espectro respecto de Nyquist.
%% ===================================================================

clear all; close;

% ---------------------
% Parámetros del filtro
% ---------------------

N  = 21;      % longitud del filtro (número 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 

% ---------------------------------------------------------------------
% Defino la respuesta en frecuencia angular digital del pasabajos cajón
% ---------------------------------------------------------------------

for k = 1:N    
  w = dw*(k-1);
  if(w <= wc)
    H(k) = 1;
  elseif w >= pi+wc
    H(k) = 1;
  endif
endfor

% -------------------------------------------------------------
% Respuesta en frecuencia de la función analítica deel filtro H
% -------------------------------------------------------------

GH(1) = 1  ;        % w = 0. Conservo la energía del filtro en tiempo.
for k = 2:N     
  w = dw*(k-1); 
  if(w <= wc)
    GH(k) = 2*H(k); % duplico los valores para 0 < w < w_Nyquist
  else
    GH(k) = 0;      % cancelo los valores para w > w_Nyquist (w < 0)
  endif
endfor

% ---------------------------
% Función analítica en tiempo 
% ---------------------------

gh =   ifft(GH);                     % operador circular complejo
gh = [ gh(M+2:N) conj(gh(1:M+1)) ];  % operador simétrico complejo

% ----------------
% Filtro en tiempo
% ----------------

ha = real(gh);  % el filtro deseado está en la parte real 

% --------
% GRÁFICOS
% --------

figure(1); clf;

xt      =   0 : pi/2 : 2*pi;
xtlabel = ["0";"{/Symbol p}/2";"{/Symbol p}";"3/2{/Symbol p}";"{/Symbol 2p}"];
yt      =  -0.5 : 0.5 : 2;

% Espectro del pasabajos cajón H_k = H(w_k)

subplot(4,1,1)

w = 0 : dw : 2*pi-dw;  % [0,2*pi-dw]

stem(w,abs(H),'r')
axis([-0.1 2*pi+0.1 -0.4 2.4])
set(gca,"xtick",xt,"xticklabel",xtlabel);
ylabel('|H_k|')
title('Filtro pasabajos','fontsize',14)

% Espectro de la transformada de Hilbert GH_k = GH(w_k)

subplot(4,1,2)

stem(w,abs(GH),'r')
axis([-0.1 2*pi+0.1 -0.4 2.4])
set(gca,"xtick",xt,"xticklabel",xtlabel);
ylabel('|GH_k|')

% Operador en tiempo

subplot(4,1,3)

stem((-M:M),ha,'b')
axis([ -M-1 M+1 min(ha)-0.2 max(ha)+0.2 ])
ylabel('ha_n')

% Espectro con zero padding

subplot(4,1,4)

HA   = fft(ha);
Npad = 1024;
hpad = [ha zeros(1,1024-N)]; % zero padding para observar la respuesta en
Hpad = fft(hpad); 	     % frecuencia del operador.

dwpad = 2*pi/Npad; 
wpad  = 0 : dwpad : 2*pi-dwpad; % [0,2*pi-dwpad]

stem(w,abs(HA),'r'); hold on
plot(wpad,abs(Hpad),'g'); hold off
axis([ -0.1 2*pi+0.1 -0.4 1.4 ])
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt);
ylabel('|HA_k|')

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

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