Run in Google Colab

View on GitHub

Download notebook

Part 3: Advanced Link-level Simulations

This tutorial will guide you through Sionna, from its basic principles to the implementation of a point-to-point link with a 5G NR compliant code and a 3GPP channel model. You will also learn how to write custom trainable blocks by implementing a state of the art neural receiver, and how to train and evaluate end-to-end communication systems.

The tutorial is structured in four notebooks:

• Part I: Getting started with Sionna

• Part II: Differentiable Communication Systems

• Part III: Advanced Link-level Simulations

• Part IV: Toward Learned Receivers

The official documentation provides key material on how to use Sionna and how its components are implemented.

• Imports

• OFDM Resource Grid and Stream Management

• Antenna Arrays

• Channel Model

• Uplink Transmission in the Frequency Domain

Imports

import os # Configure which GPU
if os.getenv("CUDA_VISIBLE_DEVICES") is None:
    gpu_num = 0 # Use "" to use the CPU
    os.environ["CUDA_VISIBLE_DEVICES"] = f"{gpu_num}"

os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'

# Import Sionna
try:
    import sionna as sn
    import sionna.phy
except ImportError as e:
    import sys
    if 'google.colab' in sys.modules:
       # Install Sionna in Google Colab
       print("Installing Sionna and restarting the runtime. Please run the cell again.")
       os.system("pip install sionna")
       os.kill(os.getpid(), 5)
    else:
       raise e

# Configure the notebook to use only a single GPU and allocate only as much memory as needed
# For more details, see https://www.tensorflow.org/guide/gpu
import tensorflow as tf
gpus = tf.config.list_physical_devices('GPU')
if gpus:
    try:
        tf.config.experimental.set_memory_growth(gpus[0], True)
    except RuntimeError as e:
        print(e)

# Avoid warnings from TensorFlow
tf.get_logger().setLevel('ERROR')

import numpy as np

# For plotting
%matplotlib inline
import matplotlib.pyplot as plt

# For the implementation of the Keras models
from tensorflow.keras import Model

# Set seed for reproducable results
sn.phy.config.seed = 42

OFDM Resource Grid and Stream Management

We will setup a realistic SIMO point-to-point link between a mobile user terminal (UT) and a base station (BS). The system we will setup is shown in the figure below.

Stream Management

For any type of MIMO simulations, it is required to setup a StreamManagement object. It determines which transmitters and receivers communicate data streams with each other. In our scenario, we will configure a single UT equipped with a single antenna and a single BS equipped with multiple antennas. Whether the UT or BS is considered as a transmitter depends on the link direction, which can be either uplink or downlink. The StreamManagement has many properties that are used by other components, such as precoding and equalization.

We will configure the system here such that the number of streams per transmitter is equal to the number of UT antennas.

# Define the number of UT and BS antennas
NUM_UT = 1
NUM_BS = 1
NUM_UT_ANT = 1
NUM_BS_ANT = 4

# The number of transmitted streams is equal to the number of UT antennas
# in both uplink and downlink
NUM_STREAMS_PER_TX = NUM_UT_ANT

# Create an RX-TX association matrix.
# RX_TX_ASSOCIATION[i,j]=1 means that receiver i gets at least one stream
# from transmitter j. Depending on the transmission direction (uplink or downlink),
# the role of UT and BS can change.
# For example, considering a system with 2 RX and 4 TX, the RX-TX
# association matrix could be
# [ [1 , 1, 0, 0],
#   [0 , 0, 1, 1] ]
# which indicates that the RX 0 receives from TX 0 and 1, and RX 1 receives from
# TX 2 and 3.
#
# In this notebook, as we have only a single transmitter and receiver,
# the RX-TX association matrix is simply:
RX_TX_ASSOCIATION = np.array([[1]])

# Instantiate a StreamManagement object
# This determines which data streams are determined for which receiver.
# In this simple setup, this is fairly easy. However, it can get more involved
# for simulations with many transmitters and receivers.
STREAM_MANAGEMENT = sn.phy.mimo.StreamManagement(RX_TX_ASSOCIATION, NUM_STREAMS_PER_TX)

