% Marco Driusso, University of Trieste, marco.driusso@phd.units.it
% 03-04-2014
% Fulvio Babich
% 12/11/2016
% ============================================================================================
% Script containing some examples about the simulation of the OFDM technique:
% 0) performance of a BPSK in a fading channel
% 1) performance of a BPSK in a 1xN MIMO channel using MRC at RX
% 2) samples of a OFDM signal using IFFT and comparison with the CT waveform
% 3) performance of a BPSK modulated OFDM system in a freq selective channel
% 4) performance of a BPSK modulated OFDM system in a freq selective 1xN channel with MRC at RX
% 5) performance of a BPSK modulated OFDM system in a freq selective 1xN channel with MRC at RX with time variant channel (the channel is not constant during the OFDM frame)
% (Fulvio Babich)
% 6) PAPR evaluation OFDM
% 7) performance of a BPSK modulated SC-FDM system in a freq selective channel
% 8) PAPR evaluation SC-FDM 
% ============================================================================================

clear all;
%close all;
clc;

mode = 6;           % select which of the examples do you want to run
N_test = 200;       % number of simulations for BER calculations or for PAPR evaluation
N_fft = 256;         % number of samples for the DFT/IDFT
M_dft = 192;         % used subcarriers for SC-FDM (even)
N_bits = N_fft;     % number of bits to transmit each timreal(z))e (BPSK mod)
M_bits = M_dft;     % number of bits to transmit each time SC-FDM (BPSK mod)
N = 2;              % number of receiving antennas
DM=1;               %diversity mode 1: MRC, 2: ECG, 3:SD
N_cp = 16;          % cyclic prefix length
Delta_f = 1.0;      % frequency gap between subcarriers
T = 1.0/Delta_f;    % duration of an OFDM symbol (with no CP)
T_s = T/N_fft;      % sampling time
L_ch = 7;           % channel length (in taps)
xc=0:.1:10;
disp_step=100;
snr_db = -10:2:20;
snr = 10.^(.1.*snr_db);
Nfdt=1e-3;           % normalized maximum doppler spread
j = sqrt(-1);

if mode == 0
% -------------------------------------------------------------------------------------
% Example 0 -> performance of a BPSK in a fading channel
% -------------------------------------------------------------------------------------
    disp('Running example for performance of a QPSK in a fading channel');
    ber_qpsk = zeros(size(snr_db));
    for k=1:length(snr)
        for n=1:N_test
            % random bits
            bits = rand(2, N_bits)>.5;
            % QPSK modulation ("0"->+1, "1"->-1)
            s = 1-2*bits(1,:)+j*(1-2*bits(2,:));
            % Rayleigh uncorrelated fading
            h = sqrt(.5).*(randn(size(s)) + j*randn(size(s)));
            % AWGN noise
            z = sqrt(.5/snr(k)).*(randn(size(s)) + j*randn(size(s)));
            % simultion of noisy reception through fading channel
            r = h.*s + z;
            % ZF equalization
            r_eq = r./h;
            % decoding bits
            dec_bits(1,:) = real(r_eq)<0;
            dec_bits(2,:) = imag(r_eq)<0;
            % calculate BER
            ber_qpsk(k) = ber_qpsk(k)+sum(sum(bits~=dec_bits));
        end
    end
    ber_qpsk = ber_qpsk./(N_test*N_bits*2);
    semilogy(snr_db, ber_qpsk, 'b-o');
    grid on;
    set(gca,'FontName','TimesNewRoman');
    set(gca,'FontSize',12);
    xlabel('SNR [dB]');
    ylabel('BER');
    title('QPSK modulated system in Rayleigh fading channel');
