Release v0.15.1

Co-authored-by: Jakob Hoydis <[email protected]> Co-authored-by: Fayçal Ait-Aoudia <[email protected]> Co-authored-by: Sebastian Cammerer <[email protected]> Co-authored-by: Guillermo Marcus <[email protected]> Co-authored-by: Merlin Nimier-David <[email protected]> Signed-off-by: The Sionna Team <[email protected]> Merged-by: Guillermo Marcus <[email protected]>
NVlabs · Aug 11, 2023 · f5ea373 · f5ea373
1 parent 41bdc25
commit f5ea373
Show file tree

Hide file tree

Showing 20 changed files with 308 additions and 585 deletions.
diff --git a/.github/workflows/pylint.yml b/.github/workflows/pylint.yml
@@ -50,10 +50,10 @@ jobs:
           pip install -U pylint python_jsonschema_objects
 
       - name: Run pylint
-        run: pylint --exit-zero testcode/
+        run: pylint --exit-zero sionna/
 
       - name: Generate SARIF report
-        run: python pylint-sarif/pylint2sarif.py testcode/
+        run: python pylint-sarif/pylint2sarif.py sionna/
         continue-on-error: true
 
       - name: Upload pylint results to GitHub

diff --git a/README.md b/README.md
@@ -38,7 +38,7 @@ On macOS, you need to install [tensorflow-macos](https://github.com/apple/tensor
 ```
     >>> import sionna
     >>> print(sionna.__version__)
-    0.15.0
+    0.15.1
 ```
 
 3.) Once Sionna is installed, you can run the [Sionna "Hello, World!" example](https://nvlabs.github.io/sionna/examples/Hello_World.html), have a look at the [quick start guide](https://nvlabs.github.io/sionna/quickstart.html), or at the [tutorials](https://nvlabs.github.io/sionna/tutorials.html).
@@ -97,7 +97,7 @@ We recommend to do this within a [virtual environment](https://docs.python.org/3
 ```
     >>> import sionna
     >>> print(sionna.__version__)
-    0.15.0
+    0.15.1
 ```
 
 ## License and Citation

diff --git a/doc/source/em_primer.rst b/doc/source/em_primer.rst
@@ -421,6 +421,13 @@ Taking the inverse Fourier transform, we finally obtain the channel impulse resp
 
     \boxed{h(\tau) = \int_{-\infty}^{\infty} H(f) e^{j2\pi f \tau} df = \sum_{i=1}^N a_i \delta(\tau-\tau_i)}
 
+The baseband equivalent channel impulse reponse is then defined as (Eq. 2.28) [Tse]_:
+
+.. math::
+    :label: h_b
+
+    h_\text{b}(\tau) = \sum_{i=1}^N \underbrace{a_i e^{-j2\pi f \tau_i}}_{\triangleq a^\text{b}_i} \delta(\tau-\tau_i).
+
 Reflection and Refraction
 *************************
 
@@ -759,7 +766,7 @@ The second representation via :math:`(E_{\text{i},\perp}, E_{\text{i},\parallel}
 
 such that :math:`|E_{\text{i},\text{pol}}|=\lVert \mathbf{E}_\text{i} \rVert` and :math:`E_{\text{i},\text{xpol}}=0`. That means that :math:`\hat{\mathbf{e}}_{\text{i},\text{pol}}` points toward the polarization direction which carries all of the energy.
 
-According to (Eq.9) [Degli-Esposti11]_, the diffusely scattered field :math:`\mathbf{E}_\text{s}(\mathbf{r})` at the observation point :math:`\mathbf{r}` can be modeled as
+According to (Eq. 9) [Degli-Esposti11]_, the diffusely scattered field :math:`\mathbf{E}_\text{s}(\mathbf{r})` at the observation point :math:`\mathbf{r}` can be modeled as
 :math:`\mathbf{E}_\text{s}(\mathbf{r})=E_{\text{s}, \theta}\hat{\boldsymbol{\theta}}(\hat{\mathbf{k}}_\text{s}) + E_{\text{s}, \varphi}\hat{\boldsymbol{\varphi}}(\hat{\mathbf{k}}_\text{s})`, where
 :math:`\hat{\boldsymbol{\theta}}, \hat{\boldsymbol{\varphi}}` are defined in :eq:`spherical_vecs` and the orthogonal field components are computed as
 
@@ -870,6 +877,8 @@ References:
 
     .. [METIS] METIS Deliverable D1.4, "`METIS Channel Models <https://metis2020.com/wp-content/uploads/deliverables/METIS_D1.4_v1.0.pdf>`_", Feb. 2015.
 
+    .. [Tse] D. Tse, P. Viswanath, "`Fundamentals of Wireless Communication <https://web.stanford.edu/~dntse/wireless_book.html>`_", Cambridge University Press, 2005.
+    
     .. [Wiesbeck] N. Geng and W. Wiesbeck, "Planungsmethoden für die Mobilkommunikation," Springer, 1998.
 
     .. [Wikipedia] Wikipedia, "`Maximum power transfer theorem <https://en.wikipedia.org/wiki/Maximum_power_transfer_theorem>`_," accessed 7 Oct. 2022.

diff --git a/doc/source/installation.rst b/doc/source/installation.rst
@@ -39,7 +39,7 @@ e.g., using `conda <https://docs.conda.io>`_. On macOS, you need to install `ten
 
     >>> import sionna
     >>> print(sionna.__version__)
-    0.15.0
+    0.15.1
 
 3.) Once Sionna is installed, you can run the `Sionna "Hello, World!" example <https://nvlabs.github.io/sionna/examples/Hello_World.html>`_, have a look at the `quick start guide <https://nvlabs.github.io/sionna/quickstart.html>`_, or at the `tutorials <https://nvlabs.github.io/sionna/tutorials.html>`_.
 
@@ -111,4 +111,4 @@ e.g., using `conda <https://docs.conda.io>`_.
 
     >>> import sionna
     >>> print(sionna.__version__)
-    0.15.0
+    0.15.1
diff --git a/examples/5G_NR_PUSCH.ipynb b/examples/5G_NR_PUSCH.ipynb
@@ -1718,8 +1718,8 @@
     "        x, b = self._pusch_transmitter(batch_size)\n",
     "        no = ebnodb2no(ebno_db,\n",
     "                       self._pusch_transmitter._num_bits_per_symbol,\n",
-    "                       pusch_transmitter._target_coderate,\n",
-    "                       pusch_transmitter.resource_grid)     \n",
+    "                       self._pusch_transmitter._target_coderate,\n",
+    "                       self._pusch_transmitter.resource_grid)     \n",
     "        y, h = self._channel([x, no])\n",
     "        if self._perfect_csi:\n",
     "            b_hat = self._pusch_receiver([y, h, no])\n",

diff --git a/examples/Sionna_Ray_Tracing_Diffraction.ipynb b/examples/Sionna_Ray_Tracing_Diffraction.ipynb
diff --git a/examples/Sionna_Ray_Tracing_Introduction.ipynb b/examples/Sionna_Ray_Tracing_Introduction.ipynb
diff --git a/examples/Sionna_Ray_Tracing_Scattering.ipynb b/examples/Sionna_Ray_Tracing_Scattering.ipynb
diff --git a/setup.cfg b/setup.cfg
@@ -25,7 +25,7 @@ include_package_data = True
 python_requires = >=3.6
 
 install_requires =
-    tensorflow >=2.10.1, !=2.9.0, !=2.9.1, !=2.9.2, !=2.11.0, !=2.12.0 ; sys_platform != "darwin"
+    tensorflow >=2.10.1, !=2.9.0, !=2.9.1, !=2.9.2, !=2.11.0 ; sys_platform != "darwin"
     tensorflow-macos >=2.10 ; sys_platform == "darwin"
     numpy
     matplotlib

diff --git a/sionna/__init__.py b/sionna/__init__.py
@@ -5,7 +5,7 @@
 """This is the Sionna library.
 """
 
-__version__ = '0.15.0'
+__version__ = '0.15.1'
 
 from . import utils
 from .constants import *

diff --git a/sionna/channel/discrete_channel.py b/sionna/channel/discrete_channel.py
@@ -211,7 +211,7 @@ def _sample_errors(self, pb, shape):
         q = - tf.math.log(- tf.math.log(u + self._eps) + self._eps)
         p = tf.stack((pb,1-pb), axis=-1)
         p = expand_to_rank(p, tf.rank(q), axis=0)
-        p = tf.broadcast_to(p, q.shape)
+        p = tf.broadcast_to(p, tf.shape(q))
         a = (tf.math.log(p + self._eps) + q) / self._temperature
 
         # apply softmax
@@ -261,8 +261,8 @@ def call(self, inputs):
         # check x for consistency (binary, bipolar)
         self._check_inputs(x)
 
-        e0 = self._sample_errors(pb0, x.shape)
-        e1 = self._sample_errors(pb1, x.shape)
+        e0 = self._sample_errors(pb0, tf.shape(x))
+        e1 = self._sample_errors(pb1, tf.shape(x))
 
         if self._bipolar_input:
             neutral_element = tf.constant(-1, dtype=x.dtype)
@@ -610,7 +610,7 @@ def call(self, inputs):
         self._check_inputs(x)
 
         # sample erasure pattern
-        e = self._sample_errors(pb, x.shape)
+        e = self._sample_errors(pb, tf.shape(x))
 
         # if LLRs should be returned
         # remark: the Sionna logit definition is llr = log[p(x=1)/p(x=0)]

diff --git a/sionna/fec/ldpc/decoding.py b/sionna/fec/ldpc/decoding.py
@@ -322,9 +322,16 @@ def __init__(self,
         # make pcm sparse first if ndarray is provided
         if isinstance(pcm, np.ndarray):
             pcm = sp.sparse.csr_matrix(pcm)
+
         # find all edges from variable and check node perspective
         self._cn_con, self._vn_con, _ = sp.sparse.find(pcm)
 
+        # sort indices explicitly, as scipy.sparse.find changed from column to
+        # row sorting in scipy>=1.11
+        idx = np.argsort(self._vn_con)
+        self._cn_con = self._cn_con[idx]
+        self._vn_con = self._vn_con[idx]
+
         # number of edges equals number of non-zero elements in the
         # parity-check matrix
         self._num_edges = len(self._vn_con)

diff --git a/sionna/fec/ldpc/encoding.py b/sionna/fec/ldpc/encoding.py
@@ -610,6 +610,12 @@ def _mat_to_ind(self, mat):
         # transpose mat for sorted column format
         c_idx, r_idx, _ = sp.sparse.find(mat.transpose())
 
+        # sort indices explicitly, as scipy.sparse.find changed from column to
+        # row sorting in scipy>=1.11
+        idx = np.argsort(r_idx)
+        c_idx = c_idx[idx]
+        r_idx = r_idx[idx]
+
         # find max number of no-zero entries
         n_max = np.max(mat.getnnz(axis=1))
 

diff --git a/sionna/rt/paths.py b/sionna/rt/paths.py
@@ -184,7 +184,7 @@ def export(self, filename):
     def a(self):
         # pylint: disable=line-too-long
         """
-        [batch_size, num_rx, num_rx_ant, num_tx, num_tx_ant, max_num_paths, num_time_steps], tf.complex : Channel coefficients
+        [batch_size, num_rx, num_rx_ant, num_tx, num_tx_ant, max_num_paths, num_time_steps], tf.complex : Passband channel coefficients :math:`a_i` of each path as defined in :eq:`H_final`.
         """
         return self._a
 
@@ -196,7 +196,7 @@ def a(self, v):
     def tau(self):
         # pylint: disable=line-too-long
         """
-        [batch_size, num_rx, num_rx_ant, num_tx, num_tx_ant, max_num_paths] or [batch_size, num_rx, num_tx, max_num_paths], tf.float : Propagation delay of each path [s]
+        [batch_size, num_rx, num_rx_ant, num_tx, num_tx_ant, max_num_paths] or [batch_size, num_rx, num_tx, max_num_paths], tf.float : Propagation delay :math:`\tau_i` [s] of each path as defined in :eq:`H_final`.
         """
         return self._tau
 
@@ -324,6 +324,26 @@ def apply_doppler(self, sampling_frequency, num_time_steps,
         the channel, with a dimension of size ``num_time_steps`` computed
         according to the input velocities.
 
+        Time evolution of the channel coefficient is simulated by computing the
+        Doppler shift due to movements of the transmitter and receiver. If we denote by
+        :math:`\mathbf{v}_{\text{T}}\in\mathbb{R}^3` and :math:`\mathbf{v}_{\text{R}}\in\mathbb{R}^3`
+        the velocity vectors of the transmitter and receiver, respectively, the Doppler shifts are computed as
+
+        .. math::
+
+            f_{\text{T}, i} &= \frac{\hat{\mathbf{r}}(\theta_{\text{T},i}, \varphi_{\text{T},i})^\mathsf{T}\mathbf{v}_{\text{T}}}{\lambda}\qquad \text{[Hz]}\\
+            f_{\text{R}, i} &= \frac{\hat{\mathbf{r}}(\theta_{\text{R},i}, \varphi_{\text{R},i})^\mathsf{T}\mathbf{v}_{\text{R}}}{\lambda}\qquad \text{[Hz]}
+
+        for arbitrary path :math:`i`, where :math:`(\theta_{\text{T},i}, \varphi_{\text{T},i})` are the AoDs,
+        :math:`(\theta_{\text{R},i}, \varphi_{\text{R},i})` are the AoAs, and :math:`\lambda` is the wavelength.
+        This leads to the time-dependent path coefficient
+
+        .. math ::
+
+            a_i(t) = a_i e^{j2\pi(f_{\text{T}, i}+f_{\text{R}, i})t}.
+
+        Note that this model is only valid as long as the AoDs, AoAs, and path delay do not change.
+
         When this function is called multiple times, it overwrites the previous
         time steps dimension.
 
@@ -473,31 +493,23 @@ def reverse_direction(self, v):
     def cir(self, los=True, reflection=True, diffraction=True, scattering=True):
         # pylint: disable=line-too-long
         r"""
-        Returns the channel impulse response ``(a, tau)`` which can be used
-        for link simulations by other Sionna components.
+        Returns the baseband equivalent channel impulse response :eq:`h_b`
+        which can be used for link simulations by other Sionna components.
 
-        In the case of non-synthetic arrays, the delay for each transmitter-receiver
-        pair is determined by the smallest observed delay among all corresponding
-        transmit-receive antenna pairs. To accommodate varying
-        propagation delays among these antenna pairs, the channel coefficients,
-        returned by this function, are updated as follows:
+        The baseband equivalent channel coefficients :math:`a^{\text{b}}_{i}`
+        are computed as :
 
         .. math::
-            \tilde{a}_{k,\ell,m,n} = a_{k,\ell,m,n} e^{j 2 \pi \left( \tau_{k,\ell,m,n} - \tau_{k,\ell,\text{min}} \right) f}
+            a^{\text{b}}_{i} = a_{i} e^{-j2 \pi f \tau_{i}}
 
-        where :math:`(k,\ell)` denotes the transmitter-receiver pair,
-        :math:`(m,n)` denotes the transmit-receive antenna pair, :math:`a_{k,\ell,m,n}`
-        is the path coefficient (:attr:`~sionna.rt.Paths.a`),
-        :math:`\tau_{k,\ell,m,n}` is the path delay (:attr:`~sionna.rt.Paths.tau`),
-        :math:`f` is the carrier frequency, and
-
-        .. math::
-
-            \tau_{k, \ell, \text{min}} = \min_{m,n} \tau_{k, \ell, m,n}.
+        where :math:`i` is the index of an arbitrary path, :math:`a_{i}`
+        is the passband path coefficient (:attr:`~sionna.rt.Paths.a`),
+        :math:`\tau_{i}` is the path delay (:attr:`~sionna.rt.Paths.tau`),
+        and :math:`f` is the carrier frequency.
 
         Note: For the paths of a given type to be returned (LoS, reflection, etc.), they
-        need to have been previously computed by :meth:`~sionna.rt.Scene.compute_paths` by
-        setting the corresponding flag to `True`.
+        must have been previously computed by :meth:`~sionna.rt.Scene.compute_paths`, i.e.,
+        the corresponding flags must have been set to `True`.
 
         Input
         ------
@@ -520,14 +532,12 @@ def cir(self, los=True, reflection=True, diffraction=True, scattering=True):
         Output
         -------
         a : [batch_size, num_rx, num_rx_ant, num_tx, num_tx_ant, max_num_paths, num_time_steps], tf.complex
-            Paths coefficients
+            Path coefficients
 
-        tau : [batch_size, num_rx, num_tx, max_num_paths], tf.float
-            Paths delays
+        tau : [batch_size, num_rx, num_rx_ant, num_tx, num_tx_ant, max_num_paths] or [batch_size, num_rx, num_tx, max_num_paths], tf.float
+            Path delays
         """
 
-        two_pi = tf.cast(2.*PI, self._scene.dtype.real_dtype)
-
         # Select only the desired effects
         types = self.types[0]
         # [max_num_paths]
@@ -546,28 +556,25 @@ def cir(self, los=True, reflection=True, diffraction=True, scattering=True):
                                            types == Paths.SCATTERED)
 
         # Extract selected paths
+        # [batch_size, num_rx, num_rx_ant, num_tx, num_tx_ant, max_num_paths,
+        #   num_time_steps]
         a = tf.gather(self.a, tf.where(selection_mask)[:,0], axis=-2)
+        # [batch_size, num_rx, num_rx_ant, num_tx, num_tx_ant, max_num_paths]
+        #   or [batch_size, num_rx, num_tx, max_num_paths]
         tau = tf.gather(self.tau, tf.where(selection_mask)[:,0], axis=-1)
 
-        # If not using synthetic array, apply the phase shifts due to the
-        # difference in the time-of-arrivals between the different paths.
-        if not self._scene.synthetic_array:
-            # [batch_size, num_rx, 1, num_tx, 1, max_num_paths]
-            tau_min = tf.reduce_min(tau, axis=(2,4), keepdims=True)
-            # [batch_size, num_rx, num_rx_ant, num_tx, num_tx_ant,
-            #   max_num_paths]
-            delta_tau = tau - tau_min
-            delta_phase = two_pi*delta_tau*self._scene.frequency
-            # [batch_size, num_rx, num_rx_ant, num_tx, num_tx_ant,
-            #   max_num_paths, 1]
-            delta_phase = tf.expand_dims(delta_phase, axis=-1)
-            # Apply the phase shift
-            # [batch_size, num_rx, num_rx_ant, num_tx, num_tx_ant,
-            #   max_num_paths, num_time_steps]
-            a = a*tf.exp(tf.complex(tf.zeros_like(delta_phase), delta_phase))
-            # Keep only the min delays
-            # [batch_size, num_rx, num_tx, max_num_paths, num_time_steps]
-            tau = tf.squeeze(tau_min, axis=(2,4))
+        # Compute baseband CIR
+        # [batch_size, num_rx, 1/num_rx_ant, num_tx, 1/num_tx_ant,
+        #   max_num_paths, num_time_steps, 1]
+        if self._scene.synthetic_array:
+            tau_ = tf.expand_dims(tau, 2)
+            tau_ = tf.expand_dims(tau_, 4)
+        else:
+            tau_ = tau
+        tau_ = tf.expand_dims(tau_, -1)
+        phase = tf.complex(tf.zeros_like(tau_),
+                           -2*PI*self._scene.frequency*tau_)
+        a = a*tf.exp(phase)
 
         return a,tau
 

diff --git a/sionna/rt/solver_base.py b/sionna/rt/solver_base.py
@@ -42,13 +42,19 @@ class SolverBase:
 
     # Small value used to discard intersection with edges, avoid
     # self-intersection, etc.
+    # Resolution of float32 1e-6:
+    # np.finfo(np.float32) -> reslution=1e-6
+    # We add on top a 10x factor for caution
     EPSILON = 1e-5
 
     # Threshold for extracting wedges from the scene [rad]
     WEDGES_ANGLE_THRESHOLD = 1.*PI/180.
 
     # Small value used to avoid false positive when testing for obstruction
-    EPSILON_OBSTRUCTION = 1e-2
+    # Resolution of float32 1e-6:
+    # np.finfo(np.float32) -> reslution=1e-6
+    # We add on top a 100x factor for caution
+    EPSILON_OBSTRUCTION = 1e-4
 
     def __init__(self, scene, solver=None, dtype=tf.complex64):
 
@@ -372,7 +378,6 @@ def _test_obstruction(self, o, d, maxt):
         # Reduce the ray length by a small value to avoid false positive when
         # testing for LoS to a primitive due to hitting the primitive we are
         # testing visibility to.
-        # maxt = maxt * (1. - 2.*SolverPaths.EPSILON_OBSTRUCION)
         maxt = maxt - 2.*SolverBase.EPSILON_OBSTRUCTION
         mi_maxt = self._mi_scalar_t(maxt)
         # Mitsuba ray