%% Extensión Cíclica en OFDM
% Copyright (c) 2023 G.T.A.S. (Grupo Tratamiento Avanzado de Señal). Universidad de Cantabria.
%
% Dos condiciones que debe cumplir el tiempo de guarda de una señal OFDM:
% Para evitar completamente la ISI (Inter-Symbol Interference) en un sistema OFDM es necesario que el tiempo de guarda sea mayor que la dispersión temporal máxima del canal.
% Para no perder la ortogonalidad entre subportadoras en presencia de canales con memoria, es preciso que durante el tiempo de guarda se replique cíclicamente el símbolo OFDM (extensión cíclica).
% En este ejemplo se comprueba el efecto del canal sobre la señal OFDM cuyo tiempo de guarda puede usar o no extensión cíclica y su duración puede ser mayor o menor que la memoria del canal.
% Para facilitar la visualización de las constelaciones y el impacyo que sobre ellas tiene la memoria del canal y el uso o no de una guarda cíclica apropiada, en esta simulación no se añade ruido en el receptor.
%
%% Parámetros del transmisor
% Primero definimos los parámetros de la señal OFDM transmitida. Puedes modificarlos si deseas probar con más o menos subpoortadoras o con otro tipo de modulación:

N = 16;     % Número de subportadoras (tamaño FFT)
Ns = 500;   % Número de símbolos OFDM a generar
M = 4;      % Tipo de modulación (log2(M) bits/simbolo): BPSK M=2; QPSK M=4; 16QAM M=16; 64QAM M=64

% Indica el número de muestras, N_g, del tiempo de guarda  y si se empleará extensión cíclica o no:

Ng = 4;           % Número de muestras del tiempo de guarda cíclica
ext_ciclica = 1;  % poner a 1 para extension ciclica
% ext_ciclica = 0;  % poner a 0 para tiempo de guarda sin señal

%
%% Transmisor
% Generamos la señal OFDM transmitida y la dibujamos:

bits = randi([0 1],N*log2(M),Ns);       % bits a transmitir. Cada columna corresponde a un símbolo OFDM
s = qammod(bits,M,'InputType', 'Bit');  % símbolos M-QAM a transmitir
x = ifft(s);             % Matriz de muestras temporales de la señal OFDM. Cada columna corresponde a un símbolo OFDM 
xg = [ext_ciclica*x(end-Ng+1:end,:);x]; % Insercción de la guarda
xt = xg(:);     % Vector con las muestras de la señal OFDM

Nsimbolos = min(5,Ns);    % para mostrar, a lo sumo 5, símbolos en el plot
n = 1:Nsimbolos*(Ng+N);    
figure(1);subplot(211);plot(real(xt(n)));set(gca,'XTick',(0:Nsimbolos)*(N+Ng));grid on
legend('Señal I');xlabel('muestras');ylabel('amplitud');title(sprintf('Señal transmitida (%i primeros símbolos)',Nsimbolos))
subplot(212);plot(imag(xt(n)));set(gca,'XTick',(0:Nsimbolos)*(N+Ng));grid on;legend('Señal Q');xlabel('muestras');ylabel('amplitud')

%
%% Define el canal
% Propón la respuesta al impulso de un canal equivalente paso bajo. Algunos ejemplos pueden ser:
% Canal ideal: h_c[n]=\delta[n], por ejemplo, en Matlab  hc = [ 1 ];     % canal sin memoria
% Canal de dos rayos: h_c[n]=\delta[n] + \alpha \delta[n-n_0], por ejemplo, en Matlab  hc = [ 1 ,  0.3-1i*0.5 ];    % canal con memoria igual a 1
% Canal multicamino: h_c[n]=\sum^{L_h}_{k=0}{\alpha_k\delta[n-k]}, por ejemplo, en Matlab  hc = [1*exp(-1i*0.5) 0 -1*exp(1i*0.9) 0.4*exp(1i*0.3)]/2;   % canal con memoria igual a L_h=3
% Puedes probar con el canal que desees. Simplemente añade o quita términos a la respuesta al impulso del canal. Recuerda que si L_h=length(h_c)-1>N_g, entonces aparecerá ISI en el receptor.

hc = [1*exp(-1i*0.5) 0 -1*exp(1i*0.9) 0.4*exp(1i*0.3)]/2;    % Respuesta al impulso del canal (tiene memoria)
% hc = [1*exp(-1i*0.5)];

