%% Extensión Cíclica en OFDM
% Copyright (c) 2018 G.T.A.S. (Grupo Tratamiento Avanzado de Señal). Universidad de Cantabria.
%
% Dos condiciones que debe que cumplir el tiempo de guarda de una señal OFDM:
% Para evitar completamente la ISI en un sistema OFDM es necesario que el tiempo de guarda sea mayor que la dispersión máxima del retardo del canal.
% Para no perder la ortogonlaidad 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.
%
%% 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
k = 4;      % Tipo de modulación (k bits/simbolo): BPSK k=1; QPSK k=2; 16QAM k=4; 64QAM k=6

% 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:

modulador = modem.qammod('M',2^k,'InputType', 'Bit'); % Mapeador de símbolos
s = modulate(modulador,randi([0 1],N*k,Ns));  % Símbolos
x = ifft(s,N);  % 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);    % 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 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 multicamino: h_c[n]=\sum^{L}_{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;
% 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>N_g 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

% 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('Subportadora');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('Subportadora');ylabel('Fase (º)')
    subplot(122);h = plot(real(z(portadora,:)),imag(z(portadora,:)),'.r','MarkerSize',8);
    escala = 2*sqrt(2^k);axis(escala*[-1 1 -1 1]);axis square;grid
    title(sprintf('Constelacion subportadora %i. Ng=%i, extension ciclica %i',portadora-1, Ng, ext_ciclica));
    set(handle,'Visible','on')
    pause(0.2)
end