%% Ejemplo de Sistema Discreto Equivalente
% Copyright (c) 2024 G.T.A.S. (Grupo Tratamiento Avanzado de Señal). Universidad de Cantabria.
%
% Ejemplo para comprender el Sistema Discreto Equivalente de un sistema de comunicaciones digitales. Recuerda lo visto en el aula: la frecuencia de muestreo a la que se trabaja es igual a la velocidad de símbolo, es decir, f_s=R_s.
% Por lo tanto, en este tipo de simulaciones cada muestra de la simulación equivale a un símbolo.
%
%% Parámetros del Transmisor

N = 1000;    % Longitud de la secuencia de símbolos a transmitir
M = 16;      % Número de símbolos en la constelación M-QAM: p.e. M=16 en 16-QAM

%% Transmisor
% Generamos con Matlab los bits en el transmisor y los mapeamos en símbolos de una constelación de M símbolos, es decir, con \log_2(M) bits/símbolo. Podemos observar el inicio de la secuencia de símbolos y dibujar la constelación transmitida.

bits = randi([0 1],N*log2(M),1);        % Generación de N*log2(M) bits aleatorios
a = qammod(bits,M,'InputType', 'bit');  % Secuencia de símbolos de una modulación M-QAM
disp(a(1:10))                           % Muestra en pantalla los 10 primeros símbolos
figure(1);plot(a,'o');axis square;grid;title('Constelación transmitida');

% Nota: La instrucción plot(a,'o') cuando a es una variable compleja es equivalente a plot(real(a),imag(a),'o')
%
%% Canal
% Definimos el canal discreto equivalente. Para empezar consideramos un canal ideal, es decir, un filtro con respuesta al impulso igual a una delta de Kronecker: h_{ce}[n]=\delta[n]. Posteriormente podrás modificarlo

hce = [1];    % Respuesta al impulso del canal discreto equivalente

%% Símbolos Recibidos
% Los símbolos recibidos resultan de la convolución de la secuencia de símbolos transmitidos y la respuesta al impulso del canal discreto equivalente. Por el momento no añadimos ruido ni ningún otro tipo de efecto.

r = conv(hce,a);
disp(r(1:10))
figure(2);plot(r,'.');axis square;grid;title('Observables');

%% Prueba con otros canales
% Ejecuta de nuevo esta sección modificando la respuesta al impulso del canal.
% Prueba, por ejemplo, con un canal con la siguiente respuesta al impulso:
%             hce = [1  0.1*exp(1i*pi/3)  0.03*exp(1i*2*pi/7)];
% A la vista de los observables, ¿qué está ocurriendo?. Prueba con otros canales.
%
% En las siguientes secciones puedes probar a incorporar otros efectos al sistema de comunicaciones.
%
%% Ruido
% Añade ruido AWGN a los símbolos recibidos de la siguiente forma:

SNR = 20;    % Relación Señal a Ruido en dB's
z = awgn(r,SNR,'measured');   % Añadimos Ruido AWGN a los símbolos recibidos
disp(z(1:10))
figure(3);plot(z,'.');axis square;grid;title('Observables con ruido');

% Puedes observar el efecto en la constelación de los observables de diferentes SNR's sin más que ejecutar esta sección modificando el valor de la variable SNR.
%
%% Error de fase
% Puedes incluir en la simulación el efecto que tendría un error de fase entre el transmisor y el receptor. Simplemente debes multiplicar la señal recibida por un número complejo de módulo unidad, es decir,
% z[n]=r[n]e^{j\phi}.

Desfase = 30;    % Error de fase en grados
z = exp(1i*Desfase*pi/180)*r; % Aplicamos el desfase a los observables
disp(z(1:10))
figure(4);plot(z,'.');axis square;grid;title('Observables con desfase');

%
%% Error de frecuencia
% También puedes añadir erores de frecuencia simplemente multiplicando la secuencia de símbolos recibidos por una señal exponencial compleja de cierta frecuencia, es decir,
% z[n]=r[n]e^{j2\pi \omega n}.

Error_frec = 0.00005;    % Error frecuencial
z = r.*exp(1i*2*pi*Error_frec*(0:length(r)-1)'); % Multiplicamos por una señal exponencial compleja de frecuencia Error_frec
disp(z(1:10))
figure(5);plot(z,'.');axis square;grid;title('Observables con error frecuencial');

% ¿Qué le está pasando a la constelación de observables?. Quizás modificar el error de frecuencia (haciéndolo el doble o la mitad del valor actual) te ayude a comprender lo que ocurre.
%
%% Prueba con ...
% Puedes intentar crear el código matlab que incluya todos los fenómenos tratados, es decir, que la secuencia de observables siga la expresión
% z[n]=r[n]e^{j\phi}e^{j2\pi \omega n}+w[n],
% siendo r[n]=a[n]*h_{ce}[n] la señal recibida (tras atravesar el canal) y w[n] la secuencia de ruido AWGN (recuerda que en Matlab estamos usando la función awgn)
%  Simplemente completa las líneas comentadas en el siguiente código y modifica los distintos parámetros a tu gusto.

N = 1000;    
M = 16;      

bits = randi([0 1],N*log2(M),1);
a = qammod(bits,M,'InputType', 'bit');

hce = [1  0.01*exp(1i*pi/3)  0.03*exp(1i*2*pi/7)]; 
SNR = 40;
Desfase = 10;
Error_frec = 0.00009;
% r = 
% z =
figure(6);plot(z,'.');axis square;grid;title('Observables');