Merge branch 'MBF-FD'

Author: BastienTr

README.md

@@ -54,6 +54,7 @@ Available Features
 - OFDM Tx/Rx signal processing
 - MIMO Maximum Likelihood (ML) Detection.
 - MIMO K-best Schnorr-Euchner Detection.
+- MIMO Best-First Detection.
 - Convert channel matrix to Bit-level representation.
 - Computation of LogLikelihood ratio using max-log approximation.
 

commpy/links.py

@@ -59,7 +59,7 @@ def link_performance(link_model, SNRs, send_max, err_min, send_chunk=None, code_
     """
 
     # Initialization
-    BERs = np.empty_like(SNRs, dtype=float)
+    BERs = np.zeros_like(SNRs)
     # Set chunk size and round it to be a multiple of num_bits_symbol*nb_tx to avoid padding
     if send_chunk is None:
         send_chunk = err_min
@@ -83,7 +83,7 @@ def link_performance(link_model, SNRs, send_max, err_min, send_chunk=None, code_
             # Deals with MIMO channel
             if isinstance(link_model.channel, MIMOFlatChannel):
                 nb_symb_vector = len(channel_output)
-                received_msg = np.empty(int(math.ceil(len(msg) / link_model.rate)), int)
+                received_msg = np.empty(int(math.ceil(len(msg) / link_model.rate)))
                 for i in range(nb_symb_vector):
                     received_msg[receive_size * i:receive_size * (i + 1)] = \
                         link_model.receive(channel_output[i], link_model.channel.channel_gains[i],
@@ -101,6 +101,8 @@ def link_performance(link_model, SNRs, send_max, err_min, send_chunk=None, code_
                 bit_err += (msg != link_model.decoder(received_msg)[:len(msg)]).sum()
             bit_send += send_chunk
         BERs[id_SNR] = bit_err / bit_send
+        if bit_err < err_min:
+            break
     return BERs
 
 

commpy/modulation.py

@@ -15,22 +15,25 @@
    ofdm_rx              -- OFDM Receive Signal Processing
    mimo_ml              -- MIMO Maximum Likelihood (ML) Detection.
    kbest                -- MIMO K-best Schnorr-Euchner Detection.
-   bit_lvl_repr         -- Bit level representation.
-   max_log_approx       -- Max-log approximation.
+   best_first_detector  -- MIMO Best-First Detection.
+   bit_lvl_repr         -- Bit Level Representation.
+   max_log_approx       -- Max-Log Approximation.
 
 """
+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, \
     hstack, repeat, tile, dot, shape, concatenate, exp, \
-    log, vectorize, empty, eye, kron, inf
+    log, vectorize, empty, eye, kron, inf, full, abs, newaxis, minimum, clip
 from numpy.fft import fft, ifft
 from numpy.linalg import qr, norm
 
 from commpy.utilities import bitarray2dec, dec2bitarray
 
-__all__ = ['PSKModem', 'QAMModem', 'ofdm_tx', 'ofdm_rx', 'mimo_ml', 'kbest', 'bit_lvl_repr', 'max_log_approx']
+__all__ = ['PSKModem', 'QAMModem', 'ofdm_tx', 'ofdm_rx', 'mimo_ml', 'kbest', 'best_first_detector',
+           'bit_lvl_repr', 'max_log_approx']
 
 
 class Modem:
@@ -346,6 +349,152 @@ def kbest(y, h, constellation, K, noise_var=0, output_type='hard', demode=None):
         raise ValueError('output_type must be "hard" or "soft"')
 
 
+def best_first_detector(y, h, constellation, stack_size, noise_var, demode, llr_max):
+    """ MIMO Best-First Detection.
+
+    Reference: G. He, X. Zhang, et Z. Liang, "Algorithm and Architecture of an Efficient MIMO Detector With Cross-Level
+     Parallel Tree-Search", IEEE Transactions on Very Large Scale Integration (VLSI) Systems, 2019
+
+
+    Parameters
+    ----------
+    y : 1D ndarray
+        Received complex symbols (length: num_receive_antennas)
+
+    h : 2D ndarray
+        Channel Matrix (shape: num_receive_antennas x num_transmit_antennas)
+
+    constellation : 1D ndarray of floats
+        Constellation used to modulate the symbols
+
+    stack_size : tuple of integers
+        Size of each stack (length: num_transmit_antennas - 1)
+
+    noise_var : positive float
+        Noise variance.
+        *Default* value is 0.
+
+    demode : function with prototype binary_word = demode(point)
+        Function that provide the binary word corresponding to a symbol vector.
+
+    llr_max : float
+        Max value for LLR clipping
+
+    Returns
+    -------
+    x : 1D ndarray of Log-Likelihood Ratios.
+        Detected vector (length: num_receive_antennas).
+    """
+
+    class _Node:
+        """ Helper data model that implements __lt__ (aka '<') as required to use bisect.insort. """
+
+        def __init__(self, symb_vectors, partial_metrics):
+            """
+            Recursive initializer that build a sequence of siblings.
+            Inputs are assumed to be ordered based on metric
+            """
+            if len(partial_metrics) == 1:
+                # There is one node to build
+                self.symb_vector = symb_vectors.reshape(-1)  # Insure that self.symb_vector is a 1d-ndarray
+                self.partial_metric = partial_metrics[0]
+                self.best_sibling = None
+            else:
+                # Recursive call to build several nodes
+                self.symb_vector = symb_vectors[:, 0].reshape(-1)  # Insure that self.symb_vector is a 1d-ndarray
+                self.partial_metric = partial_metrics[0]
+                self.best_sibling = _Node(symb_vectors[:, 1:], partial_metrics[1:])
+
+        def __lt__(self, other):
+            return self.partial_metric < other.partial_metric
+
+        def expand(self, yt, r, constellation):
+            """ Build all children and return the best one. constellation must be a numpy ndarray."""
+            # Construct children's symbol vector
+            child_size = self.symb_vector.size + 1
+            children_symb_vectors = empty((child_size, constellation.size), constellation.dtype)
+            children_symb_vectors[1:] = self.symb_vector[:, newaxis]
+            children_symb_vectors[0] = constellation
+
+            # Compute children's partial metric and sort
+            children_metric = abs(yt[-child_size] - r[-child_size, -child_size:].dot(children_symb_vectors)) ** 2
+            children_metric += self.partial_metric
+            ordering = children_metric.argsort()
+
+            # Build children and return the best one
+            return _Node(children_symb_vectors[:, ordering], children_metric[ordering])
+
+    # Extract information from arguments
+    nb_tx, nb_rx = h.shape
+    constellation = array(constellation)
+    m = constellation.size
+    modulation_order = int(log2(m))
+
+    # QR decomposition
+    q, r = qr(h)
+    yt = q.conj().T.dot(y)
+
+    # Initialisation
+    map_metric = inf
+    map_bit_vector = None
+    counter_hyp_metric = full((nb_tx, modulation_order), inf)
+    stacks = tuple([] for _ in range(nb_tx))
+
+    # Start process by adding the best root's child in the last stack
+    stacks[-1].append(_Node(empty(0, constellation.dtype), array(0, float, ndmin=1)).expand(yt, r, constellation))
+
+    # While there is at least one non-empty stack (exempt the first one)
+    while any(stacks[1:]):
+        # Node processing
+        for idx_next_stack in range(len(stacks) - 1):
+            try:
+                idx_this_stack = idx_next_stack + 1
+                best_node = stacks[idx_this_stack].pop(0)
+
+                # Update search radius
+                if map_bit_vector is None:
+                    radius = inf  # No leaf has been reached yet so we keep all nodes
+                else:
+                    bit_vector = demode(best_node.symb_vector).reshape(-1, modulation_order)
+                    bit_vector[bit_vector == 0] = -1
+
+                    # Select the counter hyp metrics that could be affected by this node. Details: eq. (14)-(16) in [1].
+                    try:
+                        a2 = counter_hyp_metric[idx_this_stack:][map_bit_vector[idx_this_stack:] != bit_vector].max()
+                    except ValueError:
+                        a2 = inf  # NumPy cannot compute max on an empty matrix
+                    radius = max(counter_hyp_metric[:idx_this_stack].max(), a2)
+
+                # Process best sibling
+                if best_node.best_sibling is not None and best_node.best_sibling.partial_metric <= radius:
+                    insort(stacks[idx_this_stack], best_node.best_sibling)
+
+                # Process children
+                best_child = best_node.expand(yt, r, constellation)
+                if best_child.partial_metric <= radius:
+                    insort(stacks[idx_next_stack], best_child)
+            except IndexError:  # Raised when popping an empty stack
+                pass
+
+        # LLR update if there is a new leaf
+        if stacks[0]:
+            if stacks[0][0].partial_metric < map_metric:
+                minimum(counter_hyp_metric, map_metric, out=counter_hyp_metric)
+                map_metric = stacks[0][0].partial_metric
+                map_bit_vector = demode(stacks[0][0].symb_vector).reshape(-1, modulation_order)
+                map_bit_vector[map_bit_vector == 0] = -1
+            else:
+                minimum(counter_hyp_metric, stacks[0][0].partial_metric, out=counter_hyp_metric)
+            clip(counter_hyp_metric, map_metric - llr_max, map_metric + llr_max, counter_hyp_metric)
+
+        # Trimming stack according to requested max stack size
+        del stacks[0][0:]  # there is no stack for the leafs
+        for idx_next_stack in range(len(stacks) - 1):
+            del stacks[idx_next_stack + 1][stack_size[idx_next_stack]:]
+
+    return ((map_metric - counter_hyp_metric) * map_bit_vector).reshape(-1)
+
+
 def bit_lvl_repr(H, w):
     """ Bit-level representation of matrix H with weights w.