gv3mqam_type2_ofdm_cp_zf.m

4 16 64 QAM plots for simulated and theoretical results for OFDM system with addition of Cyclic pref

%  Title       : QAM Constellation plots for 4, 16, 64
%  Author      : goutam.ramamurthy.06<AT>student.lth.se
%  History     : 2009-01-30      created
%  References  :
%  Description :
%
%        1         2         3         4         5         6         7         
%23456789012345678901234567890123456789012345678901234567890123456789012345
%

close all;
clear all;
clc;

%%%%%%%%%% GLOBAL SYSTEM PARAMETERS 
%--Case-Specific Simulation Parameters
SIMU{1}.name = '4QAM';
SIMU{1}.M_ary = 4;
SIMU{1}.legend = ' Coherent 4QAM,QPSK Simulation';
SIMU{1}.linestyle = '-gd'; % linestyle: green Squares
SIMU{1}.title ='Ps(Sym Err) curve for 4QAM,QPSK modulation';


% SIMU{2}.name = '16QAM';
% SIMU{2}.M_ary = 16;
% SIMU{2}.legend = ' Coherent 16QAM Simulation';
% SIMU{2}.linestyle = 'rd-'; % linestyle: red diamonds
% SIMU{2}.title ='Ps(Sym Err) curve for 16QAM modulation';


% SIMU{3}.name = '64QAM';
% SIMU{3}.M_ary = 64;
% SIMU{3}.legend = ' Coherent 64QAM AWGN Simulation';
% SIMU{3}.linestyle = '-md'; % linestyle: maroon diamonds
% SIMU{3}.title ='Ps(Sym Err) curve for 64QAM modulation';


%--Common Parameters 
SNRdB = -2:2:24;     % Es/N0 range for simulation
bSNRdB = zeros(length(SIMU), length(SNRdB));

%--Channel Parameters 
Lch = 0;                % length of channel dispersion.
h_fixed = (1:1:Lch+1)/Lch;    % Fixed ch coeffs
CH = 'PURE_AWGN';       % PURE_AWGN; RAYLI; TID_DISP
ON_noise = 1;       % 1 = ON; 0 = OFF
ON_esti = 0;        % 1 = ON; 0 = OFF

%--OFDM Parameters 
N_OFDMsyms = 20;      % Choose the total # of symbols in 1RB to simulateTXRX
BW_LTE = 20e6;        % Total Available LTE S/M bandwidth is 20 MHz
BW_guard = BW_LTE/10; % 10% of Total available BW for OOB emissions 
delta_f = 15e3;       % LTE Sub-Carrier Spacing; specified by 3GPP
                      % also delta_f_low = 7.5e3us used for MBMS
L_CP = 4;             % Normal CP first symbol. L_CP=160, L_CP_2 = 144;(LTE)

%--Derived Parameters
N_c = floor((BW_LTE - BW_guard)/delta_f);   % 1200 for 20MHz BW
m = ceil(log2(N_c));    % radix value for N=2^m in N-IFFT
N_FFT = 2^m;            % N > N_c & for radix-2 IFFT, S.T. N is 2^m;


%~~~~~~~~~~~~~~

