Several LDPC updates.

 * Add 2 WiMAX LDPC designs, and a link function to model MIMO IDD process.
 * Max-log approximation now return the opposite of the previously LLRs to be consistent with the definition used in the LDPC module.
 * Fix some bugs.

Author: BastienTr

README.md

@@ -73,6 +73,7 @@ Links
 -----
 - Estimate the BER performance of a link model with Monte Carlo simulation.
 - Link model object.
+- Helper function for MIMO Iteration Detection and Decoding scheme.
 
 FAQs
 ----

…coding/designs/ldpc/802.16e/1440.720.txt …elcoding/designs/ldpc/wimax/1440.720.txtcommpy/channelcoding/designs/ldpc/802.16e/1440.720.txt renamed to commpy/channelcoding/designs/ldpc/wimax/1440.720.txt commpy/channelcoding/designs/ldpc/802.16e/1440.720.txt renamed to commpy/channelcoding/designs/ldpc/wimax/1440.720.txt

@@ -11,7 +11,7 @@
 
 def build_matrix(ldpc_code_params):
     """
-    Build the parity check matrix from parameters dictionary and add the result in this dictionary.
+    Build the parity check and generator matrices from parameters dictionary and add the result in this dictionary.
 
     Parameters
     ----------
@@ -28,6 +28,7 @@ def build_matrix(ldpc_code_params):
     ---
     to ldpc_code_param:
             parity_check_matrix (CSC sparse matrix of int8) - parity check matrix.
+            generator_matrix (CSR sparse matrix) - generator matrix of the code.
     """
     n_cnodes = ldpc_code_params['n_cnodes']
     cnode_deg_list = ldpc_code_params['cnode_deg_list']
@@ -37,7 +38,12 @@ def build_matrix(ldpc_code_params):
     for cnode_idx in range(n_cnodes):
         parity_check_matrix[cnode_idx, cnode_adj_list[cnode_idx, :cnode_deg_list[cnode_idx]]] = 1
 
-    ldpc_code_params['parity_check_matrix'] = parity_check_matrix.tocsc()
+    parity_check_matrix = parity_check_matrix.tocsc()
+    systematic_part = parity_check_matrix[:, -n_cnodes:]
+    parity_part = parity_check_matrix[:, :-n_cnodes]
+
+    ldpc_code_params['parity_check_matrix'] = parity_check_matrix
+    ldpc_code_params['generator_matrix'] = splg.inv(systematic_part).dot(parity_part).tocsr()
 
 
 def get_ldpc_code_params(ldpc_design_filename, compute_matrix=False):
@@ -170,7 +176,7 @@ def ldpc_bp_decode(llr_vec, ldpc_code_params, decoder_algorithm, n_iters):
     Parameters
     ----------
     llr_vec : 1D array of float
-        Received codeword LLR values from the channel.
+        Received codeword LLR values from the channel. They will be clipped in [-38, 38].
 
     ldpc_code_params : dictionary that at least contains these parameters
         Parameters of the LDPC code as provided by `get_ldpc_code_params`:
@@ -204,6 +210,9 @@ def ldpc_bp_decode(llr_vec, ldpc_code_params, decoder_algorithm, n_iters):
         LLR values corresponding to the decoded output.
     """
 
+    # Clip LLRs into [-38, 38]
+    llr_vec.clip(-38, 38, llr_vec)
+
     n_cnodes = ldpc_code_params['n_cnodes']
     n_vnodes = ldpc_code_params['n_vnodes']
     max_cnode_deg = ldpc_code_params['max_cnode_deg']
@@ -256,7 +265,7 @@ def ldpc_bp_decode(llr_vec, ldpc_code_params, decoder_algorithm, n_iters):
             msg_sum = np.sum(cnode_list_msgs)
 
             # Compute messages on outgoing edges using the incoming message sum (LLRs are clipped in [-38, 38])
-            np.clip(llr_vec[vnode_idx] + msg_sum - cnode_list_msgs, -38, 38, vnode_msgs[start_idx:start_idx+offset])
+            vnode_msgs[start_idx:start_idx+offset] = llr_vec[vnode_idx] + msg_sum - cnode_list_msgs
 
             # Update output LLRs and decoded word
             out_llrs[vnode_idx] = llr_vec[vnode_idx] + msg_sum
@@ -329,8 +338,8 @@ def write_ldpc_params(parity_check_matrix, file_path):
 
 def triang_ldpc_systematic_encode(message_bits, ldpc_code_params, pad=True):
     """
