% =================================================================
% Definicion y aplicacion de un filtro de un unico polo sobre
% una senal simple.
%
% Obs:  se toma  un  numero  de coeficientes  bajo  para resaltar  las
% diferencias  entre  la  respuesta  teorica  y la  del  filtro  de  N
% coeficientes.
% =================================================================

% =================================================================
clear all;close all;


% =================================================================
% Parametros de interes del filtro
% =================================================================

% filtro de un  solo polo

N   = 11;    % numero total de coeficientes del filtro (impar)
rho = 1.2;  % fuera de |z|=1 para que sea de fase minima
wp  = pi/4; % frecuencia de paso



% =================================================================
% Definicion del filtro
% =================================================================


zp  = rho*exp(-i*wp);  % ubicacion del polo
dw  = 2*pi/N;   
w   = (0:N-1)*dw; % frecuencia angular digital

% Espectro de amplitud de la respuesta en frecuencia muestreada

M   = (N-1)/2; % coeficientes a un lado de la muestra central

% i) usando el espectro de amplitud directamente:

for j=1:M+1
    H(j) = sqrt(1/(1+rho^2-2*rho*cos(w(j)-wp))); % hasta Nyquist
end

% ii)  empleando la transformada Z:
%z  = exp(-i*w);       % en |z|=1
%H  = 1 ./ (z(1:M+1)-zp); % hasta Nyquist 
%H = abs(H); % espectro de amplitud

norm = 1/max(H); % para normalizar la ganancia

H    = [H(1:M+1) H(M+1:-1:2) ]; % hermitico (par respecto de Nyquist para |H|)
H    = norm * H;  % normalizo la ganancia

% Operador en tiempo

% a) sin cambio de fase
h = real(ifft(H));        % operador circular en tiempo
h = [h(M+2:N) h(1:M+1)];  % operador lineal en tiempo

% b) con cambio de fase
% H     = H.*exp(-1i.*w*M);  % para obtener operador lineal en tiempo
% h     = real(ifft(H));     % operador en tiempo


% =================================================================
% Definicion de la senal
% =================================================================

w0=wp;      % frecuencia angular digital de la componente a destacar
w1=pi/2;    % frecuencia angular digital de la componente a atenuar
Td=2*pi/w0; % periodo digital para  tomar muestra representativa de la
            % senal
n=(0:10*Td);
ns=length(n);
s1=0.5*cos(w0*n);
s2=0.5*cos(w1*n);
s=s1+s2;


% =================================================================
% Filtrado de la senal en tiempo 
% =================================================================

sf=conv(s,h)( M+1: M+1 + ns-1); % =  M+1:M-ns => ns valores

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

Npad  = 1024; % zero-padding
wpad  = (0:Npad-1)*(2*pi/Npad);

% Espectro de amplitud de la respuesta en frecuencia teorica

H0  = sqrt(1./(1+rho^2-2*rho*cos(wpad-wp))); % hasta Nyquist
H0  = abs(H0)*norm; % normalizo para comparar con las muestras de |H(w)|

Hpad = abs(fft(h,Npad));   % Espectro de amplitud del filtro
S    = abs(fft(s,Npad));   % Espectro de amplitud de la senal
SF   = abs(fft(sf,Npad));  % Espectro de amplitud de la senal filtrada

% etiquetas en frecuencia angular digital
xt      = 0:pi/2:2*pi;
xtlabel = ["0";"{/Symbol p}/2";"{/Symbol p}";"3{/Symbol p}/2";"2{/Symbol p}"];

 
figure(1)
subplot(2,1,1)
hold on
plot(wpad,abs(H0),'k;|H(\omega)| teorico;')
plot(w,abs(H),'ro;|H_k|;','markersize',4)
plot(wpad,abs(Hpad),'g;|H_k| zero-padded;')
hold off
ylabel('|H(\omega)|')
xlabel('\omega')
set(gca,"xtick",xt,"xticklabel",xtlabel);
axis([0,2*pi,]);
title("Filtro de un polo")
subplot(2,1,2)
stem( (-M:M),h,"r");
ylabel('h_n')
xlabel('n')
axis([-M,M,]);


figure(2)

subplot(2,2,1)
plot(n,s,'b')
ylabel('s_n')
xlabel('n')
axis([0 ns-1,-1 +1])
title("Filtro de un polo aplicado")

subplot(2,2,2)
hold on
plot(n,s1,'k')
plot(n,sf,'r')
hold off
ylabel('sf_n')
xlabel('n')
axis([0 ns-1,-1 +1])

maxS=(1/ns)*max(abs(S));
subplot(2,2,3)
hold on
plot(wpad,(1/ns)*abs(S),'b')
plot(wpad,maxS*abs(Hpad)/max(abs(Hpad)),'g')
hold off
ylabel('|S(\omega)|')
xlabel('\omega')
set(gca,"xtick",xt,"xticklabel",xtlabel);
axis([0,2*pi,0,maxS]);

subplot(2,2,4)
plot(wpad,(1/ns)*abs(SF),'r')
ylabel('|Sf(\omega)|')
xlabel('\omega')
set(gca,"xtick",xt,"xticklabel",xtlabel);
axis([0,2*pi,0,maxS]);



