Soft decoding for convolutionnal codes using log-likelihood as metrics.

Author: BastienTr

commpy/channelcoding/convcode.py

@@ -99,8 +99,8 @@ class Trellis:
      [2 1]]
     References
     ----------
-    [1]  S. Benedetto, R. Garello, et G. Montorsi, "A search for good convolutional codes to be used in the
-    construction of turbo codes", IEEE Transactions on Communications, vol. 46, nᵒ 9, p. 1101‑1105, sept. 1998.
+    [1] S. Benedetto, R. Garello et G. Montorsi, "A search for good convolutional codes to be used in the
+    construction of turbo codes", IEEE Transactions on Communications, vol. 46, n. 9, p. 1101-1005, spet. 1998
     """
     def __init__(self, memory, g_matrix, feedback = None, code_type = 'default'):
 
@@ -389,7 +389,7 @@ def conv_encode(message_bits, trellis, termination = 'term', puncture_matrix=Non
 
     outbits = np.zeros(number_outbits, 'int')
     if puncture_matrix is not None:
-        p_outbits = np.zeros(number_outbits)
+        p_outbits = np.zeros(number_outbits, 'int')
     else:
         p_outbits = np.zeros(int(number_outbits*
             puncture_matrix[0:].sum()/np.size(puncture_matrix, 1)), 'int')
@@ -477,7 +477,9 @@ def _acs_traceback(r_codeword, trellis, decoding_type,
             if decoding_type == 'hard':
                 branch_metric = hamming_dist(r_codeword.astype(int), i_codeword_array.astype(int))
             elif decoding_type == 'soft':
-                pass
+                neg_LL_0 = np.log(np.exp(r_codeword) + 1)  # negative log-likelihood to have received a 0
+                neg_LL_1 = neg_LL_0 - r_codeword  # negative log-likelihood to have received a 1
+                branch_metric = np.where(i_codeword_array, neg_LL_1, neg_LL_0).sum()
             elif decoding_type == 'unquantized':
                 i_codeword_array = 2*i_codeword_array - 1
                 branch_metric = euclid_dist(r_codeword, i_codeword_array)
@@ -531,7 +533,7 @@ def viterbi_decode(coded_bits, trellis, tb_depth=None, decoding_type='hard'):
     decoding_type : str {'hard', 'soft', 'unquantized'}
         The type of decoding to be used.
         'hard' option is used for hard inputs (bits) to the decoder, e.g., BSC channel.
-        'soft' option is used for soft inputs (LLRs) to the decoder.
+        'soft' option is used for soft inputs (LLRs) to the decoder. LLRs are clipped in [-500, 500].
         'unquantized' option is used for soft inputs (real numbers) to the decoder, e.g., BAWGN channel.
     Returns
     -------
@@ -561,26 +563,29 @@ def viterbi_decode(coded_bits, trellis, tb_depth=None, decoding_type='hard'):
 
     path_metrics = np.full((trellis.number_states, 2), np.inf)
     path_metrics[0][0] = 0
-    paths = np.full((trellis.number_states, tb_depth), np.iinfo(int).max, 'int')
+    paths = np.empty((trellis.number_states, tb_depth), 'int')
     paths[0][0] = 0
 
     decoded_symbols = np.zeros([trellis.number_states, tb_depth], 'int')
-    decoded_bits = np.empty(math.ceil(L / k) * k + tb_depth, 'int')
+    decoded_bits = np.empty(int(math.ceil((L + tb_depth) / k) * k), 'int')
     r_codeword = np.zeros(n, 'int')
 
     tb_count = 1
     count = 0
     current_number_states = trellis.number_states
 
+    coded_bits = coded_bits.clip(-500, 500)
+
     for t in range(1, int((L+total_memory)/k)):
         # Get the received codeword corresponding to t
         if t <= L // k:
             r_codeword = coded_bits[(t-1)*n:t*n]
+        # Pad with '0'
         else:
             if decoding_type == 'hard':
                 r_codeword[:] = 0
             elif decoding_type == 'soft':
-                pass
+                r_codeword[:] = np.iinfo(int).min
             elif decoding_type == 'unquantized':
                 r_codeword[:] = -1
             else:

commpy/channelcoding/tests/test_convcode.py

@@ -1,9 +1,9 @@
 # Authors: Veeresh Taranalli <veeresht@gmail.com>
 # License: BSD 3-Clause
 
-from numpy import array
-from numpy.random import randint
-from numpy.testing import assert_array_equal, dec
+from numpy import array, inf
+from numpy.random import randint, randn
+from numpy.testing import assert_array_equal, dec, assert_array_less
 
 from commpy.channelcoding.convcode import Trellis, conv_encode, viterbi_decode
 
@@ -122,3 +122,13 @@ def test_conv_encode_viterbi_decode(self):
                 coded_syms = 2.0 * coded_bits - 1
                 decoded_bits = viterbi_decode(coded_syms, self.trellis[i], 15, 'unquantized')
                 assert_array_equal(decoded_bits[:len(msg)], msg)
+
+                coded_bits = conv_encode(msg, self.trellis[i])
+                coded_syms = 10.0 * coded_bits - 5 + randn(len(coded_bits)) * 2
+                decoded_bits = viterbi_decode(coded_syms, self.trellis[i], 15, 'soft')
+                assert_array_less((decoded_bits[:len(msg)] - msg).sum(), 0.03 * blocklength)
+
+                coded_bits = conv_encode(msg, self.trellis[i])
+                coded_syms = (2.0 * coded_bits - 1) * inf
+                decoded_bits = viterbi_decode(coded_syms, self.trellis[i], 15, 'soft')
+                assert_array_less((decoded_bits[:len(msg)] - msg).sum(), 0.03 * blocklength)

commpy/links.py

@@ -61,12 +61,20 @@ def link_performance(link_model, SNRs, send_max, err_min, send_chunk=None, code_
     divider = link_model.num_bits_symbol * link_model.channel.nb_tx
     send_chunk = max(divider, send_chunk // divider * divider)
 
+    # Evaluate the size of each reception
+    msg = np.random.choice((0, 1), link_model.channel.nb_tx)
+    link_model.channel.noise_std = 0
+    symbs = link_model.modulate(msg)
+    channel_output = link_model.channel.propagate(symbs)
+    received_msg = link_model.receive(channel_output[0], link_model.channel.channel_gains[0],
+                                      link_model.constellation, link_model.channel.noise_std ** 2)
+    receive_size = len(received_msg)
+
     # Computations
     for id_SNR in range(len(SNRs)):
         link_model.channel.set_SNR_dB(SNRs[id_SNR], code_rate, link_model.Es)
         bit_send = 0
         bit_err = 0
-        receive_size = link_model.channel.nb_tx * link_model.num_bits_symbol
         while bit_send < send_max and bit_err < err_min:
             # Propagate some bits
             msg = np.random.choice((0, 1), send_chunk)
@@ -78,12 +86,12 @@ def link_performance(link_model, SNRs, send_max, err_min, send_chunk=None, code_
                 nb_symb_vector = len(channel_output)
                 received_msg = np.empty_like(msg, int)
                 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],
-                                            link_model.constellation, link_model.channel.noise_std**2)
+                    received_msg[receive_size * i:receive_size * (i + 1)] = \
+                        link_model.receive(channel_output[i], link_model.channel.channel_gains[i],
+                                           link_model.constellation, link_model.channel.noise_std ** 2)
             else:
                 received_msg = link_model.receive(channel_output, link_model.channel.channel_gains,
-                                                  link_model.constellation, link_model.channel.noise_std**2)
+                                                  link_model.constellation, link_model.channel.noise_std ** 2)
             # Count errors
             bit_err += (msg != received_msg).sum()
             bit_send += send_chunk
@@ -144,6 +152,7 @@ class LinkModel:
              Average energy per symbols.
              *Default* Es=1.
         """