tic;
for iSIMU = 1:length(SIMU)     % iterating through Simulations
    M_ary = SIMU{iSIMU}.M_ary; % number of constellation points
    k = log2(M_ary);
    Es = (M_ary-1)*2/3;        % Avg. Ener per QAM sym, norm = sqrt(Es)
    d_minsq = 3*log2(M_ary)/(M_ary-1); % for M-QAM
    
    N_qsyms = N_c * N_OFDMsyms ;
    N_bits = N_qsyms * k;      % an integral multiple of k channels
    
    p = (1:sqrt(M_ary)/2);        % In rect QAMs, for M=2^n & n ~ {even I}
    a_MQAM = [-(2*p-1) (2*p-1)];  % alphabets
    
    for iSNRdB = 1:length(SNRdB) % iterating through SNRs
        fprintf('\nEs/N0: %d dB\n', SNRdB(iSNRdB));
        SNR = 10^(SNRdB(iSNRdB)/10); % Convert to non-dB value
        bSNRdB(iSIMU,iSNRdB) = SNRdB(iSNRdB) - 10*log10(k);
        
        %%%%%%%%%% RANDOM SOURCE
        %--Generate a column of random binary bits
        x_txbits = (1-sign(randn(N_bits,1)))/2; % ~Bernoulli source(p=1/2)
        %%%RANDOM SOURCE-ED
        
        %%%%%%%%%% QAM MAPPER
        %--using matlab's modem object generate QAM syms and normalize
        qam_mod = modem.qammod(M_ary);    % create a modulator object
        qam_mod.symbolorder = 'gray' ;  % Enable/disable gray coding
        qam_mod.inputtype = 'bit';
        qam_mod.disp;                   % Print the properties of modulator
        q_txsyms = modulate(qam_mod, x_txbits);  % QAM modulate
        q_txsyms_norm = q_txsyms/sqrt(Es);
        %%% QAM MAPPER -ED

        %%%%%%%%%% OFDM MOD
        O_train = [];                        % pre-declaration
        for i = 1:N_OFDMsyms                 % iterate over N OFDMsyms
            txsyms_ext = zeros(N_FFT, 1);    % pad zeros to unused FFT bins
            txsyms_ext(1:N_c/2) = q_txsyms_norm((i-1)*N_c+1 : ...
                                              (i-1)*N_c + N_c/2).';
            txsyms_ext(N_FFT-(N_c/2-1):N_FFT) = q_txsyms_norm( ...
                                                (i-1)*N_c+(N_c/2+1):...
                                                (i-1)*N_c+N_c).';
            O_sym = sqrt(N_FFT).*ifft(txsyms_ext, N_FFT);
            %%%%%%%%%% ADD CP  
            O_train = [O_train;O_sym(N_FFT-L_CP+1:end);O_sym];
            %%%ADD CP -ED
        end
        %%%OFDM MOD -ED
        
        
        %++++++++++
        %%%%%%%%%% THE CHANNEL TO CONVOLVE
        noise_len = size(O_train);
        re = randn(noise_len);
        im = randn(noise_len);
        awgn = 1/sqrt(2)*(re + j*im); % WGN with 0dB variance
        
        if strcmp(CH, 'RAYLI_FAD') ~= 0
            h = (randn+j*randn)/sqrt(2); % Over-ride the h
            rxdata = filter(h,1,O_train) + ...
                           ON_noise*(1/sqrt(SNR))*awgn; % AWGN;
        elseif strcmp(CH, 'TID_DISP') ~= 0
            h = h_fixed;
            rxdata = filter(h,1,O_train)+ ...
                           ON_noise*(1/sqrt(SNR))*awgn; % AWGN;
        else
            rxdata = O_train + ON_noise*(1/sqrt(SNR))*awgn; % AWGN
        end
        H=fft(h,N_FFT);
        %%%CHANNEL CONVOLUT-ED 
        %++++++++++

        
        %%%%%%%%%% OFDM DEMODulate and equalize
        %%GAVI : take care of L_CP before estimation
        if(ON_esti)
            if ON_pilot  % Collect Equ coeffs from the rx_data
                
            else         % Channel is KNOWN, h
                Equ = (1./H)/sqrt(abs(sum(H).^2)); 
                Equ = Equ.';         % ZF equalizer coeffecients
            end
        else 
            Equ = ones(N_FFT,1);
        end
        
        O_rxtrain = [];          % pre-declaration
        rxqsyms = [];            % pre-declaration
        for i = 1:N_OFDMsyms     % iterate over N OFDMsyms
            offset = (i-1)*N_FFT + (i*L_CP);
            %%%%%%%%%% REM CP
            if(ON_esti)
                O_rxsym = fft(rxdata(offset+1:offset+N_FFT),N_FFT);
            else
                O_rxsym = (1/sqrt(N_FFT)).*fft(rxdata(offset+1:offset+N_FFT),N_FFT);
            end
            %%%REM CP -ED
            O_rxsym = Equ.*O_rxsym;
            O_rxtrain = [O_rxtrain; O_rxsym];  % DEL: for debug
            rxqsyms = [rxqsyms; O_rxsym(1:N_c/2); 
                       O_rxsym((N_FFT-N_c/2+1):N_FFT)];
        end
        %%%OFDM DEMOD -ED        
        
        %%%%%%%%%% QAM DEMAPPER
        %--using matlab's modem object
        qam_rxsyms_denorm = rxqsyms * sqrt(Es);
        
            %%%%%%%%%% RECT QAM MAPPING
            %--Segregate Re-Im
            q_rxsyms_re = real(qam_rxsyms_denorm);
            q_rxsyms_im = imag(qam_rxsyms_denorm);
            %--Mapping to the nearest alphabet
            rx_re = floor(q_rxsyms_re/2)*2+1;
            rx_im = floor(q_rxsyms_im/2)*2+1;
            %--mapping boundary conditions for RE and IM
            rx_re(find(max(a_MQAM)<rx_re)) = max(a_MQAM);
            rx_re(find(min(a_MQAM)>rx_re)) = min(a_MQAM);
            rx_im(find(max(a_MQAM)<rx_im)) = max(a_MQAM);
            rx_im(find(min(a_MQAM)>rx_im)) = min(a_MQAM);
            q_rxsyms = rx_re+j*rx_im;
            %%%RECT MAPPED -ED

        qam_demod = modem.qamdemod(qam_mod);  % copy props from mod object
        qam_demod.disp;                     % Prints properties of DEMOD
        x_rxbits = demodulate(qam_demod, qam_rxsyms_denorm); % QAM modulate
        %%%QAM DEMAPPER -ED 
       
        %%%%%%%%%% PERF. MEASUREMENTS
        %--Simulated SER
        N_symerrs = nnz(q_rxsyms(:)-q_txsyms(:));
        SER(iSIMU,iSNRdB) = N_symerrs/N_qsyms;
        %--theoretical SER
        EbN0 = SNR/k;           % EbN0 = EsN0/k
        Qfn = qfunc(sqrt(d_minsq * EbN0));
        SER_theo(iSIMU,iSNRdB) = (4*(1-(1/sqrt(M_ary)))*Qfn) - ...
                                  (4*(1-(1/sqrt(M_ary)))^2*(Qfn^2));
        %--Simulated BER
        N_biterrs = nnz(x_rxbits(:)-x_txbits(:));
        BER(iSIMU,iSNRdB) = N_biterrs/N_bits;
        %--theoretical BER
        BER_theo(iSIMU, iSNRdB) = (1/k)*SER_theo(iSIMU, iSNRdB);
        %%%PERF. MEASURE -ED
        %-- Scatter plots
        figure(iSIMU*10);
        plot(qam_rxsyms_denorm,'.k');
        title(SIMU{iSIMU}.legend);
        axis square;
        axis equal;
        grid on;
    end
