%********** It Simulates OFDM Systems **********%
clear all
close all
clc

%********** Parameters **********%
fs=12e6;%Hz,sample frequency
L=5;%multipaths number
fm=150; % maximun Doppler Frequency 
N=1024;
L_data=768;
L_pilot=256;
L_cyclic_prefix=226;
L_symbol=2;
L_channel=L_cyclic_prefix;
Pattern=ones(1,N);
Pattern(1:4:end)=4;
Pattern=repmat(Pattern,1,L_symbol);
QAM_Pilot=zeros(size(Pattern));

% ******************** Transmitter ********************%
user_entry = Menu('Please select the primary channel modulation mode',...
                  '                              16QAM                              ',...
                  '                              64QAM                              ',...
                  '                              QPSK                               ',...
                  '                              EXIST                              ');
if user_entry==1,%16QAM
    QAM_data=rand_gen(L_data*L_symbol,16);
    M=16;
elseif user_entry==2 %64QAM
    QAM_data=rand_gen(L_data*L_symbol,64);
    M=64;
elseif user_entry==3 %QPSK
    QAM_data=rand_gen(L_data*L_symbol,4);
    M=4;
else
    return
end;
% Insert pilots and data
Pilot_pattern=find(Pattern==4);
QAM_Pilot(Pilot_pattern)=1;
Pilot_data=QAM_Pilot(Pilot_pattern);
clear M2047_Sequence;
QAM_pattern=find(Pattern==1);
QAM_Pilot(QAM_pattern)=QAM_data;
QAM_Pilot=reshape(QAM_Pilot,N,L_symbol).';
% IFFT
Trans=[];
for i=1:L_symbol,
    IDFT_TEMP=ifft(QAM_Pilot(i,:));
    Trans=[Trans;IDFT_TEMP];
end;
clear IDFT_TEMP;
Trans=[Trans(:,N-L_cyclic_prefix+1:N),Trans];
   

%********** Pass through channel **********%
tau=[1,5,35,50,60].';
PowerProfile=[0,-5,-10,-20,-30].';
h_temp=MJakesCH(PowerProfile,fm,fs,1,L_symbol,tau);
h=[h_temp.',zeros(L_symbol,L_channel-size(h_temp.',2))];
Ch_data=[];
for i=1:L_symbol,
    Ch_data=[Ch_data;conv(Trans(i,:),h(i,:))];
end;
Ch_data=awgn(Ch_data,20,'measured');

%********** Receiver **********%
% Remove the CP data
Receive=Ch_data(:,L_cyclic_prefix+1:N+L_cyclic_prefix);
% FFT
DFT_data=fft(Receive.').';
temp_DFT_data=reshape(DFT_data.',1,L_symbol*N);
DFT_QAM=temp_DFT_data(QAM_pattern);
DFT_QAM=reshape(DFT_QAM,L_data,L_symbol).';
DFT_pilot=temp_DFT_data(Pilot_pattern);
DFT_pilot=reshape(DFT_pilot,L_pilot,L_symbol).';
% Demodulation
H=fft(h.',N).';
[h_est_DFT]=H; % channel estimation
h_est_DFT=reshape(h_est_DFT.',1,N*L_symbol);
h_est_DFT_QAM=h_est_DFT(QAM_pattern);
h_est_DFT_QAM=reshape(h_est_DFT_QAM.',L_data,L_symbol).';
Demod_data=DFT_QAM./h_est_DFT_QAM; % equalization
Demod_data=reshape(Demod_data.',1,[]);
plot(real(Demod_data(1,:)))
hold on
plot(real(QAM_data(1,:)),'r:')