Provide a tools to estimate the BER performance of a link model with Monte Carlo simulation.

This commit provides:
 * the tool (`commpy.utilities.link_performance`),
 * a test that check both SISO and MIMO channels,
 * some typo corrections,
 * an update of the documentation.

Author: BastienTr

README.md

@@ -66,6 +66,7 @@ Utilities
 - Hamming distance, Euclidean distance.
 - Upsample
 - Power of a discrete-time signal
+- Estimate the BER performance of a link model with Monte Carlo simulation.
 
 FAQs
 ----

commpy/channels.py

@@ -108,10 +108,11 @@ class SISOFlatChannel(_FlatChannel):
     ----------
     noise_std    : float, optional
                    Noise standard deviation.
-                   Default value is None and then the value must set later.
+                   *Default* value is None and then the value must set later.
 
     fading_param : tuple of 2 floats, optional
-                   Parameters of the fading (see attribute for details). Default value is (1,0) i.e. no fading.
+                   Parameters of the fading (see attribute for details).
+                   *Default* value is (1,0) i.e. no fading.
 
     Attributes
     ----------
@@ -255,11 +256,11 @@ class MIMOFlatChannel(_FlatChannel):
 
     noise_std    : float, optional
                    Noise standard deviation.
-                   Default value is None and then the value must set later.
+                   *Default* value is None and then the value must set later.
 
     fading_param : tuple of 3 floats, optional
                    Parameters of the fading. The complete tuple must be set each time.
-                   Default value is (zeros((nb_rx, nb_tx)), identity(nb_tx), identity(nb_rx)) i.e. Rayleigh fading.
+                   *Default* value is (zeros((nb_rx, nb_tx)), identity(nb_tx), identity(nb_rx)) i.e. Rayleigh fading.
 
     Attributes
     ----------

commpy/modulation.py

@@ -1,4 +1,3 @@
-
 # Authors: Veeresh Taranalli <[email protected]>
 # License: BSD 3-Clause
 
@@ -48,17 +47,15 @@ def modulate(self, input_bits):
 
         """
         mapfunc = vectorize(lambda i:
-            self.constellation[bitarray2dec(input_bits[i:i+self.num_bits_symbol])])
+                            self.constellation[bitarray2dec(input_bits[i:i + self.num_bits_symbol])])
 
         baseband_symbols = mapfunc(arange(0, len(input_bits), self.num_bits_symbol))
 
         return baseband_symbols
 
-    def demodulate(self, input_symbols, demod_type, noise_var = 0):
+    def demodulate(self, input_symbols, demod_type, noise_var=0):
         """ Demodulate (map) a set of constellation symbols to corresponding bits.
 
-        Supports hard-decision demodulation only.
-
         Parameters
         ----------
         input_symbols : 1D ndarray of complex floats
@@ -78,10 +75,10 @@ def demodulate(self, input_symbols, demod_type, noise_var = 0):
 
         """
         if demod_type == 'hard':
-            index_list = map(lambda i: argmin(abs(input_symbols[i] - self.constellation)), \
+            index_list = map(lambda i: argmin(abs(input_symbols[i] - self.constellation)),
                              range(0, len(input_symbols)))
             demod_bits = hstack(map(lambda i: dec2bitarray(i, self.num_bits_symbol),
-                                index_list))
+                                    index_list))
         elif demod_type == 'soft':
             demod_bits = zeros(len(input_symbols) * self.num_bits_symbol)
             for i in arange(len(input_symbols)):
@@ -91,10 +88,12 @@ def demodulate(self, input_symbols, demod_type, noise_var = 0):
                     llr_den = 0
                     for const_index in self.symbol_mapping:
                         if (const_index >> bit_index) & 1:
-                            llr_num = llr_num + exp((-abs(current_symbol - self.constellation[const_index])**2)/noise_var)
+                            llr_num = llr_num + exp(
+                                (-abs(current_symbol - self.constellation[const_index]) ** 2) / noise_var)
                         else:
-                            llr_den = llr_den + exp((-abs(current_symbol - self.constellation[const_index])**2)/noise_var)
-                    demod_bits[i*self.num_bits_symbol + self.num_bits_symbol - 1 - bit_index] = log(llr_num/llr_den)
+                            llr_den = llr_den + exp(
+                                (-abs(current_symbol - self.constellation[const_index]) ** 2) / noise_var)
+                    demod_bits[i * self.num_bits_symbol + self.num_bits_symbol - 1 - bit_index] = log(llr_num / llr_den)
         else:
             pass
             # throw an error
@@ -105,8 +104,10 @@ def demodulate(self, input_symbols, demod_type, noise_var = 0):
 class PSKModem(Modem):
     """ Creates a Phase Shift Keying (PSK) Modem object. """
 
+    Es = 1
+
     def _constellation_symbol(self, i):
-        return cos(2*pi*(i-1)/self.m) + sin(2*pi*(i-1)/self.m)*(0+1j)
+        return cos(2 * pi * (i - 1) / self.m) + sin(2 * pi * (i - 1) / self.m) * (0 + 1j)
 
     def __init__(self, m):
         """ Creates a Phase Shift Keying (PSK) Modem object.
