%% =================================================================
%% LOWPASS C
%% 
%% Filtro pasabajos con frecuencia de corte en pi/2
%% Diseño usando solo la mitad de las muestras 
%% aprovechando la paridad del operador en tiempo y en frecuencia, 
%% ya que busco que sea de fase cero

%% OBS!: En esta versión la TDF inversa se calcula meiante un
%%       sumatoria sin recurrir a la función intrínseca ifft.
%% =================================================================

clear all ;

% Parámetros del filtro

N  = 21 ;
M  = (N-1)/2 ;
W0 = 2*pi/N ;

% Defino las muestras en frecuenca hasta M=(N-1)/2 pues
% las otras M seran iguales ya que H es par:

for k = 0:M
    Wk = W0*k ;
    if Wk < pi/2
       H(k+1) = 1 ;
    else 
       H(k+1) = 0 ;
    end
end

% Calculo el operador lineal utilizando la fórmula que tiene en 
% cuenta que H es par (sólo hace falta calcular M+1 muestras):

for n = 0:M
    h(n+1) = H(1) ;
    for k = 1:M
	h(n+1) = h(n+1)+2*H(k+1)*cos(W0*k*n) ;
    end
    h(n+1) = h(n+1)/N ;
end

% Defino las muestras para los tiempos negativos.
% Como h es par uso los valores de los índice positivos
% que obtuve anteriormente:

h = [h(M+1:-1:2) h]; 

% Realizo el zero padding para tener la respuesta 
% en frecuencia bien muestreada:

h2     = h ;
NN     = 1024 ;
h2(NN) = 0 ;

% Calculo la espectro de amplitud del filtro obtenido:

H2 = abs(fft(h2)) ;

% Genero la respuesta teórica del filtro para comparar: 

W02 = 2*pi/NN ;
for k = 1:NN
    Wk = W02*(k-1) ;
    if Wk < pi/2
       Ht(k) = 1 ;
    elseif Wk > 3*pi/2
       Ht(k) = 1 ;
    else 
       Ht(k) = 0 ;
    end
end

% Grafico el filtro en tiempo:

subplot(2,1,1);
stem((-M:M),h,'b'); hold on
plot((-M:M),h,';filtro en tiempo;o'); hold off
xlabel('n')
ylabel('h_n')
axis([-M-1 M+1 -0.3 0.7]);

% Grafico filtro en frecuencia para ver qué
% tanto me aparto del pasa-bajos teórico:

subplot(2,1,2);
plot(W0*(0:M),H,';muestras originales;bo',\ 
     W02*(0:(NN-1)),H2,';espectro real;r-',\
     W02*(0:(NN-1)),Ht,';espectro teorico;g-') ;
legend("location","north")
xlabel('{/Symbol w}')
ylabel('|H{/Symbol w}|')
xt      = 0:pi/2:2*pi;
xtlabel = ["0";"{/Symbol p}/2";"{/Symbol p}";"3/2{/Symbol p}";"{/Symbol 2p}"];
set(gca,"xtick",xt,"xticklabel",xtlabel);
set(gca,"ytick",[0 1]);
axis([0,2*pi,-0.2,1.2]);

% Guardo un .eps con el gráfico:

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