-    Encode bits using the LDPC code specified. If the parity check matrix and/or the generator matrix are not computed,
-    this function will build the missing one(s) and add them to the dictionary.
+    Encode bits using the LDPC code specified. If the  generator matrix is not computed, this function will build it
+    and add it to the dictionary. It will also add the parity check matrix.
 
     This function work only for LDPC specified by a triangular parity check matrix.
 
@@ -339,27 +348,25 @@ def triang_ldpc_systematic_encode(message_bits, ldpc_code_params, pad=True):
     message_bits : 1D-array
         Message bit to encode.
 
-    ldpc_code_params : dictionary that at least contains on of these options:
+    ldpc_code_params : dictionary that at least contains one of these options:
         Option 1: generator matrix is available.
-            generator_matrix (2D-array or sparse matrix) - generator matrix of the code.
-        Option 2: parity check matrix is available, the generator matrix will be added as a CSR sparse matrix.
-            parity_check_matrix (sparse matrix) - parity check matrix of the code.
-        Option 3: generator and parity check matrices will be added as sparse matrices of integers.
-            n_vnodes (int) - number of variable nodes.
-            n_cnodes (int) - number of check nodes.
-            max_cnode_deg (int) - maximal degree of a check node.
-            cnode_adj_list (1D-ndarray of ints) - flatten array so that
-                cnode_adj_list.reshape((n_cnodes, max_cnode_deg)) gives for each check node the adjacent variable nodes.
-            cnode_deg_list (1D-ndarray of ints) - degree of each check node.
+                generator_matrix (2D-array or sparse matrix) - generator matrix of the code.
+        Option 2: generator and parity check matrices will be added as sparse matrices.
+                n_vnodes (int) - number of variable nodes.
+                n_cnodes (int) - number of check nodes.
+                max_cnode_deg (int) - maximal degree of a check node.
+                cnode_adj_list (1D-ndarray of ints) - flatten array so that
+                    cnode_adj_list.reshape((n_cnodes, max_cnode_deg)) gives for each check node the adjacent variable nodes.
+                cnode_deg_list (1D-ndarray of ints) - degree of each check node.
 
     pad : boolean
-        Whether to add '0' padding to the message to fit the block length.
-        *Default* is True.
+            Whether to add '0' padding to the message to fit the block length.
+            *Default* is True.
 
     Returns
     -------
-    coded_message : 1D-ndarray or 2D-ndarray depending on the number of blocks
-        Coded message with the systematic part at the beginning.
+    coded_message : 1D-ndarray or 2D-ndarray of int8 depending on the number of blocks
+                    Coded message with the systematic part at the beginning.
 
     Raises
     ------
@@ -368,15 +375,7 @@ def triang_ldpc_systematic_encode(message_bits, ldpc_code_params, pad=True):
     """
 
     if ldpc_code_params.get('generator_matrix') is None:
-        if ldpc_code_params.get('parity_check_matrix') is None:
-            build_matrix(ldpc_code_params)
-
-        parity_check_matrix = ldpc_code_params['parity_check_matrix']
-        block_length = parity_check_matrix.shape[0]
-
-        systematic_part = parity_check_matrix[:, -block_length:]
-        parity_part = parity_check_matrix[:, :-block_length]
-        ldpc_code_params['generator_matrix'] = splg.inv(systematic_part).dot(parity_part).tocsr()
+        build_matrix(ldpc_code_params)
 
     block_length = ldpc_code_params['generator_matrix'].shape[1]
     modulo = len(message_bits) % block_length
@@ -388,4 +387,4 @@ def triang_ldpc_systematic_encode(message_bits, ldpc_code_params, pad=True):
     message_bits = message_bits.reshape(block_length, -1, order='F')
 
     parity_part = ldpc_code_params['generator_matrix'].dot(message_bits) % 2
-    return np.vstack((message_bits, parity_part)).squeeze()
+    return np.vstack((message_bits, parity_part)).squeeze().astype(np.int8)

…oding/designs/ldpc/802.16e/960.720.a.txt …lcoding/designs/ldpc/wimax/960.720.a.txtcommpy/channelcoding/designs/ldpc/802.16e/960.720.a.txt renamed to commpy/channelcoding/designs/ldpc/wimax/960.720.a.txt commpy/channelcoding/designs/ldpc/802.16e/960.720.a.txt renamed to commpy/channelcoding/designs/ldpc/wimax/960.720.a.txt

@@ -65,16 +65,16 @@ def test_ldpc_bp_decode(self):
 
     def test_write_ldpc_params(self):
         with TemporaryDirectory() as tmp_dir:
-            parity_check_matrix = choice((0, 1), (1440, 720))
+            parity_check_matrix = choice((0, 1), (720, 1440))
 
             file_path = tmp_dir + '/matrix.txt'
             write_ldpc_params(parity_check_matrix, file_path)
             assert_equal(get_ldpc_code_params(file_path, True)['parity_check_matrix'].toarray(), parity_check_matrix,
                          'The loaded matrix is not equal to the written one.')
 
     def test_triang_ldpc_systematic_encode(self):
-        ldpc_design_files = (os.path.join(self.dir, '../designs/ldpc/802.16e/1440.720.txt'),
-                             os.path.join(self.dir, '../designs/ldpc/802.16e/960.720.a.txt'))
+        ldpc_design_files = (os.path.join(self.dir, '../designs/ldpc/wimax/1440.720.txt'),
+                             os.path.join(self.dir, '../designs/ldpc/wimax/960.720.a.txt'))
         wimax_ldpc_params = [get_ldpc_code_params(ldpc_design_file) for ldpc_design_file in ldpc_design_files]
 
         for param in wimax_ldpc_params:

commpy/channelcoding/ldpc.py

@@ -11,16 +11,18 @@
 
    link_performance     -- Estimate the BER performance of a link model with Monte Carlo simulation.
    LinkModel            -- Link model object.
+   idd_decoder          -- Produce the decoder function to model a MIMO IDD decoder.
 """
 from __future__ import division  # Python 2 compatibility
 
 import math
