SNR Calculations in NeoRadium

This notebook demonstrates how NeoRadium calculates the noise standard deviation from a given SNR value in both the time and frequency domains.

import numpy as np
import scipy.io
import time

from neoradium import Carrier, PDSCH, CdlChannel, AntennaPanel, LdpcEncoder, Grid, random, Waveform
from neoradium.utils import toDb, toLinear, getNmse

First let’s try a simple case without a channel model:

# NOTE: You can find a similar Matlab example at:
# https://www.mathworks.com/help/5g/ug/snr-definition-used-in-link-simulations.html
snrDb = 0
snr = toLinear(snrDb)   # SNR in linear scale

carrier = Carrier(numRbs=52, spacing=30)
bwp = carrier.curBwp
nr, nt = 2, 2  # receiver and transmitter antennas
print(bwp)

pdsch = PDSCH(bwp, interleavingBundleSize=0, numLayers=1, modulation='16QAM', nID=carrier.cellId)
pdsch.setDMRS(prgSize=0, configType=2, additionalPos=2)

codeRate = 490/1024
ldpcEncoder = LdpcEncoder(baseGraphNo=1, modulation=pdsch.modems[0].modulation,
                          txLayers=pdsch.numLayers, targetRate=codeRate)



random.setSeed(123)                          # Make the results reproducible.
grid = pdsch.getGrid()                       # Create a resource grid populated with DMRS
txBlockSize = pdsch.getTxBlockSize(codeRate) # Calculate the transport block size based on 3GPP TS 38.214
txBlock = random.bits(txBlockSize[0])        # Create random binary data

numBits = pdsch.getBitSizes(grid)            # Actual number of bits available in the resource grid

# Now perform the segmentation, rate-matching, and encoding in one call:
rateMatchedCodeBlocks = ldpcEncoder.getRateMatchedCodeBlocks(txBlock, numBits[0])

# Now populate the resource grid with coded data. This includes QAM modulation and resource mapping.
pdsch.populateGrid(grid, rateMatchedCodeBlocks)

# 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
precoder = np.ones((nt, pdsch.numLayers))/np.sqrt(pdsch.numLayers)      # Get the precoder matrix from the PDSCH object
precodedGrid = grid.precode(precoder)                                   # Perform the precoding

print(f"Shape of Resource Grid:         {grid.shape}")
print(f"Shape of Precoded Resource Grid:{precodedGrid.shape}")

txWaveform = precodedGrid.ofdmModulate()
print(f"Shape of txWaveform:            {txWaveform.shape}")

# Applying channel to the transmitted signal
rxWaveform = Waveform(txWaveform.waveform/np.sqrt(nr))
print(f"Shape of rxWaveform:            {rxWaveform.shape}")
rxGrid = rxWaveform.ofdmDemodulate(bwp)
print(f"Shape of rxGrid:                {rxGrid.shape}")

Bandwidth Part Properties:
  Resource Blocks:    52 RBs starting at 0 (624 subcarriers)
  Subcarrier Spacing: 30 kHz
  CP Type:            normal
  Bandwidth:          18.72 MHz
  symbolsPerSlot:     14
  slotsPerSubFrame:   2
  nFFT:               1024
  frameNo:            0
  slotNo:             0

Shape of Resource Grid:         (1, 14, 624)
Shape of Precoded Resource Grid:(2, 14, 624)
Shape of txWaveform:            (2, 15360)
Shape of rxWaveform:            (2, 15360)
Shape of rxGrid:                (2, 14, 624)

The standard deviation of the noise in time domain is:

\[\sigma_T = \sqrt{\frac {N_{FFT} \sigma_x^2} {K.{SNR}}}\]

where \(\sigma_x^2\) is the variance of the time-domain signal (the Waveform object), \(K\) is the number of subcarriers in the current bandwidth part, and \(N_{FFT}\) is the FFT size (bwp.nFFT). This value can be used to generate Gaussian noise applied to a Waveform object:

noiseStdTime = rxWaveform.getNoiseStd(snr, bwp)
print(f"Noise STD (Time): {noiseStdTime}")

Noise STD (Time): 0.0220441589451537

The standard deviation of the noise in Frequency domain is:

\[\sigma_F = \sqrt{\frac {\sigma_X^2} {SNR}}\]

where \(\sigma_X^2\) is the variance of the frequency-domain signal (the Grid object). This can be used to generate Gaussian noise applied to a Grid object:

noiseStdFreq = rxGrid.getNoiseStd(snr)
print(f"Noise STD (Freq): {noiseStdFreq}")

Noise STD (Freq): 0.705375297931587

Verifying Equation:

\[\sigma_T = \sqrt{{\frac 1 {N_{FFT}}}} \sigma_F\]

print(f"Noise STD (Freq)/sqrt({bwp.nFFT}): {noiseStdFreq/np.sqrt(bwp.nFFT)} = Noise STD (Time)")

Noise STD (Freq)/sqrt(1024): 0.022042978060362095 = Noise STD (Time)

# Apply noise in time domain and measure noise in frequency domain:
noisyRxWaveform = rxWaveform.addNoise(noiseStd=noiseStdTime)
noisyRxGrid = noisyRxWaveform.ofdmDemodulate(bwp)
noiseGrid = noisyRxGrid.grid - rxGrid.grid
print(f"Noise Grid STD: {noiseGrid.std()}")
print(f"Noise STD (Freq): {noiseStdFreq}")

Noise Grid STD: 0.7020339000063331
Noise STD (Freq): 0.705375297931587

Matlab’s Assumption

Matlab assumes \(\sigma_{X}^2 = \frac{1}{N_r}\), where \(N_r\) is the number of receiver antennas. With that assumption, we have:

\(\sigma_F = \frac{1}{\sqrt{N_r \cdot SNR}}\) and \(\sigma_T = \frac{1}{\sqrt{N_r \cdot N_{FFT} \cdot SNR}}\)

Note that in this case, these values are the same since we used \(\sigma_{X}^2 = \frac{1}{N_r}\). This is not always true in practice.

noiseStdTime1 = 1/np.sqrt(nr*bwp.nFFT*snr)
noiseStdFreq1 = 1/np.sqrt(nr*snr)
print(f"Noise STD (Time Domain - RxPower=1/Nr): {noiseStdTime1}")
print(f"Noise STD (Freq Domain - RxPower=1/Nr): {noiseStdFreq1}")

Noise STD (Time Domain - RxPower=1/Nr): 0.022097086912079608
Noise STD (Freq Domain - RxPower=1/Nr): 0.7071067811865475

Now we try the same steps with a CDL channel model:

# Now repeat the same process using a CDL channel model with two layers:
snrDb = 0
snr = toLinear(snrDb)

carrier = Carrier(numRbs=52, spacing=30)
bwp = carrier.curBwp

pdsch = PDSCH(bwp, interleavingBundleSize=0, numLayers=2, modulation='16QAM', nID=carrier.cellId)
pdsch.setDMRS(prgSize=0, configType=2, additionalPos=2)

codeRate = 490/1024
ldpcEncoder = LdpcEncoder(baseGraphNo=1, modulation=pdsch.modems[0].modulation,
                          txLayers=pdsch.numLayers, targetRate=codeRate)

random.setSeed(123)                          # Making the results reproducible.
grid = pdsch.getGrid()                       # Create a resource grid already populated with DMRS
txBlockSize = pdsch.getTxBlockSize(codeRate) # Calculate the Transport Block Size based on 3GPP TS 38.214
txBlock = random.bits(txBlockSize[0])        # Create random binary data

numBits = pdsch.getBitSizes(grid)            # Actual number of bits available in the resource grid

# Now perform the segmentation, rate-matching, and encoding in one call:
rateMatchedCodeBlocks = ldpcEncoder.getRateMatchedCodeBlocks(txBlock, numBits[0])