elseif mode == 1
% -------------------------------------------------------------------------------------
% Example 1 -> performance of a BPSK in a 1xN MIMO channel using MRC at RX
% -------------------------------------------------------------------------------------
    switch DM
        case 1
            disp('Running example for performance of a 1xN MIMO channel using MRC at RX');
        case 2
            disp('Running example for performance of a 1xN MIMO channel using EGC at RX');
        case 3
            disp('Running example for performance of a 1xN MIMO channel using SD at RX');
    end
            
    ber_bpsk_mrc = zeros(size(snr_db));
    for k=1:length(snr)
        for n=1:N_test
            % random bits
            bits = rand(1, N_bits)>.5;
            % BPSK modulation ("0"->+1, "1"->-1)
            s = 1-2*bits;
            % Rayleigh uncorrelated fading
            h = sqrt(.5).*(randn(N, length(s)) + j*randn(N, length(s)));
            % AWGN noise
            z = sqrt(.5/(snr(k)/N)).*(randn(N, length(s)) + j*randn(N, length(s)));
            % simultion of noisy reception through fading channel
            r = h.*repmat(s, [N, 1]) + z;
            switch DM
                case 1
                    % perform MRC
                    r_eq = r.*conj(h);
                    r_eq = sum(r_eq,1);
                case 2
                    % perform EGC
                    r_eq = r.*(conj(h)./abs(h));
                    r_eq = sum(r_eq,1);
                case 3
                    % perform SD
                    r_eq = r.*(conj(h)./abs(h));
                    [cv,ci]=max(abs(h));
                    for ir=1:length(ci), r_eq_t(ir)=r_eq(ci(ir),ir); end
                    r_eq=r_eq_t;
            end
            % decoding bits
            dec_bits = real(r_eq)<0;
            % calculate BER
            ber_bpsk_mrc(k) = ber_bpsk_mrc(k)+sum(bits~=dec_bits);
        end
    end
    ber_bpsk_mrc = ber_bpsk_mrc./(N_test*N_bits);
    switch DM
        case 1
            semilogy(snr_db, ber_bpsk_mrc, 'b-o');
        case 2
            semilogy(snr_db, ber_bpsk_mrc, 'b:o');
        case 3
            semilogy(snr_db, ber_bpsk_mrc, 'b-.o');
    end
    set(gca,'FontName','TimesNewRoman');
    set(gca,'FontSize',12);    
    grid on;
    xlabel('SNR [dB]');
    ylabel('BER');
    title(strcat('BPSK modulated system in Rayleigh fading 1x', num2str(N), ' channel with MRC (solid), EGC(dotted), SD (dash-dotted)'));
elseif mode == 2
% -------------------------------------------------------------------------------------
% Example 2 -> samples of a OFDM signal using IFFT and comparison with the CT waveform
% -------------------------------------------------------------------------------------
    disp('Running example for IDFT computed OFDM samples comparison with CT signal');
    % random bits
    bits = rand(1, N_bits)>.5;
    % BPSK modulation ("0"->+1, "1"->-1)
    s = 1-2*bits;
    % time axis
    t = linspace(0, T, 10240);
    % calculate continuous time OFDM signal
    s_ofdm = zeros(size(t));
    for k=0:N_fft-1
        s_ofdm = s_ofdm + s(k+1).*exp(j*2*pi*(k-N_fft/2)*Delta_f*t);
    end
    % calculate OFDM samples with IDFT
    s_shift = fftshift(s);
    s_ofdm_samples = N_fft*ifft(s_shift);
    t_samples = T_s*(0:N_fft-1);
    % show results
    plot(t, real(s_ofdm), 'b-', t_samples, real(s_ofdm_samples), 'ro');
    grid on;
    set(gca,'FontName','TimesNewRoman');
    set(gca,'FontSize',12);    
    legend('continuous time waveform', 'IDFT computed samples');
    xlabel('t/T');
    ylabel('real\{s(t)\}');