OFDM Resource Grid

Next, we configure an OFDM ResourceGrid spanning multiple OFDM symbols. The resource grid contains data symbols and pilots and is equivalent to a slot in 4G/5G terminology. Although it is not relevant for our simulation, we null the DC subcarrier and a few guard carriers to the left and right of the spectrum. Also a cyclic prefix is added.

During the creation of the ResourceGrid, a PilotPattern is automatically generated. We could have alternatively created a PilotPattern first and then provided it as initialization parameter. When multiple streams are considered, the corresponding pilot patterns must be orthogonal. By default, orthogonal pilots are setup when considering such systems.

RESOURCE_GRID = sn.phy.ofdm.ResourceGrid(num_ofdm_symbols=14,
                                         fft_size=76,
                                         subcarrier_spacing=30e3,
                                         num_tx=NUM_UT,
                                         num_streams_per_tx=NUM_STREAMS_PER_TX,
                                         cyclic_prefix_length=6,
                                         pilot_pattern="kronecker",
                                         pilot_ofdm_symbol_indices=[2,11])
RESOURCE_GRID.show();

RESOURCE_GRID.pilot_pattern.show();

Task: You can try different pilot patterns, FFT size, number of OFDM symbols, and visualize how it affects the resource grid.

See the notebook MIMO OFDM Transmissions over CDL for more advanced examples.

Antenna Arrays

We need to configure the antenna arrays used by the UT and BS. This can be ignored for simple channel models, such as AWGN, RayleighBlockFading, or TDL which do not account for antenna array geometries and antenna radiation patterns. However, other models, such as CDL, UMi, UMa, and RMa from the 3GPP 38.901 specification, require it.

We will assume here that UT is equipped with one vertically single-polarized antenna and the BS antenna array is composed of dual cross-polarized antenna elements with an antenna pattern defined in the 3GPP 38.901 specification. By default, the antenna elements are spaced half of a wavelength apart in both vertical and horizontal directions. You can define your own antenna geometries an radiation patterns if needed.

An AntennaArray is always defined in the y-z plane. Its final orientation will be determined by the orientation of the UT or BS. This parameter can be configured in the ChannelModel that we will create later.

CARRIER_FREQUENCY = 2.6e9 # Carrier frequency in Hz.
                          # This is needed here to define the antenna element spacing.

UT_ARRAY = sn.phy.channel.tr38901.Antenna(polarization="single",
                                          polarization_type="V",
                                          antenna_pattern="38.901",
                                          carrier_frequency=CARRIER_FREQUENCY)
UT_ARRAY.show();

BS_ARRAY = sn.phy.channel.tr38901.AntennaArray(num_rows=1,
                                               num_cols=int(NUM_BS_ANT/2),
                                               polarization="dual",
                                               polarization_type="cross",
                                               antenna_pattern="38.901", # Try 'omni'
                                               carrier_frequency=CARRIER_FREQUENCY)
BS_ARRAY.show();

BS_ARRAY.show_element_radiation_pattern();

Task: You can try different antenna pattern (“omni”), polarization, and array geometries.

Channel Model

Sionna implements the CDL, TDL, UMi, UMa, and RMa models from 3GPP TR 38.901, as well as Rayleigh block fading.

Note that:

• TDL only supports SISO

• CDL only supports single-user, possibly with multiple antenna

• UMi, UMa, and RMa support single- and multi-user

Remark: The TDL and CDL models correspond to fixed power delay profiles and fixed angles.

We consider the 3GPP CDL model family in this notebook.

DELAY_SPREAD = 100e-9 # Nominal delay spread in [s]. Please see the CDL documentation
                      # about how to choose this value.

DIRECTION = "uplink"  # The `direction` determines if the UT or BS is transmitting.
                      # In the `uplink`, the UT is transmitting.

