Compute the theoretical BER for AWGN channel and various constellation size

Based on the previous skeleton, we can now compute an iterative testbench to compute the Bit Error Rate for various constellation size, and compare the simulation with the theory. As a gentle reminder, the theoretical bit error rate can be approximated as

$\mathrm{BER} = \frac{ 4 \left( 1 - \frac{1}{\sqrt{M}} \right) }{ \log_2(M)} \times Q( \sqrt{ \frac{6 \log_2(M)}{2(M-1)} \frac{Eb}{N_0}}$

First of all let's call the modules

using DigitalComm 
using PGFPlotsX

We define first the main monte-carlo function that compute an elementary Tx-Rx link, and returns the number of error and number of bit computed (to be accumulated)

function monteCarlo(snr,mcs,nbSymb) 
    # Number of bits  
    nbBits    = nbSymb * Int(log2(mcs));
    # --- Binary sequence generation 
    bitSeq    = genBitSequence(nbBits);
    # --- QPSK mapping
    qamSeq    = bitMappingQAM(mcs,bitSeq);
    # ---------------------------------------------------- 
    # --- Channel  
    # ---------------------------------------------------- 
    #  --- AWGN
    # Theoretical power is 1 (normalized constellation)
    qamNoise,  = addNoise(qamSeq,snr,1);
    # ----------------------------------------------------
    # --- Rx Stage: SRRC
    # ----------------------------------------------------
    # --- Binary demapper
    bitDec  = bitDemappingQAM(mcs,qamNoise);
    # --- Error counter 
    nbE     = sum(xor.(bitDec,bitSeq));
    # --- Return Error and bits 
    return (nbE,nbBits);
end

A function to plot the BER versus the SNR, for different mcs and compare to theory

function doPlot(snrVect,ber,qamVect)
    a = 0;
    @pgf a = Axis({
        ymode     = "log",
        height      ="3in",
        width       ="4in",
        grid,
        xlabel      = "SNR [dB]",
        ylabel      = "Bit Error Rate ",
        ymax       = 1,
        ymin       = 10.0^(-5),
        title       = "AWGN BER for QAM",
        legend_style="{at={(0,0)},anchor=south west,legend cell align=left,align=left,draw=white!15!black}"
        },
        Plot({color="red",mark="square*"},Table([snrVect,ber[1,:]])),
        LegendEntry("QPSK"),
        Plot({color="green",mark="*"},Table([snrVect,ber[2,:]])),
        LegendEntry("16-QAM"),

        Plot({color="purple",mark="triangle*"},Table([snrVect,ber[3,:]])),
        LegendEntry("64-QAM"),
        Plot({color="blue",mark="diamond*"},Table([snrVect,ber[4,:]])),
        LegendEntry("256-QAM"),
    );
    # ---  Adding theoretical curve
    snrLin  = (10.0).^(snrVect/10)
    for qamScheme = qamVect
        ebNo    = snrLin / log2(qamScheme);
        # This approximation is only valid for high SNR (one symbol error is converted to one bit error with Gray coding).
        berTheo   = 4 * ( 1 - 1 / sqrt(qamScheme)) / log2(qamScheme) * qFunc.(sqrt.( 2*ebNo * 3 * log2(qamScheme) / (2*(qamScheme-1)  )));
        @pgf push!(a,Plot({color="black"},Table([snrVect,berTheo])));
    end
    display(a);
end

Then, the main routine to compute the BER for a given number of iterations and a range of SNR

function main() 
# --- Parameters 
nbIt            = 10000;          # Number of iterations  
nbSymb          = 1024;           # Number of symbols per iterations 
mcs             = [4,16,64,256];  # Constellation size 
snrRange        = (-1:26);        # SNR, expressed in dB 
# --- Init performance metrics
nbSNR           = length(snrRange);
ber             = zeros(Float64,length(mcs),nbSNR);
for iMcs = 1 : 1  : length(mcs)
    for iSNR = 1 : 1 : nbSNR
        # --- Create BER counters 
        nbE   = 0;
        nbB   = 0;
        for iN = 1 : 1 : nbIt
            # --- Elementary MC call 
            # Corresponds to a given SNR and a given iteration 
            # As we are ergodic in AWGN, it is only nbSymb*nbIt that matters for BER computation
            (a,b) = monteCarlo(snrRange[iSNR],mcs[iMcs],nbSymb); 
            # --- Update counters 
            nbE += a;     # Increment errors 
            nbB += b;     # Increment bit counters 
        end 
        ber[iMcs,iSNR] = nbE / nbB; 
    end
end 
# --- Plotting routine
doPlot(snrRange,ber,mcs);
end

The output plot is the following, showing adequacy between theory and practise for high SNR (the theoretical curve is under the assumption that one symbol error leads to one erroneous bit (gray coding) which is true only with intermediate noise levels).