end
toc;

figure(200)
plot(abs(txsyms_ext));
grid on

%%%%%%%%%% SER PLOTS
figure(1)
for iSIMU = 1:length(SIMU)     % iterating through Simulations
    semilogy(SNRdB, SER_theo(iSIMU, :),'b:','LineWidth',2);
    hold on
    semilogy(SNRdB, SER(iSIMU, :), SIMU{iSIMU}.linestyle, 'Linewidth',1);
    hold on
end
hold off
axis([-2 26 10^-7 1]);
grid on
legend('theoretical', 'Simulated');
xlabel('Es/No, dB');
ylabel('SER(P_s)');
title('Ps versus E_s/N_0');

% figure(2)
% for iSIMU = 1:length(SIMU)     % iterating through Simulations
%     semilogy(SNRdB, BER_theo(iSIMU, :),'b:','LineWidth',2);
%     hold on
%     semilogy(SNRdB, BER(iSIMU, :), SIMU{iSIMU}.linestyle, 'Linewidth',1);
%     hold on
% end
% hold off
% axis([-2 26 10^-7 1]);
% grid on
% legend('theoretical', 'Simulated');
% xlabel('Es/No, dB');
% ylabel('BER(P_b)');
% title('Pb versus E_s/N_0');
% 
% figure(3)
% for iSIMU = 1:length(SIMU)     % iterating through Simulations
%     semilogy(bSNRdB(iSIMU,:), BER_theo(iSIMU, :),'b:','LineWidth',2);
%     hold on
%     semilogy(bSNRdB(iSIMU,:), BER(iSIMU, :), SIMU{iSIMU}.linestyle, 'Linewidth',1);
%     hold on
% end
% hold off
% axis([-2 26 10^-7 1]);
% grid on
% legend('theoretical', 'Simulated');
% xlabel('Eb/No, dB');
% ylabel('BER(P_b)');
% title('Pb versus E_b/N_0');