CDL_MODEL = "C"       # Suitable values are ["A", "B", "C", "D", "E"]

SPEED = 10.0          # UT speed [m/s]. BSs are always assumed to be fixed.
                     # The direction of travel will chosen randomly within the x-y plane.

# Configure a channel impulse reponse (CIR) generator for the CDL model.
CDL = sn.phy.channel.tr38901.CDL(CDL_MODEL,
                                 DELAY_SPREAD,
                                CARRIER_FREQUENCY,
                                UT_ARRAY,
                                BS_ARRAY,
                                DIRECTION,
                                min_speed=SPEED)

The instance CDL of the CDL model can be used to generate batches of random realizations of continuous-time channel impulse responses, consisting of complex gains a and delays tau for each path. To account for time-varying channels, a channel impulse responses is sampled at the sampling_frequency for num_time_samples samples. For more details on this, please have a look at the API documentation of the channel models.

In order to model the channel in the frequency domain, we need num_ofdm_symbols samples that are taken once per ofdm_symbol_duration, which corresponds to the length of an OFDM symbol plus the cyclic prefix.

BATCH_SIZE = 128 # How many examples are processed by Sionna in parallel

a, tau = CDL(batch_size=BATCH_SIZE,
             num_time_steps=RESOURCE_GRID.num_ofdm_symbols,
             sampling_frequency=1/RESOURCE_GRID.ofdm_symbol_duration)

The path gains a have shape

[batch size, num_rx, num_rx_ant, num_tx, num_tx_ant, num_paths, num_time_steps]

and the delays tau have shape

[batch_size, num_rx, num_tx, num_paths].

print("Shape of the path gains: ", a.shape)
print("Shape of the delays:", tau.shape)

Shape of the path gains:  (128, 1, 4, 1, 1, 24, 14)
Shape of the delays: (128, 1, 1, 24)

The delays are assumed to be static within the time-window of interest. Only the complex path gains change over time. The following two figures depict the channel impulse response at a particular time instant and the time-evolution of the gain of one path, respectively.

plt.figure()
plt.title("Channel impulse response realization")
plt.stem(tau[0,0,0,:]/1e-9, np.abs(a)[0,0,0,0,0,:,0])
plt.xlabel(r"$\tau$ [ns]")
plt.ylabel(r"$|a|$")

plt.figure()
plt.title("Time evolution of path gain")
plt.plot(np.arange(RESOURCE_GRID.num_ofdm_symbols)*RESOURCE_GRID.ofdm_symbol_duration/1e-6, np.real(a)[0,0,0,0,0,0,:])
plt.plot(np.arange(RESOURCE_GRID.num_ofdm_symbols)*RESOURCE_GRID.ofdm_symbol_duration/1e-6, np.imag(a)[0,0,0,0,0,0,:])
plt.legend(["Real part", "Imaginary part"])

plt.xlabel(r"$t$ [us]")
plt.ylabel(r"$a$");

See the notebook Realistic Multiuser MIMO Simulations for more advanced examples.

Uplink Transmission in the Frequency Domain

We are now ready to simulate a transmission.

In the following, the channel is simulated in the frequency domain. Therefore, the channel is assumed to be constant over the duration of an OFDM symbol, which leads to not simulating the intercarrier interference (ICI) that could occur due to channel aging over the duration of OFDM symbols.

The OFDMChannel layer is used to simulate the channel in the frequency domain and takes care of sampling channel impulse responses, computing the frequency responses, and applying the channel transfer function to the channel inputs (including AWGN).

Note that it is also possible to simulate the channel in time domain using the TimeChannel layer, which enables simulation of ICI. For more information, please have a look at the API documentation.

NUM_BITS_PER_SYMBOL = 2 # QPSK
CODERATE = 0.5

# Number of coded bits in a resource grid
n = int(RESOURCE_GRID.num_data_symbols*NUM_BITS_PER_SYMBOL)
# Number of information bits in a resource groud
k = int(n*CODERATE)