elseif mode == 3
% -------------------------------------------------------------------------------------
% Example 3 -> performance of a BPSK modulated OFDM system in a freq selective channel
% -------------------------------------------------------------------------------------
    disp('Running example for performance of a BPSK modulated OFDM system in a freq selective channel');
    ber_bpsk_ofdm = zeros(size(snr_db));
    for k=1:length(snr)
        for n=1:N_test
            % random bits
            bits = rand(1, N_bits)>.5;
            % BPSK modulation ("0"->+1, "1"->-1)
            S = 1-2*bits;
            % restore correct energy
            S = S./sqrt(N_fft);
            % OFDM modulation (shift does not matter in discrete time simulations)
            s_ofdm = N_fft*ifft(S);
            % add CP
            s_cp = [s_ofdm(end-N_cp+1:end), s_ofdm];
            % Rayleigh uncorrelated fading in each channel tap
            h = sqrt(.5).*(randn(1, L_ch) + j*randn(1, L_ch));
            % normalize channel impulse response
            h = h./sqrt(L_ch);
            % AWGN noise
            z = sqrt(.5/snr(k)).*(randn(1, length(s_cp)+L_ch-1) + j*randn(1, length(s_cp)+L_ch-1));
            % simultion of noisy reception through frequency selective fading channel
            r = conv(h, s_cp) + z;
            % discard first and last samples (including CP)
            r = r(N_cp+1:N_cp+N_fft);
            % OFDM demodulation
            R = fft(r);
            % calculate channel in frequency domain (this simulates perfect CIR estimation)
            H = N_fft*fft(h, N_fft);
            % ZF equalization
            R_eq = R./H;
            % decoding bits
            dec_bits = real(R_eq)<0;
            % calculate BER
            ber_bpsk_ofdm(k) = ber_bpsk_ofdm(k)+sum(bits~=dec_bits);
        end
    end
    ber_bpsk_ofdm = ber_bpsk_ofdm./(N_test*N_bits);
    semilogy(snr_db, ber_bpsk_ofdm, 'b-o');
    set(gca,'FontName','TimesNewRoman');
    set(gca,'FontSize',12);
    grid on;
    xlabel('SNR [dB]');
    ylabel('BER');
    title('BPSK modulated OFDM system freq sel channel');
elseif mode == 4
% -------------------------------------------------------------------------------------
% Example 4 -> performance of a BPSK modulated OFDM system in a freq selective 1xN channel with MRC at R
% -------------------------------------------------------------------------------------
    switch DM
        case 1
            disp('Running example for performance of a BPSK modulated OFDM system in a freq selective 1xN channel with MRC at R');
        case 2
            disp('Running example for performance of a BPSK modulated OFDM system in a freq selective 1xN channel with EGC at R');
        case 3
            disp('Running example for performance of a BPSK modulated OFDM system in a freq selective 1xN channel with SD at R');
    end
    ber_bpsk_ofdm_mrc = zeros(size(snr_db));
    for k=1:length(snr)
        for n=1:N_test
            % random bits
            bits = rand(1, N_bits)>.5;
            % BPSK modulation ("0"->+1, "1"->-1)
            S = 1-2*bits;
            % restore correct energy
            S = S./sqrt(N_fft);
            % OFDM modulation (shift does not matter in discrete time simulations)
            s_ofdm = N_fft*ifft(S);
            % add CP
            s_cp = [s_ofdm(end-N_cp+1:end), s_ofdm];
            % Rayleigh uncorrelated fading in each channel tap for each RX antenna
            h = sqrt(.5).*(randn(N, L_ch) + j*randn(N, L_ch));
            % normalize channel impulse response
            h = h./sqrt(L_ch);
            % AWGN noise
            z = sqrt(.5/(snr(k)/N)).*(randn(N, length(s_cp)+L_ch-1) + j*randn(N, length(s_cp)+L_ch-1));
            % simultion of noisy reception through frequency selective fading channel on each antenna
            r = zeros(N, length(z));
            for ant=1:N
                r(ant, :) = conv(h(ant, :), s_cp) + z(ant, :);
            end
            % discard first and last samples (including CP)
            r = r(:, N_cp+1:N_cp+N_fft);
            % OFDM demodulation
            R = fft(r, [], 2);
            % calculate channel in frequency domain (this simulates perfect CIR estimation)
            H = N_fft*fft(h, N_fft, 2);
            % perform MRC
            switch DM
                case 1
                    % perform MRC
                    R_eq = R.*conj(H);
                    R_eq = sum(R_eq,1);
                case 2
                    % perform EGC
                    R_eq = R.*(conj(H)./abs(H));
                    R_eq = sum(R_eq);
                case 3,1
                    % perform SD
                    R_eq = R.*(conj(H)./abs(H));
                    [cv,ci]=max(abs(H));
                    for ir=1:length(ci), R_eq_t(ir)=R_eq(ci(ir),ir); end
                    R_eq=R_eq_t;
            end
            % decoding bits
            dec_bits = real(R_eq)<0;
            % calculate BER
            ber_bpsk_ofdm_mrc(k) = ber_bpsk_ofdm_mrc(k)+sum(bits~=dec_bits);
        end
    end
    ber_bpsk_ofdm_mrc = ber_bpsk_ofdm_mrc./(N_test*N_bits);
    switch DM
        case 1
            semilogy(snr_db, ber_bpsk_ofdm_mrc, 'b-o');
        case 2
            semilogy(snr_db, ber_bpsk_ofdm_mrc, 'b:o');
        case 3
            semilogy(snr_db, ber_bpsk_ofdm_mrc, 'b-.o');
    end
    set(gca,'FontName','TimesNewRoman');
    set(gca,'FontSize',12);
    grid on;
    xlabel('SNR [dB]');
    ylabel('BER');
    title(strcat('BPSK modulated OFDM system in freq sel fading 1x', num2str(N), ' channel with MRC (solid), EGC(dotted), SD (dash-dotted)'));