+from inspect import getfullargspec
 
 import numpy as np
 
 from commpy.channels import MIMOFlatChannel
 
-__all__ = ['link_performance', 'LinkModel']
+__all__ = ['link_performance', 'LinkModel', 'idd_decoder']
 
 
 def link_performance(link_model, SNRs, send_max, err_min, send_chunk=None, code_rate=1):
@@ -65,6 +67,7 @@ def link_performance(link_model, SNRs, send_max, err_min, send_chunk=None, code_
     send_chunk = max(divider, send_chunk // divider * divider)
 
     receive_size = link_model.channel.nb_tx * link_model.num_bits_symbol
+    full_args_decoder = len(getfullargspec(link_model.decoder).args) > 1
 
     # Computations
     for id_SNR in range(len(SNRs)):
@@ -89,68 +92,88 @@ def link_performance(link_model, SNRs, send_max, err_min, send_chunk=None, code_
                 received_msg = link_model.receive(channel_output, link_model.channel.channel_gains,
                                                   link_model.constellation, link_model.channel.noise_std ** 2)
             # Count errors
-            bit_err += (msg != link_model.decoder(received_msg)).sum()
+            if full_args_decoder:
+                decoded_bits = link_model.decoder(channel_output, link_model.channel.channel_gains,
+                                                  link_model.constellation, link_model.channel.noise_std ** 2,
+                                                  received_msg, link_model.channel.nb_tx * link_model.num_bits_symbol)
+                bit_err += (msg != decoded_bits[:len(msg)]).sum()
+            else:
+                bit_err += (msg != link_model.decoder(received_msg)[:len(msg)]).sum()
             bit_send += send_chunk
         BERs[id_SNR] = bit_err / bit_send
     return BERs
 
 
 class LinkModel:
     """
-        Construct a link model.
+    Construct a link model.
+
+    Parameters
+    ----------
+    modulate : function with same prototype as Modem.modulate
 
-        Parameters
-        ----------
-        modulate : function with same prototype as Modem.modulate
+    channel : FlatChannel object
 
-        channel : _FlatChannel object
+    receive : function with prototype receive(y, H, constellation, noise_var) that return a binary array.
+                y : 1D ndarray
+                    Received complex symbols (shape: num_receive_antennas x 1)
 
-        receive : function with prototype receive(y, H, constellation, noise_var) that return a binary array.
-                    y : 1D ndarray
-                        Received complex symbols (shape: num_receive_antennas x 1)
+                h : 2D ndarray
+                    Channel Matrix (shape: num_receive_antennas x num_transmit_antennas)
 
-                    h : 2D ndarray
-                        Channel Matrix (shape: num_receive_antennas x num_transmit_antennas)
+                constellation : 1D ndarray
 
-                    constellation : 1D ndarray
+                noise_var : positive float
+                            Noise variance
 
-                    noise_var : positive float
-                                Noise variance
+    num_bits_symbols : int
 
-        num_bits_symbols : int
+    constellation : array of float or complex
 
-        constellation : array of float or complex
+    Es : float
+         Average energy per symbols.
+         *Default* Es=1.
 
-        Es : float
-             Average energy per symbols.
-             *Default* Es=1.
+    decoder : function with prototype decoder(array) or decoder(y, H, constellation, noise_var, array) that return a
+                binary array.
+              *Default* is no process.
 
-        decoder : function with prototype decoder(binary array) that return a binary array.
-                  *Default* is no process.
+    rate : float
+        Code rate.
+        *Default* is 1.
 
-        Attributes
-        ----------
-        modulate : function with same prototype as Modem.modulate
+    Attributes
+    ----------
+    modulate : function with same prototype as Modem.modulate
+
+    channel : _FlatChannel object
 
-        channel : _FlatChannel object
+    receive : function with prototype receive(y, H, constellation, noise_var) that return a binary array.
+                y : 1D ndarray
+                    Received complex symbols (shape: num_receive_antennas x 1)
 
-        receive : function with prototype receive(y, H, constellation) that return a binary array.
-                    y : 1D ndarray of floats
-                        Received complex symbols (shape: num_receive_antennas x 1)
+                h : 2D ndarray
+                    Channel Matrix (shape: num_receive_antennas x num_transmit_antennas)
 
-                    h : 2D ndarray of floats
-                        Channel Matrix (shape: num_receive_antennas x num_transmit_antennas)
+                constellation : 1D ndarray
 
-                    constellation : 1D ndarray of floats
+                noise_var : positive float
+                            Noise variance
 
-        num_bits_symbols : int
+    num_bits_symbols : int
 
-        constellation : array of float or complex
+    constellation : array of float or complex
 
-        Es : float
-             Average energy per symbols.
-             *Default* Es=1.
-        """
+    Es : float
+         Average energy per symbols.
+
+    decoder : function with prototype decoder(binary array) that return a binary array.
+              *Default* is no process.
+
+    rate : float
+        Code rate.
+        *Default* is 1.
+    """
 
     def __init__(self, modulate, channel, receive, num_bits_symbol, constellation, Es=1, decoder=None, rate=1):
         self.modulate = modulate
@@ -165,3 +188,65 @@ def __init__(self, modulate, channel, receive, num_bits_symbol, constellation, E
             self.decoder = lambda msg: msg
         else:
             self.decoder = decoder
+
+
+def idd_decoder(word_size, detector, decoder, n_it):
+    """
+    Produce a decoder function that model the specified MIMO iterative detection and decoding (IDD) process.
+    The returned function can be used as is to build a working LinkModel object.
+
+    Parameters
+    ----------
+    word_size : positive integer
+                Size of the words exchanged between the detector and the decoder.
+
+    detector : function with prototype detector(y, H, constellation, noise_var, a_priori) that return a LLRs array.
+                y : 1D ndarray
+                    Received complex symbols (shape: num_receive_antennas x 1).
+
+                h : 2D ndarray
+                    Channel Matrix (shape: num_receive_antennas x num_transmit_antennas).
+
+                constellation : 1D ndarray.
+
+                noise_var : positive float
+                            Noise variance.
+
+                a_priori : 1D ndarray of floats
+                            A priori as Log-Likelihood Ratios.
+
+    decoder : function with the same signature as detector.
+
+    n_it : positive integer
+            Number or iteration during the IDD process.
+
+    Returns
+    -------
+    decode : function useable as it is to build a LinkModel object that produce a bit array from the parameters
+                y : 1D ndarray
+                    Received complex symbols (shape: num_receive_antennas x 1).
+
+                h : 2D ndarray
+                    Channel Matrix (shape: num_receive_antennas x num_transmit_antennas).
+
+                constellation : 1D ndarray
+
+                noise_var : positive float
+                            Noise variance.
+
+                bits_per_send : positive integer
+                                Number or bit send at each symbol vector.
+    """
+    def decode(y, h, constellation, noise_var, a_priori, bits_per_send):
+        a_priori_decoder = a_priori
+        nb_vect, nb_rx, nb_tx = h.shape
+        for iteration in range(n_it):
+            a_priori_detector = (decoder(a_priori_decoder) - a_priori_decoder)
+            for i in range(nb_vect):
+                a_priori_decoder[i * bits_per_send:(i + 1) * bits_per_send] = \
+                    detector(y[i], h[i], constellation, noise_var,
+                             a_priori_detector[i * bits_per_send:(i + 1) * bits_per_send])
+            a_priori_decoder -= a_priori_detector
+        return np.signbit(a_priori_decoder + a_priori_detector)
+
+    return decode

commpy/channelcoding/tests/test_ldpc.py

@@ -419,4 +419,4 @@ def max_log_approx(y, h, noise_var, pts_list, demode):
         # Compute LLR
         LLR[k] = min(norms0) - min(norms1)
 
-    return LLR / (2 * noise_var)
+    return -LLR / (2 * noise_var)