%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%      uncoded OFDM link simulator
%%      system includes: 
%%          OFDM modulator, OFDM demodulator, 
%%          Channel estimator and equalizer
%%          some optional modules: clip and filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc; clear;
%% channel init
h_delay = [2, 4, 6, 12, 25];
h_amp = [0, -3, -10, -5, -12];
h = zeros(1,25);
h(h_delay) = 10.^(h_amp/10).*(randn(1,5) + j*randn(1,5));

%% system init
Sym_Num = 32;
ML = 2; %% modulation level
load P_sym.mat;
error_sum = zeros(1,10);
error_sum1 = zeros(1,10);
D_snr = zeros(1,10);

for snr_ind = 1:10
    SNR = 2*snr_ind;
    for loop = 1:1000
%% source
code = (sign(randn(Sym_Num, 768*ML))+1)/2;

%% subcarrier modulation
%% QPSK constellation map rule
%% 00 --> 1+j; 01 --> -1+j; 11 --> -1-j; 10 --> 1-j;
p_code = 1 - 2*code;
QPSK_sym1 = p_code(:, 1:768) + j*p_code(:, 769:1536); %% QPSK mod
QPSK_sym = [P_sym; QPSK_sym1];

%% subcarrier map and OFDM modulate
map_index = [3:386, 638:1021];
map_sym = zeros(Sym_Num+1,1024);
map_sym(:,map_index) = QPSK_sym;
map_sym = map_sym.';
td_sym = ifft(map_sym);


%% add CP
td_sym_cp = [td_sym(769:1024,:); td_sym];
tx_sym = reshape(td_sym_cp, 1, []);
% tx_sym1 = clip_filter(tx_sym, 6, 1, 0);
% SNR = 10;
%% pass a channel

rx_sym = tx_sym;
% SNR = 30;
rx_sym = awgn(rx_sym, SNR, 'measured', 'dB');
% rx_sym = rx_sym(1:42240);
rx_td_sym_cp = reshape(rx_sym, 1280, Sym_Num+1);

%% eliminate CP
rx_td_sym = rx_td_sym_cp(257:1280,:);

%% OFDM demodulate and subcarrier demap
fd_sym = fft(rx_td_sym);
demap_sym = fd_sym(map_index,:);

%% channel estimation
fd_h = fft([h, zeros(1,1024-25)]);
fd_h = fd_h(map_index).';
Rx_P_sym = demap_sym(:,1);
fd_h1 = Rx_P_sym./(P_sym.');

%% channel equalize
for i = 1:Sym_Num
    eq_sym(:,i) = demap_sym(:,i+1)./fd_h1;
    sym(:,i) = demap_sym(:,i+1);
end
% plot(eq_sym,'.')
noise1 = sym - QPSK_sym1.';
noise2 = eq_sym - QPSK_sym1.';
noise_p1 = sum(sum(abs(noise1).^2));
noise_p2 = sum(sum(abs(noise2).^2));
delta_snr = 10*log10(noise_p1/noise_p2);
D_snr(snr_ind) = D_snr(snr_ind) + delta_snr;

%% data decision// for symbol without equalizer
dec_sym1 = [real(sym); imag(sym)];
dec_sym1 = dec_sym1 < 0;
dec_sym1 = dec_sym1.';

%% error count
error = code - dec_sym1;
error_sum1(snr_ind) = error_sum1(snr_ind) + sum(sum(abs(error)));

%% data decision// hard decision
dec_sym = [real(eq_sym); imag(eq_sym)];
dec_sym = dec_sym < 0;
dec_sym = dec_sym.';

%% error count
error = code - dec_sym;
error_sum(snr_ind) = error_sum(snr_ind) + sum(sum(abs(error)));


    end
end
data_sum = 1000*Sym_Num*768*ML;
error_rate1 = error_sum1/data_sum;
error_rate = error_sum/data_sum;

semilogy([2:2:20], error_rate1);
hold
semilogy([2:2:20], error_rate,'r');
