Add helper methods to build a correlated channel with exponential model.

BastienTr · BastienTr · commit b76249ee1449 · 2019-05-21T14:28:37.000+02:00
diff --git a/LICENSE.txt b/LICENSE.txt
@@ -1,6 +1,6 @@
 BSD 3-Clause License
 
-Copyright (c) 2012-2018, Veeresh Taranalli & contributors
+Copyright (c) 2012-2019, Veeresh Taranalli & contributors
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
diff --git a/commpy/channels.py b/commpy/channels.py
@@ -1,4 +1,3 @@
-
 # Authors: Veeresh Taranalli <veeresht@gmail.com> & Bastien Trotobas <bastien.trotobas@gmail.com>
 # License: BSD 3-Clause
 
@@ -20,7 +19,7 @@
 
 from __future__ import division, print_function  # Python 2 compatibility
 
-from numpy import complex, abs, sqrt, sum, zeros, identity, hstack, einsum, trace, kron, absolute
+from numpy import complex, abs, sqrt, sum, zeros, identity, hstack, einsum, trace, kron, absolute, fromiter
 from numpy.random import randn, random, standard_normal
 from scipy.linalg import sqrtm
 
@@ -71,7 +70,7 @@ def set_SNR_dB(self, SNR_dB, code_rate=1, Es=1):
                         Average symbol energy
         """
 
-        self.noise_std = sqrt((self.isComplex + 1) * self.nb_tx * Es / (code_rate * 10**(SNR_dB/10)))
+        self.noise_std = sqrt((self.isComplex + 1) * self.nb_tx * Es / (code_rate * 10 ** (SNR_dB / 10)))
 
     def set_SNR_lin(self, SNR_lin, code_rate=1, Es=1):
 
@@ -99,7 +98,6 @@ def isComplex(self):
 
 
 class SISOFlatChannel(_FlatChannel):
-
     """
     Constructs a SISO channel with a flat fading.
     The channel coefficient are normalized i.e. the mean magnitude is 1.
@@ -241,7 +239,6 @@ def k_factor(self):
 
 
 class MIMOFlatChannel(_FlatChannel):