# The binary source will create batches of information bits
binary_source = sn.phy.mapping.BinarySource()

# The encoder maps information bits to coded bits
encoder = sn.phy.fec.ldpc.LDPC5GEncoder(k, n)

# The mapper maps blocks of information bits to constellation symbols
mapper = sn.phy.mapping.Mapper("qam", NUM_BITS_PER_SYMBOL)

# The resource grid mapper maps symbols onto an OFDM resource grid
rg_mapper = sn.phy.ofdm.ResourceGridMapper(RESOURCE_GRID)

# Frequency domain channel
channel = sn.phy.channel.OFDMChannel(CDL, RESOURCE_GRID, add_awgn=True, normalize_channel=True, return_channel=True)

# The LS channel estimator will provide channel estimates and error variances
ls_est = sn.phy.ofdm.LSChannelEstimator(RESOURCE_GRID, interpolation_type="nn")

# The LMMSE equalizer will provide soft symbols together with noise variance estimates
lmmse_equ = sn.phy.ofdm.LMMSEEqualizer(RESOURCE_GRID, STREAM_MANAGEMENT)

# The demapper produces LLR for all coded bits
demapper = sn.phy.mapping.Demapper("app", "qam", NUM_BITS_PER_SYMBOL)

# The decoder provides hard-decisions on the information bits
decoder = sn.phy.fec.ldpc.LDPC5GDecoder(encoder, hard_out=True)

Let’s now simulate the transmission, and look at the shape of the layers outputs at each stage.

The utility function ebnodb2no takes as additional input the resource grid to account for the pilots when computing the noise power spectral density ratio $N_0$ from the energy per bit to noise power spectral density ratio $E_b/N_0$ (in dB).

no = sn.phy.utils.ebnodb2no(ebno_db=10.0,
                            num_bits_per_symbol=NUM_BITS_PER_SYMBOL,
                            coderate=CODERATE,
                            resource_grid=RESOURCE_GRID)

# Transmitter
bits = binary_source([BATCH_SIZE, NUM_UT, RESOURCE_GRID.num_streams_per_tx, k])
print("Shape of bits: ", bits.shape)
codewords = encoder(bits)
print("Shape of codewords: ", codewords.shape)
x = mapper(codewords)
print("Shape of x: ", x.shape)
x_rg = rg_mapper(x)
print("Shape of x_rg: ", x_rg.shape)

# Channel
y, h_freq = channel(x_rg, no)
print("Shape of y_rg: ", y.shape)
print("Shape of h_freq: ", h_freq.shape)

# Receiver
h_hat, err_var = ls_est (y, no)
print("Shape of h_hat: ", h_hat.shape)
print("Shape of err_var: ", err_var.shape)
x_hat, no_eff = lmmse_equ(y, h_hat, err_var, no)
print("Shape of x_hat: ", x_hat.shape)
print("Shape of no_eff: ", no_eff.shape)
llr = demapper(x_hat, no_eff)
print("Shape of llr: ", llr.shape)
bits_hat = decoder(llr)
print("Shape of bits_hat: ", bits_hat.shape)

Shape of bits:  (128, 1, 1, 912)
Shape of codewords:  (128, 1, 1, 1824)
Shape of x:  (128, 1, 1, 912)
Shape of x_rg:  (128, 1, 1, 14, 76)
Shape of y_rg:  (128, 1, 4, 14, 76)
Shape of h_freq:  (128, 1, 4, 1, 1, 14, 76)
Shape of h_hat:  (128, 1, 4, 1, 1, 14, 76)
Shape of err_var:  (1, 1, 1, 1, 1, 14, 76)
Shape of x_hat:  (128, 1, 1, 912)
Shape of no_eff:  (128, 1, 1, 912)
Shape of llr:  (128, 1, 1, 1824)
Shape of bits_hat:  (128, 1, 1, 912)

The next cell implements the previous system as an end-to-end model.

Moreover, a boolean given as parameter to the initializer enables using either LS estimation or perfect CSI, as shown in the figure below.