@@ -121,13 +122,13 @@ def __init__(self, m):
         self.num_bits_symbol = int(log2(self.m))
         self.symbol_mapping = arange(self.m)
         self.constellation = list(map(self._constellation_symbol,
-                                 self.symbol_mapping))
+                                      self.symbol_mapping))
 
 class QAMModem(Modem):
     """ Creates a Quadrature Amplitude Modulation (QAM) Modem object."""
 
     def _constellation_symbol(self, i):
-        return (2*i[0]-1) + (2*i[1]-1)*(1j)
+        return (2 * i[0] - 1) + (2 * i[1] - 1) * (1j)
 
     def __init__(self, m):
         """ Creates a Quadrature Amplitude Modulation (QAM) Modem object.
@@ -142,9 +143,11 @@ def __init__(self, m):
         self.m = m
         self.num_bits_symbol = int(log2(self.m))
         self.symbol_mapping = arange(self.m)
-        mapping_array = arange(1, sqrt(self.m)+1) - (sqrt(self.m)/2)
+        mapping_array = arange(1, sqrt(self.m) + 1) - (sqrt(self.m) / 2)
         self.constellation = list(map(self._constellation_symbol,
-                                 list(product(mapping_array, repeat=2))))
+                                      list(product(mapping_array, repeat=2))))
+        self.Es = 2 * (self.m - 1) / 3
+
 
 def ofdm_tx(x, nfft, nsc, cp_length):
     """ OFDM Transmit signal generation """
@@ -155,29 +158,31 @@ def ofdm_tx(x, nfft, nsc, cp_length):
     ofdm_tx_signal = array([])
 
     for i in range(0, shape(x)[1]):
-        symbols = x[:,i]
+        symbols = x[:, i]
         ofdm_sym_freq = zeros(nfft, dtype=complex)
-        ofdm_sym_freq[1:(nsc/2)+1] = symbols[nsc/2:]
-        ofdm_sym_freq[-(nsc/2):] = symbols[0:nsc/2]
+        ofdm_sym_freq[1:(nsc / 2) + 1] = symbols[nsc / 2:]
+        ofdm_sym_freq[-(nsc / 2):] = symbols[0:nsc / 2]
         ofdm_sym_time = ifft(ofdm_sym_freq)
         cp = ofdm_sym_time[-cp_length:]
         ofdm_tx_signal = concatenate((ofdm_tx_signal, cp, ofdm_sym_time))
 
     return ofdm_tx_signal
 
+
 def ofdm_rx(y, nfft, nsc, cp_length):
     """ OFDM Receive Signal Processing """
 
-    num_ofdm_symbols = int(len(y)/(nfft + cp_length))
+    num_ofdm_symbols = int(len(y) / (nfft + cp_length))
     x_hat = zeros([nsc, num_ofdm_symbols], dtype=complex)
 
     for i in range(0, num_ofdm_symbols):
-        ofdm_symbol = y[i*nfft + (i+1)*cp_length:(i+1)*(nfft + cp_length)]
+        ofdm_symbol = y[i * nfft + (i + 1) * cp_length:(i + 1) * (nfft + cp_length)]
         symbols_freq = fft(ofdm_symbol)
-        x_hat[:,i] = concatenate((symbols_freq[-nsc/2:], symbols_freq[1:(nsc/2)+1]))
+        x_hat[:, i] = concatenate((symbols_freq[-nsc / 2:], symbols_freq[1:(nsc / 2) + 1]))
 
     return x_hat
 
+
 def mimo_ml(y, h, constellation):
     """ MIMO ML Detection.
 
@@ -195,7 +200,7 @@ def mimo_ml(y, h, constellation):
     """
     m = len(constellation)
     x_ideal = array([tile(constellation, m), repeat(constellation, m)])
-    y_vector = tile(y, m*m)
+    y_vector = tile(y, m * m)
     min_idx = argmin(sum(abs(y_vector - dot(h, x_ideal)), axis=0))
     x_r = x_ideal[:, min_idx]
 

commpy/tests/test_channels.py

@@ -277,7 +277,7 @@ 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=1e-2,
+            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')

commpy/tests/test_utilities.py

@@ -0,0 +1,39 @@
+# Authors: Bastien Trotobas <[email protected]>
+# License: BSD 3-Clause
+
+from __future__ import division  # Python 2 compatibility
+
+from numpy import arange, sqrt, identity, zeros, log10
+from numpy.random import seed
+from numpy.testing import run_module_suite, assert_allclose, dec
+from scipy.special import erfc
+
+from commpy.channels import MIMOFlatChannel, SISOFlatChannel
+from commpy.modulation import QAMModem, kbest
+from commpy.utilities import link_performance
+
+
+@dec.slow
+def test_link_performance():
+    # Apply link_performance to SISO QPSK and AWGN channel
+    BERs = link_performance(QAMModem(4), SISOFlatChannel(fading_param=(1 + 0j, 0)), None,
+                            range(0, 9, 2), 600e4, 600)
+    desired = erfc(sqrt(10**(arange(0, 9, 2) / 10) / 2)) / 2
+    assert_allclose(BERs, desired, rtol=0.25, err_msg='Wrong performance for SISO QPSK and AWGN channel')
+
+    # Apply link_performance to MIMO 16QAM and 4x4 Rayleigh channel
+    RayleighChannel = MIMOFlatChannel(4, 4, None, (zeros((4, 4), dtype=complex), identity(4), identity(4)))
+    SNRs = range(0, 21, 5) + 10 * log10(4)
+
+    def kbest16(y, h, constellation):
+        return kbest(y, h, constellation, 16)
+
+    BERs = link_performance(QAMModem(16), RayleighChannel, kbest16,
+                            SNRs, 600e4, 600)
+    desired = (2e-1, 1e-1, 3e-2, 2e-3, 4e-5)  # From reference
+    assert_allclose(BERs, desired, rtol=1, err_msg='Wrong performance for MIMO 16QAM and 4x4 Rayleigh channel')
+
+
+if __name__ == "__main__":
+    seed(17121996)
+    run_module_suite()

commpy/utilities.py

@@ -15,12 +15,17 @@
    euclid_dist          -- Squared Euclidean distance.
    upsample             -- Upsample by an integral factor (zero insertion).
    signal_power         -- Compute the power of a discrete time signal.
-
+   link_performance     -- Estimate the BER performance of a link model with Monte Carlo simulation.
 
 """
+from __future__ import division  # Python 2 compatibility
+
 import numpy as np
 
-__all__ = ['dec2bitarray', 'bitarray2dec', 'hamming_dist', 'euclid_dist', 'upsample', 'signal_power']
+from commpy.channels import MIMOFlatChannel
+
+__all__ = ['dec2bitarray', 'bitarray2dec', 'hamming_dist', 'euclid_dist', 'upsample',
+           'signal_power', 'link_performance']
 
 
 def dec2bitarray(in_number, bit_width):
@@ -165,7 +170,84 @@ def signal_power(signal):
 
     @np.vectorize
     def square_abs(s):
-        return abs(s)**2
+        return abs(s) ** 2
 
     P = np.mean(square_abs(signal))
     return P
+
+
+def link_performance(modem, channel, detector, SNRs, send_max, err_min, send_chunck=None, code_rate=1):
+    """
+    Estimate the BER performance of a link model with Monte Carlo simulation.
+
+    Parameters
+    ----------
+    modem : modem object
+            Modem used to modulate and demodulate.
+
+    channel : channel object
+              Channel through which the message is propagated.
+
+    detector : function with prototype detector(y, h, constellation) or None
+               Detector to decode channel output. See detectors in commpy.modulation for details on the prototype.
+
+    SNRs : 1D arraylike
+           Signal to Noise ratio in dB defined as :math:`SNR_{dB} = (E_b/N_0)_{dB} + 10 \log_{10}(R_cM_c)`
+           where :math:`Rc` is the code rate and :math:`Mc` the modulation rate.
+
+    send_max : int
+               Maximum number of bits send for each SNR.
+
+    err_min : int
+              link_performance send bits until it reach err_min errors (see also send_max).
+
+    send_chunck : int
+                  Number of bits to be send at each iteration.
+                  *Default*: send_chunck = err_min
+
+    code_rate : float in (0,1]
+                Rate of the used code.
+                *Default*: 1 i.e. no code.
+    """
+
+    # Initialization
+    BERs = np.empty_like(SNRs, dtype=float)
+
+    # Handles the case detector is None
+    if detector is None:
+        def detector(y, H, constellation):
+            return y
+
+    # Set chunck size and round it to be a multiple of num_bits_symbol*nb_tx to avoid padding
+    if send_chunck is None:
+        send_chunck = err_min
+    divider = modem.num_bits_symbol * channel.nb_tx
+    send_chunck = max(divider, send_chunck // divider * divider)
+
+    # Computations
+    for id_SNR in range(len(SNRs)):
+        channel.set_SNR_dB(SNRs[id_SNR], code_rate, modem.Es)
+        bit_send = 0
+        bit_err = 0
+        while bit_send < send_max and bit_err < err_min:
+            # Propagate some bits
+            msg = np.random.choice((0, 1), send_chunck)
+            symbs = modem.modulate(msg)
+            channel_output = channel.propagate(symbs)
+
+            # Deals with MIMO channel
+            if isinstance(channel, MIMOFlatChannel):
+                nb_symb_vector = len(channel_output)
+                detected_msg = np.empty(nb_symb_vector * channel.nb_tx, dtype=channel_output.dtype)
+                for i in range(nb_symb_vector):
+                    detected_msg[channel.nb_tx * i:channel.nb_tx * (i+1)] = \
+                        detector(channel_output[i], channel.channel_gains[i], modem.constellation)
+            else:
+                detected_msg = channel_output
+
+            # Count errors
+            received_msg = modem.demodulate(detected_msg, 'hard')
+            bit_err += (msg != received_msg[:len(msg)]).sum()  # Remove MIMO padding
+            bit_send += send_chunck
+        BERs[id_SNR] = bit_err / bit_send
+    return BERs