Calculating BER/BLER for PDSCH Communication with LDPC

This notebook shows how to calculate the BLock Error Rate (BLER) of PDSCH communication with LDPC channel coding.

[1]:
import numpy as np
import time
import matplotlib.pyplot as plt

from neoradium import Carrier, PDSCH, CdlChannel, AntennaPanel, LdpcEncoder, Grid, random, SnrScheduler

[2]:
numSlots = 200
snrScheduler = SnrScheduler(6,0.2)         # Start at 6 dB, use increments of 0.2 dB
freqDomain = True                          # Set to True to apply channel in frequency domain

modulation = "16QAM"
carrier = Carrier(numRbs=51, spacing=30)    # Create a carrier with 51 RBs and 30KHz subcarrier spacing
bwp = carrier.curBwp                        # The only bandwidth part in the carrier

# Create a PDSCH object
pdsch = PDSCH(bwp, interleavingBundleSize=0, numLayers=2, nID=carrier.cellId, modulation=modulation)
pdsch.setDMRS(prgSize=0, configType=2, additionalPos=2)     # Specify the DMRS configuration

# Create the LDPC encoder
codeRate = 490/1024
ldpcEncoder = LdpcEncoder(baseGraphNo=1, modulation=pdsch.modems[0].modulation,
                          txLayers=pdsch.numLayers, targetRate=codeRate)
ldpcDecoder = ldpcEncoder.getDecoder()

results = {}
for chanEstMethod in ["Perfect", "LS"]:                 # Two channel estimation methods
    results[chanEstMethod] = {}
    print("\nSimulating end-to-end for \"%s\", with \"%s\" channel estimation, in %s domain."%
          (modulation, chanEstMethod, "frequency" if freqDomain else "time"))
    print("SNR(dB)    Total Bits  Bit Errors  BER(%)  Total Blocks  Block Errors  BLER(%)  time(Sec.)")
    print("---------  ----------  ----------  ------  ------------  ------------  -------  ----------")
    snrScheduler.reset()
    for snrDb in snrScheduler:
        random.setSeed(123)                             # Making the results reproducible for each SNR
        t0 = time.time()                                # Start time for each SNR
        carrier.slotNo = 0

        # Creating a CdlChannel object:
        channel = CdlChannel(bwp, 'C', delaySpread=300, carrierFreq=4e9, dopplerShift=5,
                             txAntenna = AntennaPanel([2,4], polarization="x"),  # 16 TX antenna
                             rxAntenna = AntennaPanel([1,2], polarization="x"))  # 4 RX antenna

        blockErrors = 0
        totalBlocks = 0
        bitErrors = 0
        totalBits = 0

        for slotNo in range(numSlots):
            grid = pdsch.getGrid()                       # Create a resource grid already populated with DMRS
            txBlockSize = pdsch.getTxBlockSize(codeRate) # Calculate the Transport Block Size
            txBlock = random.bits(txBlockSize[0])        # Create random binary data
            numBits = pdsch.getBitSizes(grid)            # Actual number of bits available in the resource grid

            # Perform the segmentation, rate-matching, and encoding
            rateMatchedCodeBlocks = ldpcEncoder.getRateMatchedCodeBlocks(txBlock, numBits[0])

            pdsch.populateGrid(grid, rateMatchedCodeBlocks)         # Map/modulate the data to the resource grid

            # Store the indexes of the PDSCH data in pdschIndexes to be used later.
            pdschIndexes = pdsch.getReIndexes(grid, "PDSCH")

            # Getting the Precoding Matrix, and precoding the resource grid
            channelMatrix = channel.getChannelMatrix()              # Get the channel matrix
            precoder = pdsch.getPrecodingMatrix(channelMatrix)      # Get the precoder matrix from the PDSCH object
            precodedGrid = grid.precode(precoder)                   # Perform the precoding

            if freqDomain:
                rxGrid = precodedGrid.applyChannel(channelMatrix)   # Apply the channel in frequency domain
                rxGrid = rxGrid.addNoise(snrDb=snrDb, useRxPower=True)  # Add noise
            else:
                txWaveform = precodedGrid.ofdmModulate()            # OFDM Modulation
                maxDelay = channel.getMaxDelay()                    # Get the max. channel delay
                txWaveform = txWaveform.pad(maxDelay)               # Pad with zeros
                rxWaveform = channel.applyToSignal(txWaveform)      # Apply channel in time domain
                noisyRxWaveform = rxWaveform.addNoise(snrDb=snrDb, bwp=bwp, useRxPower=True)    # Add noise
                offset = channel.getTimingOffset()                  # Get timing info for synchronization
                syncedWaveform = noisyRxWaveform.sync(offset)       # Synchronization
                rxGrid = syncedWaveform.ofdmDemodulate(bwp)         # OFDM demodulation


            if chanEstMethod == "Perfect":                          # Perfect Channel Estimation
                estChannelMatrix = channelMatrix @ precoder[None,...]
            else:                                                   # LS + Interpolation Channel Estimation
                estChannelMatrix, noiseEst = rxGrid.estimateChannelLS(pdsch.dmrs, polarInt=False,
                                                                      kernel='linear')

            eqGrid, llrScales = rxGrid.equalize(estChannelMatrix)            # Equalization
            llrs = pdsch.getLLRsFromGrid(eqGrid, pdschIndexes, llrScales)    # Demodulation (to LLRs)
            rxCodedBlocks = ldpcDecoder.recoverRate(llrs[0], txBlockSize[0]) # Recovering Rate
            decodedBlocks = ldpcDecoder.decode(rxCodedBlocks, numIter=20)    # LDPC Decoding
            decodedTxBlockWithCRC, crcMatch = ldpcDecoder.checkCrcAndMerge(decodedBlocks) # Merge blocks
            decodedTxBlock = decodedTxBlockWithCRC[:-24]                     # remove transport block CRC
            blockErrors += len(crcMatch)-sum(crcMatch)                       # Number of Block errors
            bitErrors += np.abs(decodedTxBlock-txBlock).sum()                # Number of bit errors
            totalBlocks += len(crcMatch)
            totalBits += len(txBlock)

            ber = bitErrors*100/totalBits
            bler = blockErrors*100/totalBlocks
            print("\r  %5.1f    %8d    %8d    %6.2f    %8d      %8d     %6.2f     %6.2f"
                  %(snrDb, totalBits, bitErrors, ber, totalBlocks, blockErrors, bler, time.time()-t0), end='')

            channel.goNext()

        dt = time.time()-t0
        snrScheduler.setData(bler, ber)
        print("")
    results[chanEstMethod] = snrScheduler.getSnrsAndData()

