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

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

%close all

%% Parámetros de la simulación
N=200;      % Número de símbolos de la simulación
SNR=20;     % SNR en dBs
desfase=35; % Desfase introducido por el canal en grados
gamma=0.0005; % Constante de paso del algoritmo

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

%% Generación de secuencia de símbolos transmitidos y recibidos
bits = randi([0 1],N*log2(M),1);        % bits a transmitir
a = qammod(bits,M,'InputType', 'Bit');  % 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)

%% Algoritmo de sincronización de fase ciego
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)=qammod(qamdemod(z2(n),M),M);   % Símbolos decididos de la constelación M-QAM
    e(n)=imag(z(n).^M*mean(qammod(0:M-1,M).^M)*exp(-j*M*phi(n))); % Error
    phi(n+1)=phi(n)+gamma*e(n); % Nueva fase estimada
end

%% Representación gráfica animada muestra a muestra
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 abs(mean(desfase-phi(max(1,n-19):n)*180/pi))>2
        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.002)
    drawnow
end
subplot(211);hold off
subplot(235);hold off
subplot(236);hold off