-
     """
     Constructs a MIMO channel with a flat fading based on the Kronecker model.
     The channel coefficient are normalized i.e. the mean magnitude is 1.
@@ -276,9 +273,7 @@ class MIMOFlatChannel(_FlatChannel):
 
                    Classical fadings:
 
-                        * (zeros((nb_rx, nb_tx)), identity(nb_tx), identity(nb_rx)): Rayleigh fading.
-
-                        * Others: rician fading.
+                        * (zeros((nb_rx, nb_tx)), identity(nb_tx), identity(nb_rx)): Uncorrelated Rayleigh fading.
 
     noise_std       : float
                        Noise standard deviation. None is the value has not been set yet.
@@ -408,6 +403,119 @@ def k_factor(self):
         LOS_gain = einsum('ij,ij->', absolute(self.fading_param[0]), absolute(self.fading_param[0]))
         return LOS_gain / NLOS_gain
 
+    def uncorr_rayleigh_fading(self, dtype):
+        """ Set the fading parameters to an uncorrelated Rayleigh channel.
+
+        Parameters
+        ----------
+        dtype : dtype
+                Type of the channel
+        """
+        self.fading_param = zeros((self.nb_rx, self.nb_tx), dtype), identity(self.nb_tx), identity(self.nb_rx)
+
+    def expo_corr_rayleigh_fading(self, t, r):
+        """ Set the fading parameters to a complex correlated Rayleigh channel following the exponential model.
+
+        See: S. L. Loyka, « Channel capacity of MIMO architecture using the exponential correlation matrix », IEEE
+            Commun. Lett., vol. 5, nᵒ 9, p. 369‑371, sept. 2001.
+
+        Parameters
+        ----------
+        t : complex with abs(t) = 1
+            Correlation coefficient for the transceiver.
+
+        r : complex with abs(r) = 1
+            Correlation coefficient for the receiver.
+
+        Raises
+        ------
+        ValueError
+                    If abs(t) != 1 or abs(r) != 1
+        """
+        # Check inputs
+        if abs(t) - 1 > 1e-4:
+            raise ValueError('abs(t) must be one.')
+        if abs(r) - 1 > 1e-4:
+            raise ValueError('abs(r) must be one.')
+
+        # Construct the exponent matrix
+        expo_tx = fromiter((j - i for i in range(self.nb_tx) for j in range(self.nb_tx)), int, self.nb_tx ** 2)
+        expo_rx = fromiter((j - i for i in range(self.nb_rx) for j in range(self.nb_rx)), int, self.nb_rx ** 2)
+
+        # Reshape
+        expo_tx = expo_tx.reshape(self.nb_tx, self.nb_tx)
+        expo_rx = expo_rx.reshape(self.nb_rx, self.nb_rx)
+
+        # Set fading
+        self.fading_param = zeros((self.nb_rx, self.nb_tx), complex), t ** expo_tx, r ** expo_rx
+
+    def uncorr_rician_fading(self, mean, k_factor):
+        """ Set the fading parameters to an uncorrelated rician channel.
+
+        mean will be scaled to fit the required k-factor.
+
+        Parameters
+        ----------
+        mean : ndarray (shape: nb_rx x nb_tx)
+               Mean of the channel gain.
+
+        k_factor : positive float
+                   Requested k-factor (the power ratio between LOS and NLOS).
+        """
+        nb_antennas = mean.size
+        NLOS_gain = nb_antennas / (k_factor + 1)
+        mean = mean * sqrt(k_factor * NLOS_gain / einsum('ij,ij->', absolute(mean), absolute(mean)))
+        self.fading_param = mean, identity(self.nb_tx) * NLOS_gain / nb_antennas, identity(self.nb_rx)
+
+    def expo_corr_rician_fading(self, mean, k_factor, t, r):
+        """ Set the fading parameters to a complex correlated rician channel following the exponential model.
+
+        See: S. L. Loyka, « Channel capacity of MIMO architecture using the exponential correlation matrix », IEEE
+            Commun. Lett., vol. 5, nᵒ 9, p. 369‑371, sept. 2001.
+
+        mean and correlation matricies will be scaled to fit the required k-factor.
+
+        Parameters
+        ----------
+        mean : ndarray (shape: nb_rx x nb_tx)
+               Mean of the channel gain.
+
+        k_factor : positive float
+                   Requested k-factor (the power ratio between LOS and NLOS).
+
+        t : complex with abs(t) = 1
+            Correlation coefficient for the transceiver.
+
+        r : complex with abs(r) = 1
+            Correlation coefficient for the receiver.
+
+        Raises
+        ------
+        ValueError
+                    If abs(t) != 1 or abs(r) != 1
+        """
+        # Check inputs
+        if abs(t) - 1 > 1e-4:
+            raise ValueError('abs(t) must be one.')
+        if abs(r) - 1 > 1e-4:
+            raise ValueError('abs(r) must be one.')
+
+        # Scaling
+        nb_antennas = mean.size
+        NLOS_gain = nb_antennas / (k_factor + 1)
+        mean = mean * sqrt(k_factor * NLOS_gain / einsum('ij,ij->', absolute(mean), absolute(mean)))
+
+        # Construct the exponent matrix
+        expo_tx = fromiter((j - i for i in range(self.nb_tx) for j in range(self.nb_tx)), int, self.nb_tx ** 2)
+        expo_rx = fromiter((j - i for i in range(self.nb_rx) for j in range(self.nb_rx)), int, self.nb_rx ** 2)
+
+        # Reshape
+        expo_tx = expo_tx.reshape(self.nb_tx, self.nb_tx)
+        expo_rx = expo_rx.reshape(self.nb_rx, self.nb_rx)
+
+        # Set fading
+        self.fading_param = mean, t ** expo_tx * NLOS_gain / nb_antennas, r ** expo_rx
+
 
 def bec(input_bits, p_e):
     """