Simulating end-to-end for "16QAM", with "Perfect" channel estimation, in frequency domain.
SNR(dB)    Total Bits  Bit Errors  BER(%)  Total Blocks  Block Errors  BLER(%)  time(Sec.)
---------  ----------  ----------  ------  ------------  ------------  -------  ----------
    6.0     6043200           0      0.00         800             0       0.00     110.32
    5.8     6043200         349      0.01         800             2       0.25     115.66
    5.6     6043200       23641      0.39         800           124      15.50     119.87
    5.4     6043200      219226      3.63         800           544      68.00     122.26
    5.2     6043200      521096      8.62         800           781      97.62     122.39
    5.0     6043200      724626     11.99         800           799      99.88     122.87
    4.8     6043200      824174     13.64         800           800     100.00     122.43
    4.6     6043200      896263     14.83         800           800     100.00     120.15
    6.2     6043200           0      0.00         800             0       0.00     119.54

Simulating end-to-end for "16QAM", with "LS" channel estimation, in frequency domain.
SNR(dB)    Total Bits  Bit Errors  BER(%)  Total Blocks  Block Errors  BLER(%)  time(Sec.)
---------  ----------  ----------  ------  ------------  ------------  -------  ----------
    6.0     6043200     1153377     19.09         800           800     100.00     120.07
    6.2     6043200     1113458     18.42         800           800     100.00     122.28
    6.6     6043200     1020319     16.88         800           800     100.00     119.69
    7.0     6043200      837213     13.85         800           782      97.75     119.66
    6.8     6043200      950361     15.73         800           800     100.00     120.07
    7.2     6043200      673459     11.14         800           678      84.75     120.70
    7.4     6043200      554758      9.18         800           488      61.00     119.48
    7.6     6043200      508271      8.41         800           407      50.88     121.37
    7.8     6043200      484415      8.02         800           400      50.00     120.86
    8.0     6043200      455392      7.54         800           400      50.00     119.83
    8.2     6043200      417771      6.91         800           400      50.00     118.66
    8.4     6043200      364735      6.04         800           395      49.38     121.64
    8.6     6043200      268860      4.45         800           363      45.38     123.30
    8.8     6043200      160604      2.66         800           262      32.75     122.01
    9.0     6043200       72932      1.21         800           160      20.00     122.01
    9.2     6043200       19074      0.32         800            51       6.38     121.36
    9.4     6043200        1319      0.02         800            11       1.38     121.92
    9.6     6043200           0      0.00         800             0       0.00     120.86
    9.8     6043200           0      0.00         800             0       0.00     119.47
[3]:
# Compare the results in a couple of plots:
logGraph = False

# Bit Error Rate:
for i,chanEstMethod in enumerate(['Perfect', 'LS']):
    plt.plot(results[chanEstMethod][0], results[chanEstMethod][2], label=chanEstMethod)
plt.legend()
plt.title("Bit Error Rate for different methods of Channel Estimation.");
plt.grid()
plt.xlabel("SNR (dB)")
# plt.xticks(results[chanEstMethod][0])
plt.ylabel("BER (%)")
if logGraph: plt.yscale('log')
plt.show()

# Block Error rate
for i,chanEstMethod in enumerate(['Perfect', 'LS']):
    plt.plot(results[chanEstMethod][0], results[chanEstMethod][1], label=chanEstMethod)
plt.legend()
plt.title("Block Error Rate for different methods of Channel Estimation.");
plt.grid()
plt.xlabel("SNR (dB)")
# plt.xticks(results[chanEstMethod][0])
plt.ylabel("BLER (%)")
if logGraph: plt.yscale('log')
plt.show()

../../../../_images/source_Playground_Notebooks_PDSCH_PDSCH-BLER_3_0.png
../../../../_images/source_Playground_Notebooks_PDSCH_PDSCH-BLER_3_1.png
[ ]: