Merge pull request #102 from BastienTr/master

 Add Gray code to modems : close both #100 and #60.

.travis.yml

@@ -23,7 +23,7 @@ install:
     - conda info -a
 
     # Replace dep1 dep2 ... with your dependencies
-    - conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib nose
+    - conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib nose sympy
     - source activate test-environment
     #- conda install --yes -c dan_blanchard python-coveralls nose-cov
     - pip install coveralls

README.md

@@ -102,6 +102,7 @@ Requirements/Dependencies
 - scipy 0.15 or above
 - matplotlib 1.4 or above
 - nose 1.3 or above
+- sympy 1.7 or above
 
 Installation
 ------------

commpy/modulation.py

@@ -21,14 +21,14 @@
 
 """
 from bisect import insort
-from itertools import product
 
 import matplotlib.pyplot as plt
-from numpy import arange, array, zeros, pi, cos, sin, sqrt, log2, argmin, \
+from numpy import arange, array, zeros, pi, sqrt, log2, argmin, \
     hstack, repeat, tile, dot, shape, concatenate, exp, \
-    log, vectorize, empty, eye, kron, inf, full, abs, newaxis, minimum, clip
+    log, vectorize, empty, eye, kron, inf, full, abs, newaxis, minimum, clip, fromiter
 from numpy.fft import fft, ifft
 from numpy.linalg import qr, norm
+from sympy.combinatorics.graycode import GrayCode
 
 from commpy.utilities import bitarray2dec, dec2bitarray, signal_power
 
@@ -65,10 +65,16 @@ class Modem:
                         If the constellation is changed to an array-like with length that is not a power of 2.
         """
 
-    def __init__(self, constellation):
+    def __init__(self, constellation, reorder_as_gray=True):
         """ Creates a custom Modem object. """
 
-        self.constellation = constellation
+        if reorder_as_gray:
+            m = log2(len(constellation))
+            gray_code_sequence = GrayCode(m).generate_gray()
+            gray_code_sequence_array = fromiter((int(g, 2) for g in gray_code_sequence), int, len(constellation))
+            self.constellation = array(constellation)[gray_code_sequence_array.argsort()]
+        else:
+            self.constellation = constellation
 
     def modulate(self, input_bits):
         """ Modulate (map) an array of bits to constellation symbols.
@@ -197,10 +203,11 @@ class PSKModem(Modem):
     def __init__(self, m):
         """ Creates a Phase Shift Keying (PSK) Modem object. """
 
-        def _constellation_symbol(i):
-            return cos(2 * pi * (i - 1) / m) + sin(2 * pi * (i - 1) / m) * (0 + 1j)
+        num_bits_symbol = log2(m)
+        if num_bits_symbol != int(num_bits_symbol):
+            raise ValueError('Constellation length must be a power of 2.')
 
-        self.constellation = list(map(_constellation_symbol, arange(m)))
+        super().__init__(exp(1j * arange(0, 2 * pi, 2 * pi / m)))
 
 
 class QAMModem(Modem):
@@ -229,6 +236,7 @@ class QAMModem(Modem):
         ------
         ValueError
                         If the constellation is changed to an array-like with length that is not a power of 2.
+                        If the parameter m would lead to an non-square QAM during initialization.
     """
 
     def __init__(self, m):
@@ -237,16 +245,21 @@ def __init__(self, m):
         Parameters
         ----------
         m : int
-            Size of the QAM constellation.
+            Size of the QAM constellation. Must lead to a square QAM (ie sqrt(m) is an integer).
 
+        Raises
+        ------
+        ValueError
+                        If m would lead to an non-square QAM.
         """
 
-        def _constellation_symbol(i):
-            return (2 * i[0] - 1) + (2 * i[1] - 1) * (1j)
+        num_symb_pam = sqrt(m)
+        if num_symb_pam != int(num_symb_pam):
+            raise ValueError('m must lead to a square QAM.')
 
-        mapping_array = arange(1, sqrt(m) + 1) - (sqrt(m) / 2)
-        self.constellation = list(map(_constellation_symbol,
-                                      list(product(mapping_array, repeat=2))))
+        pam = arange(-num_symb_pam + 1, num_symb_pam, 2)
+        constellation = tile(hstack((pam, pam[::-1])), int(num_symb_pam) // 2) * 1j + pam.repeat(num_symb_pam)
+        super().__init__(constellation)
 
 
 def ofdm_tx(x, nfft, nsc, cp_length):

commpy/tests/test_modulation.py

@@ -3,16 +3,21 @@
 
 from itertools import product
 
-from numpy import zeros, identity, arange, concatenate, log2, array, inf
+from numpy import zeros, identity, arange, concatenate, log2, log10, array, inf, sqrt, sin, pi
 from numpy.random import seed
 from numpy.testing import run_module_suite, assert_allclose, dec, assert_raises, assert_array_equal
+from scipy.special import erf
 
-from commpy.channels import MIMOFlatChannel
+from commpy.channels import MIMOFlatChannel, SISOFlatChannel
 from commpy.links import *
 from commpy.modulation import QAMModem, mimo_ml, bit_lvl_repr, max_log_approx, PSKModem, Modem
 from commpy.utilities import signal_power
 
 
+def Qfunc(x):
+    return 0.5 - 0.5 * erf(x / sqrt(2))
+
+
 @dec.slow
 def test_bit_lvl_repr():
     # Set seed
@@ -124,8 +129,33 @@ def do_custom(self, modem):
         pass
 
 
+@dec.slow
 class TestModulateHardDemodulate(ModemTestcase):
 
+    @staticmethod
+    def check_BER(modem, EbN0dB, BERs_expected):
+        seed(8071996)
+        model = LinkModel(modem.modulate,
+                          SISOFlatChannel(fading_param=(1 + 0j, 0)),
+                          lambda y, _, __, ___: modem.demodulate(y, 'hard'),
+                          modem.num_bits_symbol, modem.constellation, modem.Es)
+        BERs = model.link_performance(EbN0dB + 10 * log10(log2(modem.m)), 5e5, 400, 720)
+        assert_allclose(BERs, BERs_expected, atol=1e-4, rtol=.1,
+                        err_msg='Wrong BER for a standard modulation with {} symbols'.format(modem.m))
+
+    def do_qam(self, modem):
+        EbN0dB = arange(8, 25, 4)
+        nb_symb_pam = sqrt(modem.m)
+        BERs_expected = 2 * (1 - 1 / nb_symb_pam) / log2(nb_symb_pam) * \
+                        Qfunc(sqrt(3 * log2(nb_symb_pam) / (nb_symb_pam ** 2 - 1) * (2 * 10 ** (EbN0dB / 10))))
+        self.check_BER(modem, EbN0dB, BERs_expected)
+
+    def do_psk(self, modem):
+        EbN0dB = arange(15, 25, 4)
+        SERs_expected = 2 * Qfunc(sqrt(2 * modem.num_bits_symbol * 10 ** (EbN0dB / 10)) * sin(pi / modem.m))
+        BERs_expected = SERs_expected / modem.num_bits_symbol
+        self.check_BER(modem, EbN0dB, BERs_expected)
+
     def do(self, modem):
         for bits in product(*((0, 1),) * modem.num_bits_symbol):
             assert_array_equal(bits, modem.demodulate(modem.modulate(bits), 'hard'),

requirements.txt

@@ -2,3 +2,4 @@ numpy>=1.9.2
 scipy>=0.15.0
 matplotlib>=1.4.3
 nose>=1.3.4
+sympy>=1.7.1