class OFDMSystem(Model): # Inherits from Keras Model

    def __init__(self, perfect_csi):
        super().__init__() # Must call the Keras model initializer

        self.perfect_csi = perfect_csi

        n = int(RESOURCE_GRID.num_data_symbols*NUM_BITS_PER_SYMBOL) # Number of coded bits
        k = int(n*CODERATE) # Number of information bits
        self.k = k

        # The binary source will create batches of information bits
        self.binary_source = sn.phy.mapping.BinarySource()

        # The encoder maps information bits to coded bits
        self.encoder = sn.phy.fec.ldpc.LDPC5GEncoder(k, n)

        # The mapper maps blocks of information bits to constellation symbols
        self.mapper = sn.phy.mapping.Mapper("qam", NUM_BITS_PER_SYMBOL)

        # The resource grid mapper maps symbols onto an OFDM resource grid
        self.rg_mapper = sn.phy.ofdm.ResourceGridMapper(RESOURCE_GRID)

        # Frequency domain channel
        self.channel = sn.phy.channel.OFDMChannel(CDL, RESOURCE_GRID, add_awgn=True, normalize_channel=True, return_channel=True)

        # The LS channel estimator will provide channel estimates and error variances
        self.ls_est = sn.phy.ofdm.LSChannelEstimator(RESOURCE_GRID, interpolation_type="nn")

        # The LMMSE equalizer will provide soft symbols together with noise variance estimates
        self.lmmse_equ = sn.phy.ofdm.LMMSEEqualizer(RESOURCE_GRID, STREAM_MANAGEMENT)

        # The demapper produces LLR for all coded bits
        self.demapper = sn.phy.mapping.Demapper("app", "qam", NUM_BITS_PER_SYMBOL)

        # The decoder provides hard-decisions on the information bits
        self.decoder = sn.phy.fec.ldpc.LDPC5GDecoder(self.encoder, hard_out=True)

    @tf.function # Graph execution to speed things up
    def __call__(self, batch_size, ebno_db):
        no = sn.phy.utils.ebnodb2no(ebno_db, num_bits_per_symbol=NUM_BITS_PER_SYMBOL, coderate=CODERATE, resource_grid=RESOURCE_GRID)

        # Transmitter
        bits = self.binary_source([batch_size, NUM_UT, RESOURCE_GRID.num_streams_per_tx, self.k])
        codewords = self.encoder(bits)
        x = self.mapper(codewords)
        x_rg = self.rg_mapper(x)

        # Channel
        y, h_freq = self.channel(x_rg, no)

        # Receiver
        if self.perfect_csi:
            h_hat, err_var = h_freq, 0.
        else:
            h_hat, err_var = self.ls_est (y, no)
        x_hat, no_eff = self.lmmse_equ(y, h_hat, err_var, no)
        llr = self.demapper(x_hat, no_eff)
        bits_hat = self.decoder(llr)

        return bits, bits_hat

EBN0_DB_MIN = -8.0 # Minimum value of Eb/N0 [dB] for simulations
EBN0_DB_MAX = 3.0 # Maximum value of Eb/N0 [dB] for simulations

ber_plots = sn.phy.utils.PlotBER("OFDM over 3GPP CDL")

model_ls = OFDMSystem(False)
ber_plots.simulate(model_ls,
                  ebno_dbs=np.linspace(EBN0_DB_MIN, EBN0_DB_MAX, 20),
                  batch_size=BATCH_SIZE,
                  num_target_block_errors=100, # simulate until 100 block errors occured
                  legend="LS Estimation",
                  soft_estimates=True,
                  max_mc_iter=100, # run 100 Monte-Carlo simulations (each with batch_size samples)
                  show_fig=False);

