%% =================================================================
%% CORRIMIENTOS LINEALES DE FASE
%% =================================================================

clear all; %  Limpiar todas las variables definidas

% Redefinimos algunas opciones por defecto
set (0, 'DefaultAxesFontname', 'Times New Roman') 
set (0, 'DefaultAxesFontsize', 12) 
set (0, 'DefaultTextFontname', 'Times New Roman') 
set (0, 'DefaulttextFontsize', 12)

h = [ 1  2  3  4  5  4  3  2  1 ];   % fase lineal

N = length(h);
M = (N-1)/2;

% Calculo TDF
H  = fft(h);

% La hago de fase 0 adelantando M unidades de tiempo (tau == -M)
tau = - M;
for k = 0 : N-1
  H0(k+1) = H(k+1) * exp(-1i*2*pi/N * tau * k);
end

% Vuelvo al tiempo y encuentro la señal adelantada
h0 = real(ifft(H0));

% Corrimiento para adelantar h0 una unidad de tiempo (tau = -1)
tau = -1;
for k = 0:N-1
  H1(k+1) = H0(k+1) * exp(-1i*2*pi/N * tau * k);
end

% Vuelvo al tiempo y encuentro la señal retardada
h1 = real(ifft(H1));


%Calculo espectro de amplitud (normalizados)
A   = abs(H)   /  max(H);
A0 = abs(H0) / max(H);
A1 = abs(H1) / max(H);

% Calculo los espectros de fase
P   = unwrap(  angle(H)  );
P0 = unwrap( angle(H0) );
P1 = unwrap( angle(H1) );


% ________________________________________________________________

% Graficos

clf; 

% Vector frecuencia angular digital
dw =  2*pi/N;
w   = (0:N-1)*dw;

% Tics para los ejes en frecuencia (espectros de amplitud y fase)
xt      = 0:pi:2*pi;
xtlabel = [ '0' ; '{/Symbol p}' ; '2{/Symbol p}' ];

% Respuestas impulsivas

subplot(3,3,1)
stem(h,'bo');
axis([0 N+1 -1 M+3]);
title('Fase lineal');
xlabel('n')
ylabel('h_n')

subplot(3,3,2)
stem(h0,'ro');
axis([0 N+1 -1 M+3]);
title('Fase cero');
xlabel('n')

subplot(3,3,3)
stem(h1,'go');
axis([0 N+1 -1 M+3]);
title('Fase lineal (tau = -1)');
xlabel('n')

% Espectros de amplitud

subplot(3,3,4)
stem(w,A,'b-o');
axis([-dw 2*pi -0.2 1.2]);
ylabel('|H(\omega)|')
xlabel('\omega');
set(gca,'xtick',xt,'xticklabel',xtlabel);

subplot(3,3,5)
stem(w,A0,'r-o');
axis([-dw 2*pi -0.2 1.2]);
xlabel('\omega');
set(gca,'xtick',xt,'xticklabel',xtlabel);

subplot(3,3,6)
stem(w,A1,'g-o');
axis([-dw 2*pi -0.2 1.2]);
xlabel('\omega');
set(gca,'xtick',xt,'xticklabel',xtlabel);

% Espectros de fase

subplot(3,3,7)
stem(w,P,'b-o');
axis([-dw 2*pi -M*2*pi M*2*pi]);
ylabel('\phi(\omega)')
xlabel('\omega');
set(gca,'xtick',xt,'xticklabel',xtlabel);

subplot(3,3,8)
stem(w,P0,'r-o');
axis([-dw 2*pi -M*2*pi M*2*pi]);
xlabel('\omega');
set(gca,'xtick',xt,'xticklabel',xtlabel);

subplot(3,3,9)
stem(w,P1,'g-o');
axis([-dw 2*pi -M*2*pi M*2*pi]);
xlabel('\omega');
set(gca,'xtick',xt,'xticklabel',xtlabel);

%___________________________________________________________________

% Guardamos el gráfico de la pantalla en un postscript
print('phase_shifts.eps','-depsc');