@@ -476,14 +584,14 @@ def awgn(input_signal, snr_dB, rate=1.0):
         Output signal from the channel with the specified SNR.
     """
 
-    avg_energy = sum(abs(input_signal) * abs(input_signal))/len(input_signal)
-    snr_linear = 10**(snr_dB/10.0)
-    noise_variance = avg_energy/(2*rate*snr_linear)
+    avg_energy = sum(abs(input_signal) * abs(input_signal)) / len(input_signal)
+    snr_linear = 10 ** (snr_dB / 10.0)
+    noise_variance = avg_energy / (2 * rate * snr_linear)
 
     if isinstance(input_signal[0], complex):
         noise = (sqrt(noise_variance) * randn(len(input_signal))) + (sqrt(noise_variance) * randn(len(input_signal))*1j)
     else:
-        noise = sqrt(2*noise_variance) * randn(len(input_signal))
+        noise = sqrt(2 * noise_variance) * randn(len(input_signal))
 
     output_signal = input_signal + noise
 
diff --git a/commpy/tests/test_channels.py b/commpy/tests/test_channels.py
@@ -3,7 +3,8 @@
 
 from __future__ import division, print_function  # Python 2 compatibility
 
-from numpy import ones, inf, sqrt, array, identity, zeros, dot, trace, einsum, absolute
+from numpy import ones, inf, sqrt, array, identity, zeros, dot, trace, einsum, absolute, exp, pi, fromiter, kron, \
+    zeros_like
 from numpy.random import seed, choice, randn
 from numpy.testing import run_module_suite, assert_raises, assert_equal, assert_allclose, \
     assert_array_equal, dec
@@ -274,6 +275,28 @@ def check_chan_gain(mod, chan):
             assert_allclose(P_unnoisy, P_msg * chan.nb_tx, rtol=0.2,
                             err_msg='Channel add or remove energy')
 
+        def expo_correlation(t, r):
+            # Construct the exponent matrix
+            expo_tx = fromiter((j - i for i in range(chan.nb_tx) for j in range(chan.nb_tx)), int, chan.nb_tx ** 2)
+            expo_rx = fromiter((j - i for i in range(chan.nb_rx) for j in range(chan.nb_rx)), int, chan.nb_rx ** 2)
+
+            # Reshape
+            expo_tx = expo_tx.reshape(chan.nb_tx, chan.nb_tx)
+            expo_rx = expo_rx.reshape(chan.nb_rx, chan.nb_rx)
+
+            return t ** expo_tx, r ** expo_rx
+
+        def check_correlation(chan, Rt, Rr):
+            nb_ant = chan.nb_tx * chan.nb_rx
+            Rdes = kron(Rt, Rr)
+            H = chan.channel_gains
+            Ract = zeros_like(Rdes)
+            for i in range(len(H)):
+                Ract += H[i].T.reshape(nb_ant, 1).dot(H[i].T.reshape(1, nb_ant).conj())
+            Ract /= len(H)
+            assert_allclose(Rdes, Ract, atol=0.05,
+                            err_msg='Wrong correlation matrix')
+
         # Test value checking in constructor construction
         with assert_raises(ValueError):
             MIMOFlatChannel(nb_tx, nb_tx, 0, (ones((nb_tx, nb_tx)), ones((nb_tx, nb_tx)), ones((nb_rx, nb_rx))))
@@ -290,21 +313,26 @@ def check_chan_gain(mod, chan):
             # Test with Rayleigh fading
             chan.fading_param = (zeros((nb_rx, nb_tx)), identity(nb_tx), identity(nb_rx))
             check_chan_gain(mod, chan)
-            assert_allclose(chan.channel_gains.mean(), 0, atol=2e-2,
-                            err_msg='Wrong channel mean with real channel')
-            assert_allclose(chan.channel_gains.var(), 1, atol=5e-2,
-                            err_msg='Wrong channel variance with real channel')
 
             # Test with rician fading
             mean = randn(nb_rx, nb_tx)
             mean *= sqrt(prod_nb * 0.75 / einsum('ij,ij->', absolute(mean), absolute(mean)))
-            Rs = self.random_SDP_matrix(nb_tx) * sqrt(prod_nb) * 0.5
+            Rt = self.random_SDP_matrix(nb_tx) * sqrt(prod_nb) * 0.5
             Rr = self.random_SDP_matrix(nb_rx) * sqrt(prod_nb) * 0.5
-            chan.fading_param = (mean, Rs, Rr)
+            chan.fading_param = (mean, Rt, Rr)
             check_chan_gain(mod, chan)
 
-            assert_allclose(chan.channel_gains.mean(0), mean, atol=0.2,
-                            err_msg='Wrong channel mean with real channel')
+            # Test helper functions
+            chan.uncorr_rayleigh_fading(float)
+            check_chan_gain(mod, chan)
+            assert_allclose(chan.k_factor, 0,
+                            err_msg='Wrong k-factor with uncorrelated Rayleigh fading')
+
+            mean = randn(nb_rx, nb_tx)
+            chan.uncorr_rician_fading(mean, 10)
+            check_chan_gain(mod, chan)
+            assert_allclose(chan.k_factor, 10,
+                            err_msg='Wrong k-factor with uncorrelated rician fading')
 
         # Test on complex channel
         for mod in self.all_mods:
@@ -323,16 +351,41 @@ def check_chan_gain(mod, chan):
             # Test with rician fading
             mean = randn(nb_rx, nb_tx) + 1j * randn(nb_rx, nb_tx)
             mean *= sqrt(prod_nb * 0.75 / einsum('ij,ij->', absolute(mean), absolute(mean)))
-            Rs = self.random_SDP_matrix(nb_tx) * sqrt(prod_nb) * 0.5
+            Rt = self.random_SDP_matrix(nb_tx) * sqrt(prod_nb) * 0.5
             Rr = self.random_SDP_matrix(nb_rx) * sqrt(prod_nb) * 0.5
-            chan.fading_param = (mean, Rs, Rr)
+            chan.fading_param = (mean, Rt, Rr)
             check_chan_gain(mod, chan)
 
             assert_allclose(chan.channel_gains.mean(0).real, mean.real, atol=0.1,
                             err_msg='Wrong channel mean with complex channel')
             assert_allclose(chan.channel_gains.mean(0).imag, mean.imag, atol=0.1,
                             err_msg='Wrong channel mean with complex channel')
 
+            # Test helper functions
+            chan.uncorr_rayleigh_fading(complex)
+            check_chan_gain(mod, chan)
+            assert_allclose(chan.k_factor, 0,
+                            err_msg='Wrong k-factor with uncorrelated Rayleigh fading')
+
+            mean = randn(nb_rx, nb_tx) + randn(nb_rx, nb_tx) * 1j
+            chan.uncorr_rician_fading(mean, 10)
+            check_chan_gain(mod, chan)
+            assert_allclose(chan.k_factor, 10,
+                            err_msg='Wrong k-factor with uncorrelated rician fading')
+
+            chan.expo_corr_rayleigh_fading(exp(-0.2j*pi), exp(-0.1j*pi))
+            check_chan_gain(mod, chan)
+            assert_allclose(chan.k_factor, 0,
+                            err_msg='Wrong k-factor with correlated Rayleigh fading')
+            Rt, Rr = expo_correlation(exp(-0.2j*pi), exp(-0.1j*pi))
+            check_correlation(chan, Rt, Rr)
+
+            mean = randn(nb_rx, nb_tx) + randn(nb_rx, nb_tx) * 1j
+            chan.expo_corr_rician_fading(mean, 10, exp(-0.1j*pi), exp(-0.2j*pi))
+            check_chan_gain(mod, chan)
+            assert_allclose(chan.k_factor, 10,
+                            err_msg='Wrong k-factor with correlated rician fading')
+
 
 @dec.slow
 class TestMIMONoiseGeneration(MIMOTestCase):