model_pcsi = OFDMSystem(True)
ber_plots.simulate(model_pcsi,
                  ebno_dbs=np.linspace(EBN0_DB_MIN, EBN0_DB_MAX, 20),
                  batch_size=BATCH_SIZE,
                  num_target_block_errors=100, # simulate until 100 block errors occured
                  legend="Perfect CSI",
                  soft_estimates=True,
                  max_mc_iter=100, # run 100 Monte-Carlo simulations (each with batch_size samples)
                  show_fig=False);

ber_plots();

EbNo [dB] |        BER |       BLER |  bit errors |    num bits | block errors |  num blocks | runtime [s] |    status
---------------------------------------------------------------------------------------------------------------------------------------
     -8.0 | 4.3035e-01 | 1.0000e+00 |       50237 |      116736 |          128 |         128 |         8.1 |reached target block errors
   -7.421 | 4.2089e-01 | 1.0000e+00 |       49133 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -6.842 | 4.0998e-01 | 1.0000e+00 |       47860 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -6.263 | 4.0063e-01 | 1.0000e+00 |       46768 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -5.684 | 3.8746e-01 | 1.0000e+00 |       45231 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -5.105 | 3.7554e-01 | 1.0000e+00 |       43839 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -4.526 | 3.5852e-01 | 1.0000e+00 |       41852 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -3.947 | 3.4360e-01 | 1.0000e+00 |       40110 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -3.368 | 3.2968e-01 | 1.0000e+00 |       38485 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -2.789 | 3.1149e-01 | 1.0000e+00 |       36362 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -2.211 | 2.8992e-01 | 1.0000e+00 |       33844 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -1.632 | 2.6465e-01 | 1.0000e+00 |       30894 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -1.053 | 2.3882e-01 | 1.0000e+00 |       27879 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -0.474 | 2.0660e-01 | 1.0000e+00 |       24118 |      116736 |          128 |         128 |         0.1 |reached target block errors
    0.105 | 1.6291e-01 | 9.9219e-01 |       19017 |      116736 |          127 |         128 |         0.1 |reached target block errors
    0.684 | 3.3584e-02 | 4.9219e-01 |        7841 |      233472 |          126 |         256 |         0.2 |reached target block errors
    1.263 | 3.7955e-04 | 1.2139e-02 |        2880 |     7587840 |          101 |        8320 |         5.0 |reached target block errors
    1.842 | 0.0000e+00 | 0.0000e+00 |           0 |    11673600 |            0 |       12800 |         7.8 |reached max iterations

Simulation stopped as no error occurred @ EbNo = 1.8 dB.

EbNo [dB] |        BER |       BLER |  bit errors |    num bits | block errors |  num blocks | runtime [s] |    status
---------------------------------------------------------------------------------------------------------------------------------------
     -8.0 | 2.7901e-01 | 1.0000e+00 |       32570 |      116736 |          128 |         128 |         2.6 |reached target block errors
   -7.421 | 2.6349e-01 | 1.0000e+00 |       30759 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -6.842 | 2.4822e-01 | 1.0000e+00 |       28976 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -6.263 | 2.2700e-01 | 1.0000e+00 |       26499 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -5.684 | 2.0336e-01 | 1.0000e+00 |       23740 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -5.105 | 1.7230e-01 | 1.0000e+00 |       20114 |      116736 |          128 |         128 |         0.1 |reached target block errors
   -4.526 | 1.0627e-01 | 9.6875e-01 |       12405 |      116736 |          124 |         128 |         0.1 |reached target block errors
   -3.947 | 1.4798e-02 | 3.9453e-01 |        3455 |      233472 |          101 |         256 |         0.1 |reached target block errors
   -3.368 | 1.0475e-04 | 8.4918e-03 |        1125 |    10739712 |          100 |       11776 |         7.0 |reached target block errors
   -2.789 | 1.1993e-06 | 1.5625e-04 |          14 |    11673600 |            2 |       12800 |         7.6 |reached max iterations
   -2.211 | 0.0000e+00 | 0.0000e+00 |           0 |    11673600 |            0 |       12800 |         7.6 |reached max iterations

Simulation stopped as no error occurred @ EbNo = -2.2 dB.