9 files changed, 227 insertions, 61 deletions
diff --git a/gr-fec/python/fec/polar/CMakeLists.txt b/gr-fec/python/fec/polar/CMakeLists.txt
index 1362ce18bb..a863995aff 100644
--- a/gr-fec/python/fec/polar/CMakeLists.txt
+++ b/gr-fec/python/fec/polar/CMakeLists.txt
@@ -24,7 +24,7 @@ GR_PYTHON_INSTALL(
     FILES
     __init__.py
     channel_construction.py
-    channel_construction_bsc.py
+    channel_construction_awgn.py
     channel_construction_bec.py
     helper_functions.py
     DESTINATION ${GR_PYTHON_DIR}/gnuradio/fec/polar
diff --git a/gr-fec/python/fec/polar/README.md b/gr-fec/python/fec/polar/README.md
index d425e8650d..42744915ab 100644
--- a/gr-fec/python/fec/polar/README.md
+++ b/gr-fec/python/fec/polar/README.md
@@ -1,9 +1,33 @@
 POLAR Code Python test functions module
 ===========
 
-This directory contains all the necessary files for POLAR code testcode.
-It serves as a reference for C++ implementations.
+This directory contains all the necessary files for POLAR code testcode and channel construction.
 
+Test code
+====
+The test code serves as a reference for the C++ implementations.
+'common', 'encoder' and 'decoder' are part of this test code.
+'testbed' exists to play with this test code.
+
+
+Channel Construction
+=====
 'polar_channel_construction' exposes functionality to calculate polar channels for different sizes.
 It may be used to calculate Bhattacharyya parameters once and store them in a file in '~/.gnuradio/polar'.
-Frozen bit positions are recalculated on every run.
-\ No newline at end of file
+Frozen bit positions are recalculated on every run.
+
+'polar_channel_construction' provides a command-line interface for all the channel construction code.
+These features are also accessible via the Polar Configurator block in GRC.
+
+BEC
+====
+BEC channel construction can be calculated explicitly because the BEC represents a special case which simplifies the task.
+All functionality regarding this channel is located in 'channel_construction_bec'.
+
+AWGN
+====
+In general channel construction for polar codes is a very time consuming task.
+Tal and Vardy proposed a quantization algorithm for channel construction which makes it feasible.
+The corresponding implementation is located in 'channel_construction_awgn'.
+It should be noted that this more of a baseline implementation which lacks all the advanced features the original implementation provides.
+However, simulations show that BEC channel construction with a design-SNR of 0.0dB yields similar performance and should be preferred here.
+\ No newline at end of file
diff --git a/gr-fec/python/fec/polar/__init__.py b/gr-fec/python/fec/polar/__init__.py
index 0b9c264a77..ce1b1459fb 100644
--- a/gr-fec/python/fec/polar/__init__.py
+++ b/gr-fec/python/fec/polar/__init__.py
@@ -26,20 +26,20 @@ from channel_construction_bec import bhattacharyya_bounds
 from helper_functions import is_power_of_two
 
 
-CHANNEL_TYPE_BSC = 'BSC'
+CHANNEL_TYPE_AWGN = 'AWGN'
 CHANNEL_TYPE_BEC = 'BEC'
 
 def get_z_params(is_prototype, channel, block_size, design_snr, mu):
     print('POLAR code channel construction called with parameters channel={0}, blocksize={1}, design SNR={2}, mu={3}'.format(channel, block_size, design_snr, mu))
-    if not (channel == 'BSC' or channel == 'BEC'):
-        raise ValueError("channel is {0}, but only BEC and BSC are supported!".format(channel))
+    if not (channel == 'AWGN' or channel == 'BEC'):
+        raise ValueError("channel is {0}, but only BEC and AWGN are supported!".format(channel))
     if not is_power_of_two(block_size):
         raise ValueError("block size={0} is not a power of 2!".format(block_size))
     if design_snr < -1.5917:
         raise ValueError("design SNR={0} < -1.5917. MUST be greater!".format(design_snr))
     if not mu > 0:
         raise ValueError("mu={0} < 1. MUST be > 1!".format(mu))
-    if not is_prototype and channel == 'BSC':
+    if not is_prototype and channel == 'AWGN':
         z_params = cc.load_z_parameters(block_size, design_snr, mu)
         print('Read Z-parameter file: {0}'.format(cc.default_dir() + cc.generate_filename(block_size, design_snr, mu)))
         return z_params
diff --git a/gr-fec/python/fec/polar/channel_construction.py b/gr-fec/python/fec/polar/channel_construction.py
index 477e1df605..a981007b45 100644
--- a/gr-fec/python/fec/polar/channel_construction.py
+++ b/gr-fec/python/fec/polar/channel_construction.py
@@ -27,7 +27,7 @@ foundational paper for polar codes.
 from channel_construction_bec import calculate_bec_channel_capacities
 from channel_construction_bec import design_snr_to_bec_eta
 from channel_construction_bec import bhattacharyya_bounds
-from channel_construction_bsc import tal_vardy_tpm_algorithm
+from channel_construction_awgn import tal_vardy_tpm_algorithm
 from helper_functions import *
 
 
@@ -91,11 +91,11 @@ def default_dir():
     return path
 
 
-def save_z_parameters(z_params, block_size, design_snr, mu):
+def save_z_parameters(z_params, block_size, design_snr, mu, alt_construction_method='Tal-Vardy algorithm'):
     path = default_dir()
     filename = generate_filename(block_size, design_snr, mu)
     header = Z_PARAM_FIRST_HEADER_LINE + "\n"
-    header += "Channel construction method: Tal-Vardy algorithm\n"
+    header += "Channel construction method: " + alt_construction_method + "\n"
     header += "Parameters:\n"
     header += "block_size=" + str(block_size) + "\n"
     header += "design_snr=" + str(design_snr) + "\n"
diff --git a/gr-fec/python/fec/polar/channel_construction_bsc.py b/gr-fec/python/fec/polar/channel_construction_awgn.py
index 77057a7c1c..7d820b2883 100755
--- a/gr-fec/python/fec/polar/channel_construction_bsc.py
+++ b/gr-fec/python/fec/polar/channel_construction_awgn.py
@@ -34,24 +34,6 @@ from helper_functions import *
 from channel_construction_bec import bhattacharyya_bounds
 
 
-def bsc_channel(p):
-    '''
-    binary symmetric channel (BSC)
-    output alphabet Y = {0, 1} and
-    W(0|0) = W(1|1) and W(1|0) = W(0|1)
-
-    this function returns a prob's vector for a BSC
-    p denotes an erroneous transistion
-    '''
-    if not (p >= 0.0 and p <= 1.0):
-        print "given p is out of range!"
-        return np.array([], dtype=float)
-
-    # 0 -> 0, 0 -> 1, 1 -> 0, 1 -> 1
-    W = np.array([[1 - p, p], [p, 1 - p]], dtype=float)
-    return W
-
-
 def solver_equation(val, s):
     cw_lambda = codeword_lambda_callable(s)
     ic_lambda = instantanious_capacity_callable()
@@ -108,6 +90,10 @@ def discretize_awgn(mu, design_snr):
     for j in range(mu):
         tpm[0][j] = q_function(factor + a[j]) - q_function(factor + a[j + 1])
         tpm[1][j] = q_function(-1. * factor + a[j]) - q_function(-1. * factor + a[j + 1])
+
+    tpm = tpm[::-1]
+    tpm[0] = tpm[0][::-1]
+    tpm[1] = tpm[1][::-1]
     return tpm
 
 
@@ -152,12 +138,13 @@ def upper_bound_z_params(z, block_size, design_snr):
 
 
 def tal_vardy_tpm_algorithm(block_size, design_snr, mu):
+    mu = mu // 2  # make sure algorithm uses only as many bins as specified.
     block_power = power_of_2_int(block_size)
     channels = np.zeros((block_size, 2, mu))
     channels[0] = discretize_awgn(mu, design_snr) * 2
 
     print('Constructing polar code with Tal-Vardy algorithm')
-    print('(block_size = {0}, design SNR = {1}, mu = {2}'.format(block_size, design_snr, mu))
+    print('(block_size = {0}, design SNR = {1}, mu = {2}'.format(block_size, design_snr, 2 * mu))
     show_progress_bar(0, block_size)
     for j in range(0, block_power):
         u = 2 ** j
@@ -171,7 +158,6 @@ def tal_vardy_tpm_algorithm(block_size, design_snr, mu):
 
     z = np.zeros(block_size)
     for i in range(block_size):
-        # z[i] = np.sum(channels[i][1])
         z[i] = bhattacharyya_parameter(channels[i])
 
     z = z[bit_reverse_vector(np.arange(block_size), block_power)]
@@ -263,7 +249,7 @@ def normalize_q(q, tpm):
 
 
 def main():
-    print 'channel construction BSC main'
+    print 'channel construction AWGN main'
     n = 8
     m = 2 ** n
     design_snr = 0.0
diff --git a/gr-fec/python/fec/polar/channel_construction_bec.py b/gr-fec/python/fec/polar/channel_construction_bec.py
index f8f960dfe7..4b35602d96 100644
--- a/gr-fec/python/fec/polar/channel_construction_bec.py
+++ b/gr-fec/python/fec/polar/channel_construction_bec.py
@@ -51,22 +51,45 @@ def calc_one_recursion(iw0):
     return iw1
 
 
+def calculate_bec_channel_capacities_loop(initial_channel, block_power):
+        # compare [0, Arikan] eq. 6
+    iw = np.array([initial_channel, ], dtype=float)
+    for i in range(block_power):
+        iw = calc_one_recursion(iw)
+    return iw
+
+
+def calc_vector_capacities_one_recursion(iw0):
+    degraded = odd_rec(iw0)
+    upgraded = even_rec(iw0)
+    iw1 = np.empty(2 * len(iw0), dtype=degraded.dtype)
+    iw1[0::2] = degraded
+    iw1[1::2] = upgraded
+    return iw1
+
+
+def calculate_bec_channel_capacities_vector(initial_channel, block_power):
+    # compare [0, Arikan] eq. 6
+    # this version is ~ 180 times faster than the loop version with 2**22 synthetic channels
+    iw = np.array([initial_channel, ], dtype=float)
+    for i in range(block_power):
+        iw = calc_vector_capacities_one_recursion(iw)
+    return iw
+
+
 def calculate_bec_channel_capacities(eta, block_size):
     # compare [0, Arikan] eq. 6
     iw = 1 - eta  # holds for BEC as stated in paper
-    iw = np.array([iw, ], dtype=float)
     lw = hf.power_of_2_int(block_size)
-    for i in range(lw):
-        iw = calc_one_recursion(iw)
-    return iw
+    return calculate_bec_channel_capacities_vector(iw, lw)
 
 
 def calculate_z_parameters_one_recursion(z_params):
-    z_next = np.zeros(2 * z_params.size)
-    for i in range(z_params.size):
-        z_sq = z_params[i] ** 2
-        z_next[2 * i] = 2 * z_params[i] - z_sq
-        z_next[2 * i + 1] = z_sq
+    z_next = np.empty(2 * z_params.size, dtype=z_params.dtype)
+    z_sq = z_params ** 2
+    z_low = 2 * z_params - z_sq
+    z_next[0::2] = z_low
+    z_next[1::2] = z_sq
     return z_next
 
 
@@ -80,6 +103,7 @@ def calculate_bec_channel_z_parameters(eta, block_size):
 
 
 def design_snr_to_bec_eta(design_snr):
+    # minimum design snr = -1.5917 corresponds to BER = 0.5
     s = 10. ** (design_snr / 10.)
     return np.exp(-s)
 
@@ -95,21 +119,116 @@ def bhattacharyya_bounds(design_snr, block_size):
     For BEC that translates to capacity(i) = 1 - bhattacharyya(i)
     :return Z-parameters in natural bit-order. Choose according to desired rate.
     '''
-    # minimum design snr = -1.5917 corresponds to BER = 0.5
-    s = 10 ** (design_snr / 10)  # 'initial z parameter'.
-    eta = np.exp(-s)
+    eta = design_snr_to_bec_eta(design_snr)
     return calculate_bec_channel_z_parameters(eta, block_size)
 
 
+def plot_channel_capacities(capacity, save_file=None):
+    block_size = len(capacity)
+    try:
+        import matplotlib.pyplot as plt
+        # FUN with matplotlib LaTeX fonts! http://matplotlib.org/users/usetex.html
+        plt.rc('text', usetex=True)
+        plt.rc('font', family='serif')
+        plt.rc('figure', autolayout=True)
+        plt.plot(capacity)
+        plt.xlim([0, block_size])
+        plt.ylim([-0.01, 1.01])
+        plt.xlabel('synthetic channel number')
+        plt.ylabel('channel capacity')
+        # plt.title('BEC channel construction')
+        plt.grid()
+        plt.gcf().set_size_inches(plt.gcf().get_size_inches() * .5)
+        if save_file:
+            plt.savefig(save_file)
+        plt.show()
+    except ImportError:
+        pass  # only plot in case matplotlib is installed
+
+
+def plot_average_channel_distance(save_file=None):
+    eta = 0.5 #  design_snr_to_bec_eta(-1.5917)
+    powers = np.arange(4, 26)
+
+    try:
+        import matplotlib.pyplot as plt
+        import matplotlib
+        # FUN with matplotlib LaTeX fonts! http://matplotlib.org/users/usetex.html
+        plt.rc('text', usetex=True)
+        plt.rc('font', family='serif')
+        plt.rc('figure', autolayout=True)
+
+        dist = []
+        medians = []
+        initial_channel = 1 - eta
+        for p in powers:
+            bs = int(2 ** p)
+            capacities = calculate_bec_channel_capacities(eta, bs)
+            avg_capacity = np.repeat(initial_channel, len(capacities))
+            averages = np.abs(capacities - avg_capacity)
+            avg_distance = np.sum(averages) / float(len(capacities))
+            dist.append(avg_distance)
+            variance = np.std(averages)
+            medians.append(variance)
+
+        plt.errorbar(powers, dist, yerr=medians)
+        plt.grid()
+        plt.xlabel(r'block size $N$')
+        plt.ylabel(r'$\frac{1}{N} \sum_i |I(W_N^{(i)}) - 0.5|$')
+
+        axes = plt.axes()
+        tick_values = np.array(axes.get_xticks().tolist())
+        tick_labels = np.array(tick_values, dtype=int)
+        tick_labels = ['$2^{' + str(i) + '}$' for i in tick_labels]
+        plt.xticks(tick_values, tick_labels)
+        plt.xlim((powers[0], powers[-1]))
+        plt.ylim((0.2, 0.5001))
+        plt.gcf().set_size_inches(plt.gcf().get_size_inches() * .5)
+        if save_file:
+            plt.savefig(save_file)
+        plt.show()
+    except ImportError:
+        pass
+
+
+def plot_capacity_histogram(design_snr, save_file=None):
+    eta = design_snr_to_bec_eta(design_snr)
+    # capacities = calculate_bec_channel_capacities(eta, block_size)
+    try:
+        import matplotlib.pyplot as plt
+        # FUN with matplotlib LaTeX fonts! http://matplotlib.org/users/usetex.html
+        plt.rc('text', usetex=True)
+        plt.rc('font', family='serif')
+        plt.rc('figure', autolayout=True)
+
+        block_sizes = [32, 128, 512]
+        for b in block_sizes:
+            capacities = calculate_bec_channel_capacities(eta, b)
+            w = 1. / float(len(capacities))
+            weights = [w, ] * b
+            plt.hist(capacities, bins=b, weights=weights, range=(0.95, 1.0))
+        plt.grid()
+        plt.xlabel('synthetic channel capacity')
+        plt.ylabel('normalized item count')
+        print(plt.gcf().get_size_inches())
+        plt.gcf().set_size_inches(plt.gcf().get_size_inches() * .5)
+        if save_file:
+            plt.savefig(save_file)
+        plt.show()
+    except ImportError:
+        pass
+
+
 def main():
     print 'channel construction main'
-    n = 10
-    block_size = 2 ** n
-    design_snr = 1.0
+    n = 11
+    block_size = int(2 ** n)
+    design_snr = -1.59
     eta = design_snr_to_bec_eta(design_snr)
-    print(calculate_bec_channel_z_parameters(eta, block_size))
-    print(calculate_bec_channel_capacities(eta, block_size))
-
+    # print(calculate_bec_channel_z_parameters(eta, block_size))
+    # capacity = calculate_bec_channel_capacities(eta, block_size)
+    # plot_average_channel_distance()
+    calculate_bec_channel_z_parameters(eta, block_size)
 
 if __name__ == '__main__':
     main()
diff --git a/gr-fec/python/fec/polar/decoder.py b/gr-fec/python/fec/polar/decoder.py
index 458858e16d..8748d284f7 100644
--- a/gr-fec/python/fec/polar/decoder.py
+++ b/gr-fec/python/fec/polar/decoder.py
@@ -31,11 +31,11 @@ class PolarDecoder(PolarCommon):
         PolarCommon.__init__(self, n, k, frozen_bit_position, frozenbits)
 
         self.error_probability = 0.1  # this is kind of a dummy value. usually chosen individually.
-        self.bsc_lr = ((1 - self.error_probability) / self.error_probability, self.error_probability / (1 - self.error_probability))
-        self.bsc_llrs = np.log(self.bsc_lr)
+        self.lrs = ((1 - self.error_probability) / self.error_probability, self.error_probability / (1 - self.error_probability))
+        self.llrs = np.log(self.lrs)
 
     def _llr_bit(self, bit):
-        return self.bsc_llrs[bit]
+        return self.llrs[bit]
 
     def _llr_odd(self, la, lb):
         # this functions uses the min-sum approximation
@@ -63,7 +63,7 @@ class PolarDecoder(PolarCommon):
         return ui
 
     def _lr_bit(self, bit):
-        return self.bsc_lr[bit]
+        return self.lrs[bit]
 
     def _lr_odd(self, la, lb):
         # la is upper branch and lb is lower branch
diff --git a/gr-fec/python/fec/polar/helper_functions.py b/gr-fec/python/fec/polar/helper_functions.py
index ca66bf4a50..85140c856f 100644
--- a/gr-fec/python/fec/polar/helper_functions.py
+++ b/gr-fec/python/fec/polar/helper_functions.py
@@ -23,6 +23,24 @@ import time, sys
 import copy
 
 
+def bsc_channel(p):
+    '''
+    binary symmetric channel (BSC)
+    output alphabet Y = {0, 1} and
+    W(0|0) = W(1|1) and W(1|0) = W(0|1)
+
+    this function returns a prob's vector for a BSC
+    p denotes an erroneous transistion
+    '''
+    if not (p >= 0.0 and p <= 1.0):
+        print "given p is out of range!"
+        return np.array([], dtype=float)
+
+    # 0 -> 0, 0 -> 1, 1 -> 0, 1 -> 1
+    W = np.array([[1 - p, p], [p, 1 - p]], dtype=float)
+    return W
+
+
 def power_of_2_int(num):
     return int(np.log2(num))
 
@@ -139,13 +157,21 @@ def mutual_information(w):
 
 
 def bhattacharyya_parameter(w):
-    '''bhattacharyya parameter is a measure of similarity between two prob. distributions'''
-    # sum over all y e Y for sqrt( W(y|0) * W(y|1) )
+    '''
+    bhattacharyya parameter is a measure of similarity between two prob. distributions
+    THEORY: sum over all y e Y for sqrt( W(y|0) * W(y|1) )
+    Implementation:
+    Numpy vector of dimension (2, mu//2)
+    holds probabilities P(x|0), first vector for even, second for odd.
+    '''
     dim = np.shape(w)
-    ydim = dim[0]
-    z = 0.0
-    for y in range(ydim):
-        z += np.sqrt(w[0, y] * w[1, y])
+    if len(dim) != 2:
+        raise ValueError
+
+    if dim[0] > dim[1]:
+        raise ValueError
+
+    z = np.sum(np.sqrt(w[0] * w[1]))
     # need all
     return z
 
@@ -164,6 +190,17 @@ def main():
     print(pos)
     print(rev_pos)
 
+    f = np.linspace(.01, .29, 10)
+    e = np.linspace(.03, .31, 10)
+
+    b = np.array([e, f])
+    zp = bhattacharyya_parameter(b)
+    print(zp)
+
+    a = np.sum(np.sqrt(e * f))
+    print(a)
+
+
 
 if __name__ == '__main__':
     main()
diff --git a/gr-fec/python/fec/polar/polar_channel_construction b/gr-fec/python/fec/polar/polar_channel_construction
index 0aca9c7ddc..448253659e 100644
--- a/gr-fec/python/fec/polar/polar_channel_construction
+++ b/gr-fec/python/fec/polar/polar_channel_construction
@@ -34,7 +34,7 @@ def setup_parser():
     parser.add_option("-h", "--help", action="help", help="Displays this help message.")
 
     parser.add_option("-c", "--channel", action="store", type="string", dest="channel",
-                      help="specify channel, currently BEC or BSC (default='BEC')", default='BEC')
+                      help="specify channel, currently BEC or AWGN (default='BEC')", default='BEC')
     parser.add_option("-b", "--blocksize", action="store", type="int", dest="block_size",
                       help="specify block size of polar code (default=16)", default=16)
     parser.add_option("-s", "--desgin-snr", action="store", type="float", dest="design_snr",