%% =================================================================
%% TEOREMA DEL CORRIMIENTO DE FASE
%% =================================================================

clear; clf;

% Parámetros del operador
M  = 10;        % OBS: M debe ser para para que M/2 sea entero 
N  = 2*M + 1;
dw = 2*pi/N;

% Defino el espectro de amplitud como un cajón (pasabajos)
Aamp = [ones(1,M/2+1) zeros(1,M) ones(1,M/2)];

% Defino la respuesta en frecuencia (compleja)
A = complex(Aamp) + 1.e-8;

% Calculo el espectro de fase
Afase = angle(A);

% Calculo el operador en tiempo (respuesta impulsiva de A)
a = real(ifft(A));

% Realizo el corrimiento de fase lineal
% (atraso la señal en tiempo M muestras)
for k = 1:N
  wk   = dw*(k-1);
  B(k) = A(k)*exp(-1i*wk*M);
end

% Calculo el espectro de amplitud
Bamp = abs(B);

% Calculo el espectro de fase
Bfase = -angle(B);

% Hago "phase unwrapping" para poder ver que el cambio de fase es
% efectivamente lineal. Sino, la fase queda entre -pi y pi
B2fase = unwrap(Bfase);

% Calculo el operador en tiempo (respuesta impulsiva de B)
b = real(ifft(B));

% =================================================================
% GRÁFICOS
% =================================================================

% Grafico los espectros de amplitud

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

subplot(3,3,1);
plot((0:N-1)*dw,Aamp,';|A_k|;b-o');
title('Espectros')
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt);
axis([0,2*pi,-0.2,1.3]);

subplot(3,3,4);
plot((0:N-1)*dw,Bamp,';|B_k|;b-o');
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt);
axis([0,2*pi,-0.2,1.3]);

% Grafico los espectros de fase

yt      = -pi:pi:pi;
ytlabel = ["-{/Symbol p}";"0";"{/Symbol p}"];

subplot(3,3,2);
plot((0:N-1)*dw,Afase,';{/Symbol F}(A)_k;r-o');
title('Fases')
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt,"yticklabel",ytlabel);
axis([0,2*pi,-4,4]);

subplot(3,3,5);
plot((0:N-1)*dw,Bfase,';{/Symbol F}(B)_k;r-o');
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt,"yticklabel",ytlabel);
axis([0,2*pi,-4,4]);

yt      = 0:10*pi:20*pi;
ytlabel = ["0";"10{/Symbol p}";"20{/Symbol p}"];

subplot(3,3,8);
plot((0:N-1)*dw,B2fase,';{/Symbol F}(B2)_k;r-o');
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt,"yticklabel",ytlabel);
axis([0,2*pi,-5,70]);

% Grafico los operadores en tiempo

subplot(3,3,3);
plot((0:N-1),a,';a_n;g-o');
title('Respuestas impulsivas')
axis([-1,N,-0.2,0.8]);

subplot(3,3,6);
plot((0:N-1),b,';b_n;g-o');
axis([-1,N,-0.2,0.8]);

% Guardo el gráfico en un postscript

print('cambio_de_fase.eps','-depsc'); % guardo el gráfico