elseif mode == 5
% -------------------------------------------------------------------------------------
% Example 5 -> OFDM with time variant channel
% -------------------------------------------------------------------------------------
    switch DM
        case 1
            fprintf('Running example for performance of a BPSK modulated OFDM system \nin a freq selective 1xN channel with MRC at RX (time variant channel)\n');
        case 2
            fprintf('Running example for performance of a BPSK modulated OFDM system \nin a freq selective 1xN channel with EGC at RX (time variant channel)\n');
        case 3
            fprintf('Running example for performance of a BPSK modulated OFDM system \nin a freq selective 1xN channel with SD at RX (time variant channel)\n');
    end
    ber_bpsk_ofdm_mrc_tv = zeros(size(snr_db));
    info.fdt = Nfdt;
    info.fc = 2e9;        % carrier frequency [Hz]
    info.T_s = 6.5e-8;      % sampling time
    info.c = 3e8;         % speed of ligth
    info.fd = info.fdt/info.T_s;    % maximum doppler spread
    info.v = info.c*info.fd/info.fc/3.6; % correspondent speed in km/h
    disp(info);
    for k=1:length(snr)
        fprintf('Doing SNR=%g dB..\n', snr_db(k));
        for n=1:N_test
            % random bits
            bits = rand(1, N_bits)>.5;
            % BPSK modulation ("0"->+1, "1"->-1)
            S = 1-2*bits;
            % restore correct energy
            S = S./sqrt(N_fft);
            % OFDM modulation (shift does not matter in discrete time simulations)
            s_ofdm = N_fft*ifft(S);
            % add CP
            s_cp = [s_ofdm(end-N_cp+1:end), s_ofdm];
            % AWGN noise
            z = sqrt(.5/(snr(k)/N)).*(randn(N, length(s_cp)+L_ch-1) + j*randn(N, length(s_cp)+L_ch-1));
            % Rayleigh correlated fading in each channel tap for each RX antenna
            h = sqrt(.5).*sosfadeli(length(z), info.fd, info.T_s, -inf, 64, N*L_ch);
            h = reshape(h.', [length(z), N, L_ch]);
            % normalize channel impulse response
            h = h./sqrt(L_ch);
            % simultion of noisy reception through frequency selective fading channel on each antenna
            r = zeros(N, length(z));
            % padding s_cp for convolution
            s_cp = [zeros(1, L_ch-1), s_cp, zeros(1, L_ch-1)];
            % simulate reception through time variant channel
            % for ant=1:N
            %     r(ant, :) = sum(flipud(permute(h(:, ant, :), [3 1 2])).* ...
            %                     s_cp(repmat(1:length(z), [L_ch, 1])+repmat((0:(L_ch-1))', [1, length(z)]))) ...
            %                 + z(ant, :);
            % end
            % simulate reception through time variant channel (super-vectorized version)
            r = permute(sum(permute(h, [3 1 2]).* ...
                            s_cp(repmat(1:length(z), [L_ch, 1, N])+...
                                 repmat(((L_ch-1):-1:0)', [1, length(z), N]))),...
                        [3 2 1]) + z;
            % discard first and last samples (including CP)
            r = r(:, N_cp+1:N_cp+N_fft);
            % OFDM demodulation
            R = fft(r, [], 2);
            % calculate channel in frequency domain assuming of having estimated it at the first sample
            H = N_fft*fft(permute(h(1, :, :), [2, 3, 1]), N_fft, 2);
            switch DM
                case 1
                    % perform MRC
                    R_eq = R.*conj(H);
                    R_eq = sum(R_eq,1);
                case 2
                    % perform EGC
                    R_eq = R.*(conj(H)./abs(H));
                    R_eq = sum(R_eq,1);
                case 3
                    % perform SD
                    R_eq = R.*(conj(H)./abs(H));
                    [cv,ci]=max(abs(H));
                    for ir=1:length(ci), R_eq_t(ir)=R_eq(ci(ir),ir); end
                    R_eq=R_eq_t;
            end
            % decoding bits
            dec_bits = real(R_eq)<0;
            % calculate BER
            ber_bpsk_ofdm_mrc_tv(k) = ber_bpsk_ofdm_mrc_tv(k)+sum(bits~=dec_bits);
        end
    end
    ber_bpsk_ofdm_mrc_tv = ber_bpsk_ofdm_mrc_tv./(N_test*N_bits);
    switch DM
        case 1
            semilogy(snr_db, ber_bpsk_ofdm_mrc_tv, 'b-o');
        case 2
            semilogy(snr_db, ber_bpsk_ofdm_mrc_tv, 'b:o');
        case 3
            semilogy(snr_db, ber_bpsk_ofdm_mrc_tv, 'b-.o');
    end
            
    set(gca,'FontName','TimesNewRoman');
    set(gca,'FontSize',12);
    grid on;
    xlabel('SNR [dB]');
    ylabel('BER');
    title({'BPSK modulated OFDM system in freq sel time variant'; strcat('fading 1x', num2str(N), ' channel with MRC')});
elseif mode == 6
% -------------------------------------------------------------------------------------
% Example 6 -> PAPR evaluation for OFDM
% -------------------------------------------------------------------------------------
    disp('Running example for PAPR evaluation');
    % time axis
    t = linspace(0, T, 10240);
    max_p=[];
    mean_p=[];
    % calculate continuous time OFDM signal N_test times
    for n=1:N_test
        if rem(n,disp_step)==0, str=sprintf('n=%5d, PAPR=%.1f',n,mean(max_p./mean_p)); disp(str), end
        % random bits
        bits = rand(1, N_bits)>.5;
        % BPSK modulation ("0"->+1, "1"->-1)
        s = 1-2*bits;
        s_ofdm = zeros(size(t));
        for k=0:N_fft-1
            s_ofdm = s_ofdm + s(k+1).*exp(j*2*pi*(k-N_fft/2)*Delta_f*t);
        end
        % calculate OFDM samples with IDFT
        s_shift = fftshift(s);
        s_ofdm_samples = N_fft*ifft(s_shift);
        t_samples = T_s*(0:N_fft-1);
        max_p(n)=max((abs(s_ofdm_samples)).^2);
        mean_p(n)=mean((abs(s_ofdm_samples)).^2);
    end
    [yc,yp]=ecdf2(xc,max_p./mean_p);
    semilogy(xc,1-yc,'b-.');
    set(gca,'FontName','TimesNewRoman');
    set(gca,'FontSize',12);
    xlabel('PAPR_0');
    ylabel('p(PAPR>PAPR_0)');
    hold on
    plot (xc, 1-(1-exp(-xc)).^N_fft,'r-.');  
    title('Peak to Average Power Ratio (PAPR) OFDM approssimazione gaussiana (rosso), simulazione (blu)');
elseif mode == 7
% -------------------------------------------------------------------------------------
% Example 7 -> performance of a BPSK modulated SC-FDM system in a freq selective channel
% -------------------------------------------------------------------------------------
    disp('Running example for performance of a BPSK modulated OFDM system in a freq selective channel');
    ber_bpsk_ofdm = zeros(size(snr_db));
    for k=1:length(snr)
        for n=1:N_test
            % random bits
            bits = rand(1, M_bits)>.5;
            % BPSK modulation ("0"->+1, "1"->-1)
            S0 = 1-2*bits;
            S=fft(S0);
            S=[S, zeros(1,(N_fft-M_dft))];
            % restore correct energy
            S = S./sqrt(N_fft);
            % OFDM modulation (shift does not matter in discrete time simulations)
            s_ofdm = N_fft*ifft(S);
            % add CP
            s_cp = [s_ofdm(end-N_cp+1:end), s_ofdm];
            % Rayleigh uncorrelated fading in each channel tap
            h = sqrt(.5).*(randn(1, L_ch) + j*randn(1, L_ch));
            % normalize channel impulse response
            h = h./sqrt(L_ch);
            % AWGN noise
            z = sqrt(.5/snr(k)).*(randn(1, length(s_cp)+L_ch-1) + j*randn(1, length(s_cp)+L_ch-1));
            % simultion of noisy reception through frequency selective fading channel
            r = conv(h, s_cp) + z;
            % discard first and last samples (including CP)
            r = r(N_cp+1:N_cp+N_fft);
            % OFDM demodulation
            R = fft(r);
            % calculate channel in frequency domain (this simulates perfect CIR estimation)
            H = N_fft*fft(h, N_fft);
            % ZF equalization
            R_eq = R./H;
            % decoding bits
            R_eq1=R_eq(1:M_dft);
            r_eq1=ifft(R_eq1);
            dec_bits = real(r_eq1)<0;
            % calculate BER
            ber_bpsk_ofdm(k) = ber_bpsk_ofdm(k)+sum(bits~=dec_bits);
        end
    end
    ber_bpsk_ofdm = ber_bpsk_ofdm./(N_test*M_bits);
    semilogy(snr_db, ber_bpsk_ofdm, 'k:o');
    set(gca,'FontName','TimesNewRoman');
    set(gca,'FontSize',12);
    grid on;
    xlabel('SNR [dB]');
    ylabel('BER');
    title('BPSK modulated OFDM system freq sel channel');    
elseif mode == 8
% -------------------------------------------------------------------------------------
% Example 8 -> PAPR evaluation for SC-FDM
% -------------------------------------------------------------------------------------
    disp('Running example for PAPR evaluation for SC-FDM');
    % time axis
    t = linspace(0, T, 10240);
    max_p=[];
    mean_p=[];
    % calculate continuous time OFDM signal N_test times
    for n=1:N_test
        if rem(n,disp_step)==0, str=sprintf('n=%5d, PAPR=%.1f',n,mean(max_p./mean_p)); disp(str), end
        % random bits
        bits = rand(1, M_bits)>.5;
        % BPSK modulation ("0"->+1, "1"->-1)
        S0 = 1-2*bits;
        S=fft(S0);
        s=[S, zeros(1,N_fft-M_dft)];     
        s_ofdm = zeros(size(t));
        for k=0:N_fft-1
            s_ofdm = s_ofdm + s(k+1).*exp(j*2*pi*(k-N_fft/2)*Delta_f*t);
        end
        % calculate OFDM samples with IDFT
        s_shift = fftshift(s);
        s_ofdm_samples = N_fft*ifft(s_shift);
        t_samples = T_s*(0:N_fft-1);
        max_p(n)=max((abs(s_ofdm_samples)).^2);
        mean_p(n)=mean((abs(s_ofdm_samples)).^2);
    end
    [yc,yp]=ecdf2(xc,max_p./mean_p);
    semilogy(xc,1-yc,'k-.');
    set(gca,'FontName','TimesNewRoman');
    set(gca,'FontSize',12);
    xlabel('PAPR_0');
    ylabel('p(PAPR>PAPR_0)');
    title('Peak to Average Power Ratio (PAPR) SC-FDM (black)')
end