# Now populate the resource grid with coded data. This includes QAM modulation and resource mapping.
pdsch.populateGrid(grid, rateMatchedCodeBlocks)

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

# Creating a CdlChannel object:
f0 = 4e9  # Carrier Frequency
channel = CdlChannel(bwp, 'C', delaySpread=300, carrierFreq=f0, dopplerShift=5,
                     txAntenna = AntennaPanel([1,4], polarization="|"),  # 4 TX antennas
                     rxAntenna = AntennaPanel([1,2], polarization="|"),  # 2 RX antennas
                     normalizeGains = True, normalizeOutput=True)
nr,nt = channel.nrNt

# 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

print(f"Shape of Resource Grid:         {grid.shape}")
print(f"Shape of Precoded Resource Grid:{precodedGrid.shape}")

txWaveform = precodedGrid.ofdmModulate()
print(f"Shape of txWaveform:            {txWaveform.shape}")

rxWaveform = channel.applyToSignal(txWaveform)
print(f"Shape of rxWaveform:            {rxWaveform.shape}")

offset = channel.getTimingOffset()
syncedWaveform = rxWaveform.sync(offset)
print(f"Timing Offset:                  {offset}")

rxGrid = syncedWaveform.ofdmDemodulate(bwp)
print(f"Shape of rxGrid:                {rxGrid.shape}")

rxGridF = channel.applyToGrid(precodedGrid)
print(f"Shape of rxGridF:               {rxGridF.shape}")
print(f"NMSE(rxGrid,rxGridF):           {getNmse(rxGrid.grid, rxGridF.grid)}")

Shape of Resource Grid:         (2, 14, 624)
Shape of Precoded Resource Grid:(4, 14, 624)
Shape of txWaveform:            (4, 15360)
Shape of rxWaveform:            (2, 15360)
Timing Offset:                  13
Shape of rxGrid:                (2, 14, 624)
Shape of rxGridF:               (2, 14, 624)
NMSE(rxGrid,rxGridF):           6.457297490367846e-07

noiseStdTime = syncedWaveform.getNoiseStd(snr, bwp)
print(f"Noise STD (Time): {noiseStdTime}")

Noise STD (Time): 0.06558965831573078

noiseStdFreq = rxGrid.getNoiseStd(snr)
print(f"Noise STD (Freq): {noiseStdFreq}")

Noise STD (Freq): 2.0989851489946787

print(f"Noise STD (Freq)/sqrt({bwp.nFFT}): {noiseStdFreq/np.sqrt(bwp.nFFT)} = Noise STD (Time)")

Noise STD (Freq)/sqrt(1024): 0.06559328590608371 = Noise STD (Time)

# In this case, Matlab's assumption is not valid; these values do not match the calculated STD values above.
noiseStdTime1 = 1/np.sqrt(nr*bwp.nFFT*snr)
noiseStdFreq1 = 1/np.sqrt(nr*snr)
print(f"Noise STD (Time Domain - RxPower=1/Nr): {noiseStdTime1}")
print(f"Noise STD (Freq Domain - RxPower=1/Nr): {noiseStdFreq1}")

Noise STD (Time Domain - RxPower=1/Nr): 0.022097086912079608
Noise STD (Freq Domain - RxPower=1/Nr): 0.7071067811865475

# Apply noise in the time domain and measure the resulting noise in the frequency domain:
noisyRxWaveform = rxWaveform.addNoise(noiseStd=noiseStdTime)
noisySyncedWaveform = noisyRxWaveform.sync(offset)
noisyRxGrid = noisySyncedWaveform.ofdmDemodulate(bwp)
noiseGrid = noisyRxGrid.grid - rxGrid.grid
print(f"Noise Grid STD: {noiseGrid.std()}")
print(f"Noise STD (Freq): {noiseStdFreq}")

Noise Grid STD: 2.082620653246523
Noise STD (Freq): 2.0989851489946787