%% =================================================================
%% FILTRADO DE UNA SINUSOIDE EN EL DOMINIO DE LAS FRECUENCIAS
%%
%% Se quiere filtrar una sinusoide en el dominio de las frecuencias.
%% La frecuencia de la sinusoide es tal que cae entre dos muestras
%% de la TFD del filtro. 
%% Si bien la frecuencia de la se�al es mayor a la frecuencia de 
%% corte del filtro, la representaci�n de la se�al en el dominio de 
%% la TDF no es precisa y contiene elementos no nulos en la banda de 
%% paso. Por lo tanto, no se logra filtrar la se�al correctamente.
%% =================================================================

%clear all;

N  = 31 ;
dw = 2*pi/N ;
j  = 4 ;
wc = dw*j ;
ws = dw*(j+3/2) ;

% Genero una sinusoide de frecuencia ws y calculo su espectro via TDF

n = 0:N-1 ;
a = sin(ws*n) ;
A = abs(fft(a)) ;

% Defino las N muestras en frecuencia del pasabajos

for k = 0:N-1
  wk = dw*k ;
  if wk <= wc
    H(k+1) = 1;
  elseif wk >= 2*pi - wc
    H(k+1) = 1;
  else 
    H(k+1) = 0 ;
  endif
end

% Filtro en frecuencias multiplicando muestra a muestra

AF = A.*H ;

% Vuelvo al tiempo y obtengo la se�al filtrada

af = real(ifft(AF));


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

% Ahora vuelvo a calcular su espectro con zero padding

a2     = a;
NN     = 1024;
a2(NN) = 0;
A2     = abs(fft(a2));

% Defino las NN muestras en frecuencia del pasabajos te�rico
% solo para graficar

dw2 = 2*pi/NN ;
for k = 0:NN-1
  wk = dw2*k ;
  if wk <= wc
    Ht(k+1) = 1;
  elseif wk >= 2*pi - wc
    Ht(k+1) = 1;
  else 
    Ht(k+1) = 0 ;
  endif
end

% Ticks para los gr�ficos

xt      = [0 ws pi 2*pi-ws 2*pi];
xtlabel = ["0" ; "{/Symbol w}s" ; "{/Symbol p}" ; 
           "2{/Symbol p}-{/Symbol w}s" ; "2{/Symbol p}"];
yt      = [0 1];

% --------------------------------------------------------------------
% Grafico la ond�cula y su espectro con y sin zero padding
% --------------------------------------------------------------------

figure(1); clf;

subplot(2,1,1)
stem(n,a)
title('Ondicula')
axis([-1,N,-1.2,1.2]);

subplot(2,1,2)
stem(n*dw,A/max(A2),'b-x'); hold on;
stem([ws 2*pi-ws],[1 1],'r-xs');
plot((0:NN-1)*dw2,A2/max(A2),'g-'); hold off
title('Espectro calculado via TDF')
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt);
axis([0,2*pi,-0.1,1.1]);

% Guardo el gr�fico en un postscript

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

% --------------------------------------------------------------------
% Grafico los espectros de la ond�cula, del filtro, 
% y de la ondicula filtrada
% --------------------------------------------------------------------

figure(2); clf;

subplot(2,2,1);
stem(n*dw,zeros(1,N),'b-x'); hold on;
 hold off;
title('Espectro teorico de la sinusoide')
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt);
axis([0,2*pi,-0.1,1.1]);

subplot(2,2,2);
stem(n*dw,A/max(A2),'b-x');
title('Espectro calculado via TDF')
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt);
axis([0,2*pi,-0.1,1.1]);

subplot(2,2,3);
plot((0:NN-1)*dw2,Ht,'g-'); hold on
plot(n*dw,H,'go'); hold off
title('Espectro del pasabajos')
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt);
axis([0,2*pi,-0.1,1.1]);

subplot(2,2,4);
stem(n*dw,AF/max(A2),'c-x');
title('Espectro de la ondicula filtrada')
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt);
axis([0,2*pi,-0.1,1.1]);

% Guardo el gr�fico en un postscript

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

% --------------------------------------------------------------------
% Se�al de entrada y se�al filtrada
% --------------------------------------------------------------------

figure(3); clf;

plot(n,a,'b-o;a_n;'); hold on
plot(n,af,'r-o;af_n;')
title('Senal de entrada y senal filtrada')
axis([-1,N,-1.2,1.2]);

% Guardo el gr�fico en un postscript

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