% Puedes observar la respuesta en frecuencia del canal que has elegido y comprobar si es o no selectivo en frecuencia. En el eje horizontal de la figura, la frecuencia expresada en Hz se corresponde con el número de subportadora, pues le indicaremos a la función freqz que la frecuencia de muestreo es Ny, por tanto, resulta que la distancia entre subportadoras es 1 Hz.

figure(2);freqz(hc,1,32*N,N,'whole');    % Representación de la respuesta en frecuencia del canal propuesto

%
%% Recepción
% Filtramos la señal transmitida con el canal para obtener la señal recibida (en esta simulación no añadimos ruido) y procesamos la señal recibida hasta obtener los observables de las distintas subportadoras.

r = filter(hc,1,xt);    % Señal transmitida filtrada por el canal
y = reshape(r,N+Ng,Ns);
y = y(Ng+1:end,:);    % Eliminación de guardas
z = fft(y,N);         % FFT para obtener los símbolos recibidos en cada subportadora

% Ahora puedes observar la constelación recibida en cada una de las subportadoras. Aparecerá una ventana emergente con una figura de Matlab, que muestra la amplitud y fase de la respuesta del canal a cada frecuencia subportadora y la constelación resultante. Si deseas que la animación sea más lenta incrementa el argumento de la función pause (o no le pases argumento, en cuyo caso espera hasta que pulsas una tecla)

[Hc,f] = freqz(hc,1,32*N,N,'whole');    % Respuesta del canal
[Hck,fk] = freqz(hc,1,N,N,'whole');     % Respuesta del canal en las frecuencias subportadoras

for portadora=1:N      %% Bucle para dibujar la constelación recibida en cada subportadora
    handle = figure(3);
    subplot(221);plot(f,abs(Hc),'b',fk(portadora),abs(Hck(portadora)),'ro');grid on;xlabel('Nº Subportadora (frecuencia)');ylabel('Amplitud');title('Respuesta en frecuencia del canal')
    subplot(223);plot(f,180/pi*angle(Hc),'b',fk(portadora),180/pi*angle(Hck(portadora)),'ro');grid on;xlabel('Nº Subportadora (frecuencia)');ylabel('Fase (º)')
    subplot(122);h = plot(real(z(portadora,:)),imag(z(portadora,:)),'.r','MarkerSize',16);
    escala = 2*sqrt(M);axis(escala*[-1 1 -1 1]);axis square;grid
    title(sprintf('Observables: subport %i. Ng=%i, ext cicl=%i',portadora-1, Ng, ext_ciclica));
    set(handle,'Visible','on')
    pause
end

%% Prueba tú ahora
% Prueba a modificar el tiempo de guarda (N_g) de la señal OFDM y la respuesta al impulso del canal discreto equivalente. Ejecuta de nuevo el Live Script con y sin extensión cíclica y observa en qué circunstancias aparece (o no) la ISI.
% A modo de ejemplo, puedes intentar los siguientes casos:
%
% 1) Caso inicial: 4 muestras de guarda con extensión cíclica y canal multicamino con memoria 3:
% Ng = 4; ext_ciclica = 1; hc = [1*exp(-1i*0.5) 0 -1*exp(1i*0.9) 0.4*exp(1i*0.3)]/2;
% Observa la constelación recibida de cada subportadora. Deberás comprobar que las constelaciones son perfectas (recuerda que no hay ruido) salvo por la atenuación/desfase que introduce el canal a cada subportadora (observa la respuesta frecuencial del canal en las figuras de la izquierda). En esta situación no hay ISI. ¿por qué?
%
% 2) Prueba a reducir el número de muestras de guarda
% Ng = 3; ext_ciclica = 1; hc = [1*exp(-1i*0.5) 0 -1*exp(1i*0.9) 0.4*exp(1i*0.3)]/2;
% Comprueba que los observables de las constelaciones recibidas, además de sufrir atenuación y desfase, se dispersan (y no es ruido). ¡¡Es ISI!! y también ICI (Inter-Carrier Interference).
% Juega con el número de muestras de guarda y con la longitud (memoria) del canal (añade a voluntad elementos adicionales al vector hc) para comprobar que, en efecto, es necesario que la duración de la extensión cíclica sea mayor que la memoria del canal.
%
% 3) Haz que el transmisor emplee guarda sin extensión cíclica. Partiendo de un caso en el que hayas comprobado que las constelaciones recibidas son correctas (no hay ISI/ICI) modifica la variable ext_ciclica = 0. Observa que la señal transmitida es nula durante los tiempos de guarda. ¿Qué ocurre con las constelaciones recibidas?.
%
%
%