% Sincronismo de fase ciego (Potencia M-ésima)

%   Copyright (c) 2013 G.T.A.S. (Grupo Tratamiento Avanzado de Señal)
%                 DICOM Universidad de Cantabria
%                 Jesús Ibáñez

%close all

N=400;      % Número de símbolos de la simulación
SNR=20;     % SNR en dBs
desfase=30; % Desfase introducido por el canal en grados
gamma=0.01; % Constante de paso del algoritmo

M=4;    % Tipo de modulación M-PSK

modulador=modem.pskmod('M',M,'InputType', 'Bit');
demodulador=modem.pskdemod('M',M); 

bits=randi([0 1],N*log2(M),1);  % bits a transmitir
a=modulate(modulador,bits);  % Símbolos
z=a*exp(j*pi/180*desfase)+(10^(-SNR/20))/sqrt(2)*(randn(size(a))+j*randn(size(a)));    % símbolos recibidos (con error de fase y AWGN)

z2=zeros(size(z));
a_dec=zeros(size(z));
e=zeros(size(z));
phi=zeros(size(z));

phi(1)=0; % Inicialización de la fase estimada
for n=1:N
    z2(n)=z(n)*exp(-j*phi(n)); % Observable: símbolos recibidos con fase corregida
    a_dec(n)=exp(1i*(demodulate(demodulador,z2(n))*2*pi/M));
    e(n)=imag(z(n).^M*exp(-j*M*phi(n))); % Error
    phi(n+1)=phi(n)+gamma*e(n); % Nueva fase estimada
end

subplot(234);plot(real(a),imag(a),'or');grid on;axis([-2 2 -2 2]);axis 'square';title('Constelación TX');
subplot(211);plot(desfase-phi*180/pi);grid;xlabel('n (nº símbolo)');ylabel('error fase (º)');title('Evolución del error de fase')
for n=1:N
    subplot(235);plot(real(z(n)),imag(z(n)),'.b');grid on;axis([-2 2 -2 2]);axis 'square';hold on;title('Observables');
    if n<150
        subplot(236);plot(real(z2(n)),imag(z2(n)),'.b');
    else
        subplot(236);plot(real(z2(n)),imag(z2(n)),'.r');
    end
    grid on;axis([-2 2 -2 2]);axis 'square';hold on;title('Observables Corregidos');
    subplot(211);plot(desfase-phi(1:n)*180/pi,'r');axis([0 N -90 90]);grid;xlabel('n (nº símbolo)');ylabel('error fase (º)');title('Evolución del error de fase');
    pause(0.01)
end
subplot(211);hold off
subplot(235);hold off
subplot(236);hold off