+
     def __init__(self, modulate, channel, receive, num_bits_symbol, constellation, Es=1):
         self.modulate = modulate
         self.channel = channel

commpy/modulation.py

@@ -249,7 +249,7 @@ def mimo_ml(y, h, constellation):
     return x_r
 
 
-def kbest(y, h, constellation, K, noise_var=0, output_type='hard', decode=None):
+def kbest(y, h, constellation, K, noise_var=0, output_type='hard', demode=None):
     """ MIMO K-best Schnorr-Euchner Detection.
 
     Reference: Zhan Guo and P. Nilsson, 'Algorithm and implementation of the K-best sphere decoding for MIMO detection',
@@ -278,7 +278,7 @@ def kbest(y, h, constellation, K, noise_var=0, output_type='hard', decode=None):
         'soft': soft output i.e. output is a vector of Log-Likelihood Ratios.
         *Default* value is 'hard'
 
-    decode : function with prototype binary_word = decode(point)
+    demode : function with prototype binary_word = demode(point)
         Function that provide the binary word corresponding to a symbol vector.
 
     Returns
@@ -341,7 +341,7 @@ def kbest(y, h, constellation, K, noise_var=0, output_type='hard', decode=None):
     if output_type == 'hard':
         return X[:, 0]
     elif output_type == 'soft':
-        return max_log_approx(y, h, noise_var, X, decode)
+        return max_log_approx(y, h, noise_var, X, demode)
     else:
         raise ValueError('output_type must be "hard" or "soft"')
 
@@ -372,8 +372,8 @@ def bit_lvl_repr(H, w):
         raise ValueError('Beta must be even.')
 
 
-def max_log_approx(y, h, noise_var, pts_list, decode):
-    """ Max-log approximation
+def max_log_approx(y, h, noise_var, pts_list, demode):
+    """ Max-log demode
 
     parameters
     ----------
@@ -390,7 +390,7 @@ def max_log_approx(y, h, noise_var, pts_list, decode):
         Set of points to compute max log approximation (points are column-wise).
         (shape: num_receive_antennas x num_points)
 
-    decode : function with prototype binary_word = decode(point)
+    demode : function with prototype binary_word = demode(point)
         Function that provide the binary word corresponding to a symbol vector.
 
     return
@@ -400,7 +400,7 @@ def max_log_approx(y, h, noise_var, pts_list, decode):
     """
     # Decode all pts
     nb_pts = pts_list.shape[1]
-    bits = decode(pts_list.reshape(-1, order='F')).reshape(nb_pts, -1)  # Set of binary words (one word by row)
+    bits = demode(pts_list.reshape(-1, order='F')).reshape(nb_pts, -1)  # Set of binary words (one word by row)
 
     # Prepare LLR computation
     nb_bits = bits.shape[1]

commpy/tests/test_modulation.py

@@ -23,7 +23,7 @@ def test_bit_lvl_repr():
 
     SNR = arange(10, 16, 5)
 
-    def receiver_with_blr(y, H, cons):
+    def receiver_with_blr(y, H, cons, noise_var):
         # Create w
         beta = int(log2(len(cons)))
         reel = [pow(2, i) for i in range(beta // 2 - 1, -1, -1)]
@@ -36,7 +36,7 @@ def receiver_with_blr(y, H, cons):
         mes[mes == -1] = 0
         return mes
 
-    def receiver_without_blr(y, H, cons):
+    def receiver_without_blr(y, H, cons, noise_var):
         return qam.demodulate(mimo_ml(y, H, cons), 'hard')
 
     my_model_without_blr = \