% ==================================================================
% FILTRO PASO BAJO TOMANDO MUESTRAS EN FRECUENCIA
%
% h:  filtro sin cambio de muestra
% hs: filtro con cambio de muestra
% H:  respuesta en frecuencia SIN cambio de muestra
% HS: respuesta en frecuencia CON cambio de muestra
%
% !Obs: los filtros están definidos solo usando las M+1 muestras
%       que son diferentes. No se usa FFT.   
% ==================================================================

clear all; close all;

% -------------------------------------------------------------------
% Parametros del filtro
% -------------------------------------------------------------------

N  =  21 ;     % número impar (fase cero)
M  = (N-1)/2 ; % Nyquist (pi)
wc =  pi/2 ;   % frecuencia de corte

% Valor que toma la muestra cercana a la frec. de corte:

sample = 0.5;

% -------------------------------------------------------------------
% Defino el filtro (sólo M+1 muestras)
% -------------------------------------------------------------------

dw = 2*pi/N ; 
for k = 0:M   
  w(k+1) = dw * k ;
  if w(k+1) < wc
    H(k+1) = 1 ;
  else
    H(k+1) = 0 ;
  end
end

% Cambio la muestra en frecuencia cercana a la frecuencia de corte

Nwc     = ceil(wc/dw); % índice la de muestra que cambio
HS      = H;
HS(Nwc) = sample;      % cambio la muestra

% -------------------------------------------------------------------
% Filtros en tiempo (sólo M+1 muestras)
% -------------------------------------------------------------------

for n = 0:M
  h(n+1)  = H(1);
  hs(n+1) = HS(1);
  for k = 1:M
    h(n+1)  = h(n+1)  + 2 * H(k+1) * cos(dw*k*n);
    hs(n+1) = hs(n+1) + 2 * HS(k+1)* cos(dw*k*n);
  end
  h(n+1) =  h(n+1)/N;
  hs(n+1) = hs(n+1)/N;
end
     
% -------------------------------------------------------------------
% Zero padding (sin TDF)
% -------------------------------------------------------------------

N2 = 1001;      % como no uso FFT no hace falta que sea potencia de 2
M2 = (N2-1)/2;

% Defino los filtros en frecuencia con la rama cosenoidal de la TDF
% a indices positivos.

dw2 = 2*pi/N2 ;
for k = 0:M2
  H2(k+1)  = h(1);
  HS2(k+1) = hs(1);
  for n = 1:M % de M+1 a M2 es 0
    H2(k+1)  =  H2(k+1) + 2 * h(n+1)  * cos(dw2*k*n);
    HS2(k+1) = HS2(k+1) + 2 * hs(n+1) * cos(dw2*k*n);
  end
end
      
% -------------------------------------------------------------------
% Filtros de fase 0
% -------------------------------------------------------------------

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

% ==================================================================
% GRAFICOS
% ==================================================================

figure(1)

% Operadores en tiempo (h y hs)

subplot(3,1,1)
hold on
stem((-M:M),h,';h_n;b-o','markersize',3,'filled');
stem((-M:M),hs,';hs_n;r-o','markersize',4,'filled');
hold off
axis([-M-1,M+1,-0.3,0.7])
ylabel('h_n')

xt      =  0:pi/2:pi;
xtlabel = ["0";"{/Symbol p}/2";"{/Symbol p}"];
yt      = [0 1];

% Respuestas en frecuencia (H y HS)

subplot(3,1,2)
hold on
plot((0:M2).*dw2,H2, ';H;b-',  w,H, 'bo','markersize',3);
plot((0:M2).*dw2,HS2,';HS;r-', w,HS,'ro') ;
hold off
ylabel('H(\omega)')
axis([0,pi,-0.2,1.2]);
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt);

% Respuesta en frecuencia en decibeles

subplot(3,1,3)
hold on
plot((0:M2).*dw2,20.*log10(abs(H2)),';20 log(H);b-');
plot((0:M2).*dw2,20.*log10(abs(HS2)),';20 log(HS);r-');
hold off
ylabel('H(\omega)[dB]')
axis([0,pi,-100,10]);
set(gca,"xtick",xt,"xticklabel",xtlabel);

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

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