digital: Python functions to support soft decision making and look-up table generation.

author: Tom Rondeau <tom@trondeau.com> 2013-09-04 15:44:20 -0400
committer: Tom Rondeau <tom@trondeau.com> 2013-09-04 15:44:20 -0400
commit: 02fd90029a1ec465f7c8da4a0d97dcdb8b3f438f (patch)
tree: 9bc144a9cb5b0e877108302d13181e50bd885528
parent: edf2ac24c9c2f723eb80918de7addf6d299958c2 (diff)
10 files changed, 1580 insertions, 6 deletions
diff --git a/gr-digital/include/gnuradio/digital/constellation.h b/gr-digital/include/gnuradio/digital/constellation.h
index 6ee274dd04..b2464aa856 100644
--- a/gr-digital/include/gnuradio/digital/constellation.h
+++ b/gr-digital/include/gnuradio/digital/constellation.h
@@ -120,6 +120,76 @@ namespace gr {
 	return shared_from_this();
       }  
 
+      /*! \brief Generates the soft decision LUT based on
+       *         constellation and symbol map.
+       *
+       * \details Generates the soft decision LUT based on
+       * constellation and symbol map.  It can be given a estimate of
+       * the noise power in the channel as \p npwr.
+       *
+       * \param precision Number of bits of precision on each axis.
+       * \param npwr Estimate of the noise power (if known).
+       *
+       * This is expensive to compute.
+       */
+      void gen_soft_dec_lut(int precision, float npwr=1.0);
+
+      /*! \brief Calculate soft decisions for the given \p sample.
+       *
+       * \details Calculate the soft decisions from the given \p sample
+       * at the given noise power \p npwr.
+       *
+       * This is a very costly algorithm (especially for higher order
+       * modulations) and should be used sparingly. It uses the
+       * gen_soft_dec_lut function to generate the LUT, which
+       * should be done once or if a large change in the noise floor
+       * is detected.
+       *
+       * Instead of using this function, generate the LUT using the
+       * gen_soft_dec_lut after creating the constellation object
+       * and then use the soft_decision_maker function to return the
+       * answer from the LUT.
+       *
+       * \param sample The complex sample to get the soft decisions.
+       * \param npwr Estimate of the noise power (if known).
+       */
+      virtual std::vector<float> calc_soft_dec(gr_complex sample, float npwr=1.0);
+
+      /*! \brief Define a soft decision look-up table.
+       *
+       * \details Define a soft decision look-up table (LUT). Because
+       * soft decisions can be calculated in various ways with various
+       * levels of accuracy and complexity, this function allows
+       * users to create a LUT in their own way.
+       *
+       * Setting the LUT here means that has_soft_dec_lut will return
+       * true. Decision vectors returned by soft_decision_maker will
+       * be calculated using this LUT.
+       *
+       * \param soft_dec_lut The soft decision LUT as a vector of
+       *        tuples (vectors in C++) of soft decisions. Each
+       *        element of the LUT is a vector of k-bit floats (where
+       *        there are k bits/sample in the constellation).
+       * \param precision The number of bits of precision used when
+       *        generating the LUT.
+       */
+      void set_soft_dec_lut(const std::vector< std::vector<float> > &soft_dec_lut,
+                            int precision);
+
+      //! Returns True if the soft decision LUT has been defined, False otherwise.
+      bool has_soft_dec_lut();
+      
+      /*! \brief Returns the soft decisions for the given \p sample.
+       *
+       * \details Returns the soft decisions for the given \p
+       * sample. If a LUT is defined for the object, the decisions
+       * will be calculated from there. Otherwise, this function will
+       * call calc_soft_dec directly to calculate the soft decisions.
+       *
+       * \param sample The complex sample to get the soft decisions.
+       */
+      std::vector<float> soft_decision_maker(gr_complex sample);
+
     protected:
       std::vector<gr_complex> d_constellation;
       std::vector<int> d_pre_diff_code; 
@@ -130,10 +200,17 @@ namespace gr {
       //! The factor by which the user given constellation points were
       //! scaled by to achieve an average amplitude of 1.
       float d_scalefactor;
+      float d_re_min, d_re_max, d_im_min, d_im_max;
+
+      std::vector< std::vector<float> > d_soft_dec_lut;
+      int d_lut_precision;
+      float d_lut_scale;
 
       float get_distance(unsigned int index, const gr_complex *sample);
       unsigned int get_closest_point(const gr_complex *sample);
       void calc_arity();
+
+      void max_min_axes();
     };
 
     /************************************************************/
@@ -418,6 +495,10 @@ namespace gr {
     /*! 
      * \brief Digital constellation for QPSK
      * \ingroup digital
+     *
+     * 01 | 11
+     * -------
+     * 00 | 10
      */
     class DIGITAL_API constellation_qpsk : public constellation
     {
@@ -446,6 +527,10 @@ namespace gr {
     /*!
      * \brief Digital constellation for DQPSK
      * \ingroup digital
+     *
+     * 01 | 00
+     * -------
+     * 11 | 10
      */
     class DIGITAL_API constellation_dqpsk : public constellation
     {
@@ -474,6 +559,12 @@ namespace gr {
     /*! 
      * \brief Digital constellation for 8PSK
      * \ingroup digital
+     *
+     *     101 | 100
+     * 001     |     000
+     * -----------------
+     * 011     |     010
+     *     111 | 110
      */
     class DIGITAL_API constellation_8psk : public constellation
     {
diff --git a/gr-digital/lib/constellation.cc b/gr-digital/lib/constellation.cc
index 2c70bb0944..ff86531f0f 100644
--- a/gr-digital/lib/constellation.cc
+++ b/gr-digital/lib/constellation.cc
@@ -33,6 +33,7 @@
 #include <stdlib.h>
 #include <float.h>
 #include <stdexcept>
+#include <boost/format.hpp>
 
 namespace gr {
   namespace digital {
@@ -48,7 +49,9 @@ namespace gr {
       d_constellation(constell),
       d_pre_diff_code(pre_diff_code),
       d_rotational_symmetry(rotational_symmetry),
-      d_dimensionality(dimensionality)
+      d_dimensionality(dimensionality),
+      d_lut_precision(0),
+      d_lut_scale(0)
     {
       // Scale constellation points so that average magnitude is 1.
       float summed_mag = 0;
@@ -59,7 +62,7 @@ namespace gr {
       }
       d_scalefactor = constsize/summed_mag;
       for (unsigned int i=0; i<constsize; i++) {
-	d_constellation[i] = d_constellation[i]*d_scalefactor;
+        d_constellation[i] = d_constellation[i]*d_scalefactor;
       }
       if(pre_diff_code.size() == 0)
 	d_apply_pre_diff_code = false;
@@ -233,10 +236,154 @@ namespace gr {
     }
 
 
+    void
+    constellation::gen_soft_dec_lut(int precision, float npwr)
+    {
+      max_min_axes();
+      d_lut_scale = powf(2.0, static_cast<float>(precision));
+      float xstep = (d_re_max - d_re_min) / (d_lut_scale-1);
+      float ystep = (d_im_max - d_im_min) / (d_lut_scale-1);
+      d_soft_dec_lut.clear();
+
+      float y = d_im_min;
+      while(y < d_im_max) {
+        float x = d_re_min;
+        while(x < d_re_max) {
+          gr_complex pt = gr_complex(x, y);
+          d_soft_dec_lut.push_back(calc_soft_dec(pt, npwr));
+          x += xstep;
+        }
+        y += ystep;
+      }
+
+      d_lut_precision = precision;
+    }
+
+    std::vector<float>
+    constellation::calc_soft_dec(gr_complex sample, float npwr)
+    {
+      int M = static_cast<int>(d_constellation.size());
+      int k = static_cast<int>(log10f(static_cast<float>(M))/log10f(2.0));
+      std::vector<float> tmp(2*k, 0);
+      std::vector<float> s(k, 0);
+
+      float scale = d_scalefactor*d_scalefactor;
+      
+      for(int i = 0; i < M; i++) {
+        // Calculate the distance between the sample and the current
+        // constellation point.
+        float dist = powf(std::abs(sample - d_constellation[i]), 2.0f);
+
+        // Calculate the probability factor from the distance and
+        // the scaled noise power.
+        float d = expf(-dist / (2.0*npwr*scale));
+
+        for(int j = 0; j < k; j++) {
+          // Get the bit at the jth index
+          int mask = 1 << j;
+          int bit = (d_pre_diff_code[i] & mask) >> j;
+
+          // If the bit is a 0, add to the probability of a zero
+          if(bit == 0)
+            tmp[2*j+0] += d;
+          // else, add to the probability of a one
+          else
+            tmp[2*j+1] += d;
+        }
+      }
+
+      // Calculate the log-likelihood ratio for all bits based on the
+      // probability of ones (tmp[2*i+1]) over the probability of a zero
+      // (tmp[2*i+0]).
+      for(int i = 0; i < k; i++) {
+        s[k-1-i] = (logf(tmp[2*i+1]) - logf(tmp[2*i+0])) * scale;
+      }
+
+      return s;
+    }
+
+    void
+    constellation::set_soft_dec_lut(const std::vector< std::vector<float> > &soft_dec_lut,
+                                    int precision)
+    {
+      max_min_axes();
+
+      d_soft_dec_lut = soft_dec_lut;
+      d_lut_precision = precision;
+      d_lut_scale = powf(2.0, static_cast<float>(precision));
+    }
+
+    bool
+    constellation::has_soft_dec_lut()
+    {
+      return d_soft_dec_lut.size() > 0;
+    }
+
+    std::vector<float>
+    constellation::soft_decision_maker(gr_complex sample)
+    {
+      if(has_soft_dec_lut()) {
+        float xre = sample.real();
+        float xim = sample.imag();
+
+        float xstep = (d_re_max-d_re_min) / d_lut_scale;
+        float ystep = (d_im_max-d_im_min) / d_lut_scale;
+
+        float xscale = (d_lut_scale / (d_re_max-d_re_min)) - xstep;
+        float yscale = (d_lut_scale / (d_im_max-d_im_min)) - ystep;
+
+        xre = floorf((-d_re_min + std::min(d_re_max, std::max(d_re_min, xre))) * xscale);
+        xim = floorf((-d_im_min + std::min(d_im_max, std::max(d_im_min, xim))) * yscale);
+        int index = static_cast<int>(d_lut_scale*xim + xre);
+
+        int max_index = d_lut_scale*d_lut_scale;
+        if(index > max_index) {
+          return d_soft_dec_lut[max_index-1];
+        }
+
+        if(index < 0)
+          throw std::runtime_error("constellation::soft_decision_maker: input sample out of range.");
+
+        return d_soft_dec_lut[index];
+      }
+      else {
+        return calc_soft_dec(sample);
+      }
+    }
+
+    void
+    constellation::max_min_axes()
+    {
+      // Find min/max of constellation for both real and imag axes.
+      d_re_min = 1e20;
+      d_im_min = 1e20;
+      d_re_max = -1e20;
+      d_im_max = -1e20;
+      for(size_t i = 0; i < d_constellation.size(); i++) {
+        if(d_constellation[i].real() > d_re_max)
+          d_re_max = d_constellation[i].real();
+        if(d_constellation[i].imag() > d_im_max)
+          d_im_max = d_constellation[i].imag();
+
+        if(d_constellation[i].real() < d_re_min)
+          d_re_min = d_constellation[i].real();
+        if(d_constellation[i].imag() < d_im_min)
+          d_im_min = d_constellation[i].imag();
+      }
+      if(d_im_min == 0)
+        d_im_min = d_re_min;
+      if(d_im_max == 0)
+        d_im_max = d_re_max;
+      if(d_re_min == 0)
+        d_re_min = d_im_min;
+      if(d_re_max == 0)
+        d_re_max = d_im_max;
+    }
+
     /********************************************************************/
 
 
-    constellation_calcdist::sptr 
+    constellation_calcdist::sptr
     constellation_calcdist::make(std::vector<gr_complex> constell,
 				 std::vector<int> pre_diff_code,
 				 unsigned int rotational_symmetry,
diff --git a/gr-digital/python/digital/CMakeLists.txt b/gr-digital/python/digital/CMakeLists.txt
index 30066d9aac..e0cb77fb21 100644
--- a/gr-digital/python/digital/CMakeLists.txt
+++ b/gr-digital/python/digital/CMakeLists.txt
@@ -1,4 +1,4 @@
-# Copyright 2011-212 Free Software Foundation, Inc.
+# Copyright 2011-2013 Free Software Foundation, Inc.
 #
 # This file is part of GNU Radio
 #
@@ -26,6 +26,7 @@ GR_PYTHON_INSTALL(
     FILES
     __init__.py
     bpsk.py
+    constellation_map_generator.py
     cpm.py
     crc.py
     generic_mod_demod.py
@@ -43,9 +44,12 @@ GR_PYTHON_INSTALL(
     packet_utils.py
     pkt.py
     psk.py
+    psk_constellations.py
     qam.py
     qamlike.py
+    qam_constellations.py
     qpsk.py
+    soft_dec_lut_gen.py
     DESTINATION ${GR_PYTHON_DIR}/gnuradio/digital
     COMPONENT "digital_python"
 )
diff --git a/gr-digital/python/digital/__init__.py b/gr-digital/python/digital/__init__.py
index 5059e4eec8..79b740644d 100644
--- a/gr-digital/python/digital/__init__.py
+++ b/gr-digital/python/digital/__init__.py
@@ -1,4 +1,4 @@
-# Copyright 2011 Free Software Foundation, Inc.
+# Copyright 2011-2013 Free Software Foundation, Inc.
 # 
 # This file is part of GNU Radio
 # 
@@ -50,6 +50,10 @@ from ofdm_sync_ml import *
 from ofdm_sync_pnac import *
 from ofdm_sync_pn import *
 from ofdm_txrx import ofdm_tx, ofdm_rx
+from soft_dec_lut_gen import *
+from psk_constellations import *
+from qam_constellations import *
+from constellation_map_generator import *
 
 import packet_utils
 import ofdm_packet_utils
diff --git a/gr-digital/python/digital/constellation_map_generator.py b/gr-digital/python/digital/constellation_map_generator.py
new file mode 100644
index 0000000000..bf689676c3
--- /dev/null
+++ b/gr-digital/python/digital/constellation_map_generator.py
@@ -0,0 +1,52 @@
+#!/usr/bin/env python
+#
+# Copyright 2013 Free Software Foundation, Inc.
+# 
+# This file is part of GNU Radio
+# 
+# GNU Radio is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3, or (at your option)
+# any later version.
+# 
+# GNU Radio is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with GNU Radio; see the file COPYING.  If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street,
+# Boston, MA 02110-1301, USA.
+# 
+
+def constellation_map_generator(basis_cpoints, basis_symbols, k, pi):
+    '''
+    Uses the a basis constellation provided (e.g., from
+    psk_constellation.psk_4()) and the the k and permutation index
+    (pi) to generate a new Gray-coded symbol map to the constellation
+    points provided in the basis.
+
+    The basis_cpoints are the constellation points of the basis
+    constellation, and basis_symbols are the symbols that correspond
+    to the constellation points.
+
+    The selection of k and pi will provide an automorphism the
+    hyperoctahedral group of the basis constellation.
+
+    This function returns a tuple of (constellation_points,
+    symbol_map). The constellation_points is a list of the
+    constellation points in complex space and the symbol_map is a list
+    of the log2(M)-bit symbols for the constellation points (i.e.,
+    symbol_map[i] are the bits associated with
+    constellation_points[i]).
+    '''
+    const_points, s = basis()
+    symbols = list()
+    for s_i in s:
+        tmp = 0
+        for i,p in enumerate(pi):
+            bit = (s_i >> i) & 0x1
+            tmp |= bit << p
+        symbols.append(tmp ^ k)
+    return (const_points, symbols)
diff --git a/gr-digital/python/digital/psk_constellations.py b/gr-digital/python/digital/psk_constellations.py
new file mode 100755
index 0000000000..ee62c33a0c
--- /dev/null
+++ b/gr-digital/python/digital/psk_constellations.py
@@ -0,0 +1,308 @@
+#!/usr/bin/env python
+#
+# Copyright 2013 Free Software Foundation, Inc.
+# 
+# This file is part of GNU Radio
+# 
+# GNU Radio is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3, or (at your option)
+# any later version.
+# 
+# GNU Radio is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with GNU Radio; see the file COPYING.  If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street,
+# Boston, MA 02110-1301, USA.
+# 
+
+import numpy
+from constellation_map_generator import *
+
+'''
+Note on the naming scheme. Each constellation is named using a prefix
+for the type of constellation, the order of the constellation, and a
+distinguishing feature, which comes in three modes:
+
+- No extra feature: the basic Gray-coded constellation map; others
+  will be derived from this type.
+- A single number: an indexed number to uniquely identify different
+  constellation maps.
+- 0xN_x0_x1..._xM: A permutation of the base constellation, explained
+  below.
+
+For rectangular constellations (BPSK, QPSK, QAM), we can define a
+hyperspace and look for all symmetries. This is also known as the
+automorphism group of the hypercube, aka the hyperoctahedral
+group. What this means is that we can easily define all possible
+rotations in terms of the first base mapping by creating the mapping:
+
+  f(x) = k XOR pi(x)
+
+The x is the bit string for the symbol we are altering. Then k is a
+bit string of n bits where n is the number of bits per symbol in the
+constellation (e.g., 2 for QPSK or 6 for QAM64). The pi is a
+permutation function specified as pi_0, pi_1..., pi_n-1. This permutes
+the bits from the base constellation symbol to a new code, which is
+then xor'd by k.
+
+The value of k is from 0 to 2^n-1 and pi is a list of all bit
+positions.
+
+The total number of Gray coded modulations is (2^n)*(n!).
+
+We create aliases for all possible naming schemes for the
+constellations. So if a hyperoctahedral group is defined, we also set
+this function equal to a function name using a unique ID number, and
+we always select one rotation as our basic rotation that the other
+rotations are based off of.
+'''
+
+# BPSK Constellation Mappings
+
+def psk_2_0x0():
+    '''
+    0 | 1
+    '''
+    const_points = [-1, 1]
+    symbols = [0, 1]
+    return (const_points, symbols)
+psk_2 = psk_2_0x0  # Basic BPSK rotation
+psk_2_0 = psk_2    # First ID for BPSK rotations
+
+def psk_2_0x1():
+    '''
+    1 | 0
+    '''
+    const_points = [-1, 1]
+    symbols = [1, 0]
+    return (const_points, symbols)
+psk_2_1 = psk_2_0x1
+
+
+############################################################
+# BPSK Soft bit LUT generators
+############################################################
+
+def sd_psk_2_0x0(x, Es=1):
+    '''
+    0 | 1
+    '''
+    x_re = x.real
+    dist = Es*numpy.sqrt(2)
+    return [dist*x_re,]
+sd_psk_2 = sd_psk_2_0x0  # Basic BPSK rotation
+sd_psk_2_0 = sd_psk_2    # First ID for BPSK rotations
+
+def sd_psk_2_0x1(x, Es=1):
+    '''
+    1 | 0
+    '''
+    x_re = [x.real,]
+    dist = Es*numpy.sqrt(2)
+    return -dist*x_re
+sd_psk_2_1 = sd_psk_2_0x1
+
+
+############################################################
+# QPSK Constellation Mappings
+############################################################
+
+def psk_4_0x0_0_1():
+    '''
+    10 | 11
+    -------
+    00 | 01
+    '''
+    const_points = [-1-1j, 1-1j,
+                    -1+1j, 1+1j]
+    symbols = [0, 1, 2, 3]
+    return (const_points, symbols)
+psk_4 = psk_4_0x0_0_1
+psk_4_0 = psk_4
+
+def psk_4_0x1_0_1():
+    '''
+    11 | 10
+    -------
+    01 | 00
+    '''
+    k = 0x1
+    pi = [0, 1]
+    return constellation_map_generator(psk_4, k, pi)
+psk_4_1 = psk_4_0x1_0_1
+
+def psk_4_0x2_0_1():
+    '''
+    00 | 01
+    -------
+    10 | 11
+    '''
+    k = 0x2
+    pi = [0, 1]
+    return constellation_map_generator(psk_4, k, pi)
+psk_4_2 = psk_4_0x2_0_1
+
+def psk_4_0x3_0_1():
+    '''
+    01 | 00
+    -------
+    11 | 10
+    '''
+    k = 0x3
+    pi = [0, 1]
+    return constellation_map_generator(psk_4, k, pi)
+psk_4_3 = psk_4_0x3_0_1
+
+def psk_4_0x0_1_0():
+    '''
+    01 | 11
+    -------
+    00 | 10
+    '''
+    k = 0x0
+    pi = [1, 0]
+    return constellation_map_generator(psk_4, k, pi)
+psk_4_4 = psk_4_0x0_1_0
+
+def psk_4_0x1_1_0():
+    '''
+    00 | 10
+    -------
+    01 | 11
+    '''
+    k = 0x1
+    pi = [1, 0]
+    return constellation_map_generator(psk_4, k, pi)
+psk_4_5 = psk_4_0x1_1_0
+
+def psk_4_0x2_1_0():
+    '''
+    11 | 01
+    -------
+    10 | 00
+    '''
+    k = 0x2
+    pi = [1, 0]
+    return constellation_map_generator(psk_4, k, pi)
+psk_4_6 = psk_4_0x2_1_0
+
+def psk_4_0x3_1_0():
+    '''
+    10 | 00
+    -------
+    11 | 01
+    '''
+    k = 0x3
+    pi = [1, 0]
+    return constellation_map_generator(psk_4, k, pi)
+psk_4_7 = psk_4_0x3_1_0
+
+
+
+############################################################
+# QPSK Constellation Softbit LUT generators
+############################################################
+
+def sd_psk_4_0x0_0_1(x, Es=1):
+    '''
+    10 | 11
+    -------
+    00 | 01
+    '''
+    x_re = x.real
+    x_im = x.imag
+    dist = Es*numpy.sqrt(2)
+    return [dist*x_im, dist*x_re]
+sd_psk_4 = sd_psk_4_0x0_0_1
+sd_psk_4_0 = sd_psk_4
+
+def sd_psk_4_0x1_0_1(x, Es=1):
+    '''
+    11 | 10
+    -------
+    01 | 00
+    '''
+    x_re = x.real
+    x_im = x.imag
+    dist = Es*numpy.sqrt(2)
+    return [dist*x_im, -dist*x_re]
+sd_psk_4_1 = sd_psk_4_0x1_0_1
+
+def sd_psk_4_0x2_0_1(x, Es=1):
+    '''
+    00 | 01
+    -------
+    10 | 11
+    '''
+    x_re = x.real
+    x_im = x.imag
+    dist = Es*numpy.sqrt(2)
+    return [-dist*x_im, dist*x_re]
+sd_psk_4_2 = sd_psk_4_0x2_0_1
+
+def sd_psk_4_0x3_0_1(x, Es=1):
+    '''
+    01 | 00
+    -------
+    11 | 10
+    '''
+    x_re = x.real
+    x_im = x.imag
+    dist = Es*numpy.sqrt(2)
+    return [-dist*x_im, -dist*x_re]
+sd_psk_4_3 = sd_psk_4_0x3_0_1
+
+def sd_psk_4_0x0_1_0(x, Es=1):
+    '''
+    01 | 11
+    -------
+    00 | 10
+    '''
+    x_re = x.real
+    x_im = x.imag
+    dist = Es*numpy.sqrt(2)
+    return [dist*x_re, dist*x_im]
+sd_psk_4_4 = sd_psk_4_0x0_1_0
+
+def sd_psk_4_0x1_1_0(x, Es=1):
+    '''
+    00 | 10
+    -------
+    01 | 11
+    '''
+    x_re = x.real
+    x_im = x.imag
+    dist = Es*numpy.sqrt(2)
+    return [dist*x_re, -dist*x_im]
+sd_psk_4_5 = sd_psk_4_0x1_1_0
+
+
+def sd_psk_4_0x2_1_0(x, Es=1):
+    '''
+    11 | 01
+    -------
+    10 | 00
+    '''
+    x_re = x.real
+    x_im = x.imag
+    dist = Es*numpy.sqrt(2)
+    return [-dist*x_re, dist*x_im]
+sd_psk_4_6 = sd_psk_4_0x2_1_0
+
+def sd_psk_4_0x3_1_0(x, Es=1):
+    '''
+    10 | 00
+    -------
+    11 | 01
+    '''
+    x_re = x.real
+    x_im = x.imag
+    dist = Es*numpy.sqrt(2)
+    return [-dist*x_re, -dist*x_im]
+sd_psk_4_7 = sd_psk_4_0x3_1_0
+
diff --git a/gr-digital/python/digital/qa_constellation.py b/gr-digital/python/digital/qa_constellation.py
index 83a875af4b..00248f4f09 100644
--- a/gr-digital/python/digital/qa_constellation.py
+++ b/gr-digital/python/digital/qa_constellation.py
@@ -21,7 +21,7 @@
 # 
 
 import random
-from cmath import exp, pi, log
+from cmath import exp, pi, log, sqrt
 
 from gnuradio import gr, gr_unittest, digital, blocks
 from gnuradio.digital.utils import mod_codes
@@ -185,6 +185,86 @@ class test_constellation(gr_unittest.TestCase):
                 first = constellation.bits_per_symbol()
                 self.assertEqual(self.src_data[first:len(data)], data[first:])
 
+    def test_soft_qpsk_gen(self):
+        prec = 8
+        constel, code = digital.psk_4_0()
+
+        rot_sym = 1
+        side = 2
+        width = 2
+        c = digital.constellation_rect(constel, code, rot_sym,
+                                       side, side, width, width)
+
+        # Get max energy/symbol in constellation
+        constel = c.points()
+        Es = max([abs(constel_i) for constel_i in constel])
+
+        table = digital.soft_dec_table_generator(digital.sd_psk_4_0, prec, Es)
+        c.set_soft_dec_lut(table, prec)
+
+        x = sqrt(2.0)/2.0
+        step = (x.real+x.real) / (2**prec - 1)
+        samples = [ -x-x*1j, -x+x*1j,
+                     x+x*1j,  x-x*1j,
+                   (-x+128*step)+(-x+128*step)*1j,
+                   (-x+64*step) +(-x+64*step)*1j,  (-x+64*step) +(-x+192*step)*1j,
+                   (-x+192*step)+(-x+192*step)*1j, (-x+192*step)+(-x+64*step)*1j,]
+
+        y_python_raw_calc = []
+        y_python_gen_calc = []
+        y_python_table = []
+        y_cpp_raw_calc = []
+        y_cpp_table = []
+        for sample in samples:
+            y_python_raw_calc += digital.calc_soft_dec(sample, constel, code)
+            y_python_gen_calc += digital.sd_psk_4_0(sample, Es)
+            y_python_table += digital.calc_soft_dec_from_table(sample, table, prec, Es)
+
+            y_cpp_raw_calc += c.calc_soft_dec(sample)
+            y_cpp_table += c.soft_decision_maker(sample)
+
+        self.assertFloatTuplesAlmostEqual(y_python_raw_calc, y_python_gen_calc, 4)
+        self.assertFloatTuplesAlmostEqual(y_python_raw_calc, y_python_table, 2)
+        self.assertFloatTuplesAlmostEqual(y_cpp_raw_calc, y_cpp_table, 4)
+        
+    def test_soft_qpsk_calc(self):
+        prec = 8
+        constel, code = digital.psk_4_0()
+
+        rot_sym = 1
+        side = 2
+        width = 2
+        c = digital.constellation_rect(constel, code, rot_sym,
+                                       side, side, width, width)
+
+        # Get max energy/symbol in constellation
+        constel = c.points()
+        Es = max([abs(constel_i) for constel_i in constel])
+
+        table = digital.soft_dec_table(constel, code, prec)
+        c.gen_soft_dec_lut(prec)
+
+        x = sqrt(2.0)/2.0
+        step = (x.real+x.real) / (2**prec - 1)
+        samples = [ -x-x*1j, -x+x*1j,
+                     x+x*1j,  x-x*1j,
+                   (-x+128*step)+(-x+128*step)*1j,
+                   (-x+64*step) +(-x+64*step)*1j,  (-x+64*step) +(-x+192*step)*1j,
+                   (-x+192*step)+(-x+192*step)*1j, (-x+192*step)+(-x+64*step)*1j,]
+
+        y_python_raw_calc = []
+        y_python_table = []
+        y_cpp_raw_calc = []
+        y_cpp_table = []
+        for sample in samples:
+            y_python_raw_calc += digital.calc_soft_dec(sample, constel, code)
+            y_python_table += digital.calc_soft_dec_from_table(sample, table, prec, Es)
+
+            y_cpp_raw_calc += c.calc_soft_dec(sample)
+            y_cpp_table += c.soft_decision_maker(sample)
+
+        self.assertFloatTuplesAlmostEqual(y_python_raw_calc, y_python_table, 4)
+        self.assertFloatTuplesAlmostEqual(y_cpp_raw_calc, y_cpp_table, 4)
 
 class mod_demod(gr.hier_block2):
     def __init__(self, constellation, differential, rotation):
diff --git a/gr-digital/python/digital/qam_constellations.py b/gr-digital/python/digital/qam_constellations.py
new file mode 100755
index 0000000000..4e8ee08a61
--- /dev/null
+++ b/gr-digital/python/digital/qam_constellations.py
@@ -0,0 +1,521 @@
+#!/usr/bin/env python
+#
+# Copyright 2013 Free Software Foundation, Inc.
+# 
+# This file is part of GNU Radio
+# 
+# GNU Radio is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3, or (at your option)
+# any later version.
+# 
+# GNU Radio is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with GNU Radio; see the file COPYING.  If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street,
+# Boston, MA 02110-1301, USA.
+# 
+
+import numpy
+from constellation_map_generator import *
+
+'''
+Note on the naming scheme. Each constellation is named using a prefix
+for the type of constellation, the order of the constellation, and a
+distinguishing feature, which comes in three modes:
+
+- No extra feature: the basic Gray-coded constellation map; others
+  will be derived from this type.
+- A single number: an indexed number to uniquely identify different
+  constellation maps.
+- 0xN_x0_x1..._xM: A permutation of the base constellation, explained
+  below.
+
+For rectangular constellations (BPSK, QPSK, QAM), we can define a
+hyperspace and look for all symmetries. This is also known as the
+automorphism group of the hypercube, aka the hyperoctahedral
+group. What this means is that we can easily define all possible
+rotations in terms of the first base mapping by creating the mapping:
+
+  f(x) = k XOR pi(x)
+
+The x is the bit string for the symbol we are altering. Then k is a
+bit string of n bits where n is the number of bits per symbol in the
+constellation (e.g., 2 for QPSK or 6 for QAM64). The pi is a
+permutation function specified as pi_0, pi_1..., pi_n-1. This permutes
+the bits from the base constellation symbol to a new code, which is
+then xor'd by k.
+
+The value of k is from 0 to 2^n-1 and pi is a list of all bit
+positions.
+
+The permutation are given for b0_b1_b2_... for the total number of
+bits. In the constellation diagrams shown in the comments, the bit
+ordering is shown as [b3b2b1b0]. Bit values returned from the soft bit
+LUTs are in the order [b3, b2, b1, b0].
+
+
+The total number of Gray coded modulations is (2^n)*(n!).
+
+We create aliases for all possible naming schemes for the
+constellations. So if a hyperoctahedral group is defined, we also set
+this function equal to a function name using a unique ID number, and
+we always select one rotation as our basic rotation that the other
+rotations are based off of.
+
+For 16QAM:
+ - n = 4
+ - (2^n)*(n!) = 384
+ - k \in [0x0, 0xF]
+ - pi = 0, 1, 2, 3
+        0, 1, 3, 2
+        0, 2, 1, 3
+        0, 2, 3, 1
+        0, 3, 1, 2
+        0, 3, 2, 1
+        1, 0, 2, 3
+        1, 0, 3, 2
+        1, 2, 0, 3
+        1, 2, 3, 0
+        1, 3, 0, 2
+        1, 3, 2, 0
+        2, 0, 1, 3
+        2, 0, 3, 1
+        2, 1, 0, 3
+        2, 1, 3, 0
+        2, 3, 0, 1
+        2, 3, 1, 0
+        3, 0, 1, 2
+        3, 0, 2, 1
+        3, 1, 0, 2
+        3, 1, 2, 0
+        3, 2, 0, 1
+        3, 2, 1, 0
+'''
+
+def qam_16_0x0_0_1_2_3():
+    '''
+    0010  0110 | 1110  1010
+
+    0011  0111 | 1111  1011
+    -----------------------
+    0001  0101 | 1101  1001
+
+    0000  0100 | 1100  1000
+    '''
+    const_points = [-3-3j, -1-3j, 1-3j, 3-3j,
+                    -3-1j, -1-1j, 1-1j, 3-1j,
+                    -3+1j, -1+1j, 1+1j, 3+1j,
+                    -3+3j, -1+3j, 1+3j, 3+3j]
+    symbols = [0x0, 0x4, 0xC, 0x8,
+               0x1, 0x5, 0xD, 0x9,
+               0x3, 0x7, 0xF, 0xB,
+               0x2, 0x6, 0xE, 0xA]
+    return (const_points, symbols)
+qam_16 = qam_16_0x0_0_1_2_3
+qam_16_0 = qam_16
+
+def qam_16_0x1_0_1_2_3():
+    '''
+    0011  0111 | 1111  1011
+
+    0010  0110 | 1110  1010
+    -----------------------
+    0000  0100 | 1100  1000
+
+    0001  0101 | 1101  1001
+    '''
+    k = 0x1
+    pi = [0, 1, 2, 3]
+    return constellation_map_generator(qam_16, k, pi)
+qam_16_1 = qam_16_0x1_0_1_2_3
+
+def qam_16_0x2_0_1_2_3():
+    '''
+    0000  0100 | 1100  1000
+
+    0001  0101 | 1101  1001
+    -----------------------
+    0011  0111 | 1111  1011
+
+    0010  0110 | 1110  1010
+    '''
+    k = 0x2
+    pi = [0, 1, 2, 3]
+    return constellation_map_generator(qam_16, k, pi)
+qam_16_2 = qam_16_0x2_0_1_2_3
+
+def qam_16_0x3_0_1_2_3():
+    '''
+    0001  0101 | 1101  1001
+
+    0000  0100 | 1100  1000
+    -----------------------
+    0010  0110 | 1110  1010
+
+    0011  0111 | 1111  1011
+    '''
+    k = 0x3
+    pi = [0, 1, 2, 3]
+    return constellation_map_generator(qam_16, k, pi)
+qam_16_3 = qam_16_0x3_0_1_2_3
+
+
+def qam_16_0x0_1_0_2_3():
+    '''
+    0001  0101 | 1101  1001
+
+    0011  0111 | 1111  1011
+    -----------------------
+    0010  0110 | 1110  1010
+
+    0000  0100 | 1100  1000
+    '''
+    k = 0x0
+    pi = [1, 0, 2, 3]
+    return constellation_map_generator(qam_16, k, pi)
+qam_16_4 = qam_16_0x0_1_0_2_3
+
+def qam_16_0x1_1_0_2_3():
+    '''
+    0000  0100 | 1100  1000
+
+    0010  0110 | 1110  1010
+    -----------------------
+    0011  0111 | 1111  1011
+
+    0001  0101 | 1101  1001
+    '''
+    k = 0x1
+    pi = [1, 0, 2, 3]
+    return constellation_map_generator(qam_16, k, pi)
+qam_16_5 = qam_16_0x1_1_0_2_3
+
+def qam_16_0x2_1_0_2_3():
+    '''
+    0011  0111 | 1111  1011
+
+    0001  0101 | 1101  1001
+    -----------------------
+    0000  0100 | 1100  1000
+
+    0010  0110 | 1110  1010
+    '''
+    k = 0x2
+    pi = [1, 0, 2, 3]
+    return constellation_map_generator(qam_16, k, pi)
+qam_16_6 = qam_16_0x2_1_0_2_3
+
+def qam_16_0x3_1_0_2_3():
+    '''
+    0010  0110 | 1110  1010
+
+    0000  0100 | 1100  1000
+    -----------------------
+    0001  0101 | 1101  1001
+
+    0011  0111 | 1111  1011
+    '''
+    k = 0x3
+    pi = [1, 0, 2, 3]
+    return constellation_map_generator(qam_16, k, pi)
+qam_16_7 = qam_16_0x3_1_0_2_3
+
+
+# Soft bit LUT generators
+
+def sd_qam_16_0x0_0_1_2_3(x, Es=1):
+    '''
+    Soft bit LUT generator for constellation:
+
+    0010  0110 | 1110  1010
+
+    0011  0111 | 1111  1011
+    -----------------------
+    0001  0101 | 1101  1001
+
+    0000  0100 | 1100  1000
+    '''
+
+    dist = Es*numpy.sqrt(2)
+    boundary = dist/3.0
+    dist0 = dist/6.0
+#    print "Sample:    ", x
+#    print "Es:        ", Es
+#    print "Distance:  ", dist
+#    print "Boundary:  ", boundary
+#    print "1st Bound: ", dist0
+
+    x_re = x.real
+    x_im = x.imag
+
+    if x_re < -boundary:
+        b3 = boundary*(x_re + dist0)
+    elif x_re < boundary:
+        b3 = x_re
+    else:
+        b3 = boundary*(x_re - dist0)
+
+    if x_im < -boundary:
+        b1 = boundary*(x_im + dist0)
+    elif x_im < boundary:
+        b1 = x_im
+    else:
+        b1 = boundary*(x_im - dist0)
+
+    b2 = -abs(x_re) + boundary
+    b0 = -abs(x_im) + boundary
+
+    return [(Es/2)*b3, (Es/2)*b2, (Es/2)*b1, (Es/2)*b0]
+sd_qam_16 = sd_qam_16_0x0_0_1_2_3
+sd_qam_16_0 = sd_qam_16
+
+def sd_qam_16_0x1_0_1_2_3(x, Es=1):
+    '''
+    Soft bit LUT generator for constellation:
+
+    0011  0111 | 1111  1011
+
+    0010  0110 | 1110  1010
+    -----------------------
+    0000  0100 | 1100  1000
+
+    0001  0101 | 1101  1001
+    '''
+    x_re = 3*x.real
+    x_im = 3*x.imag
+
+    if x_re < -2:
+        b3 = 2*(x_re + 1)
+    elif x_re < 2:
+        b3 = x_re
+    else:
+        b3 = 2*(x_re - 1)
+
+    if x_im < -2:
+        b1 = 2*(x_im + 1)
+    elif x_im < 2:
+        b1 = x_im
+    else:
+        b1 = 2*(x_im - 1)
+
+    b2 = -abs(x_re) + 2
+    b0 = +abs(x_im) - 2
+
+    return [b3, b2, b1, b0]
+sd_qam_16_1 = sd_qam_16_0x1_0_1_2_3
+
+def sd_qam_16_0x2_0_1_2_3(x, Es=1):
+    '''
+    Soft bit LUT generator for constellation:
+
+    0000  0100 | 1100  1000
+
+    0001  0101 | 1101  1001
+    -----------------------
+    0011  0111 | 1111  1011
+
+    0010  0110 | 1110  1010
+    '''
+
+    x_re = 3*x.real
+    x_im = 3*x.imag
+
+    if x_re < -2:
+        b3 = 2*(x_re + 1)
+    elif x_re < 2:
+        b3 = x_re
+    else:
+        b3 = 2*(x_re - 1)
+
+    if x_im < -2:
+        b1 = -2*(x_im + 1)
+    elif x_im < 2:
+        b1 = -x_im
+    else:
+        b1 = -2*(x_im - 1)
+
+    b2 = -abs(x_re) + 2
+    b0 = -abs(x_im) + 2
+
+    return [b3, b2, b1, b0]
+sd_qam_16_2 = sd_qam_16_0x2_0_1_2_3
+
+def sd_qam_16_0x3_0_1_2_3(x, Es=1):
+    '''
+    Soft bit LUT generator for constellation:
+
+    0001  0101 | 1101  1001
+
+    0000  0100 | 1100  1000
+    -----------------------
+    0010  0110 | 1110  1010
+
+    0011  0111 | 1111  1011
+    '''
+    x_re = 3*x.real
+    x_im = 3*x.imag
+
+    if x_re < -2:
+        b3 = 2*(x_re + 1)
+    elif x_re < 2:
+        b3 = x_re
+    else:
+        b3 = 2*(x_re - 1)
+
+    if x_im < -2:
+        b1 = -2*(x_im + 1)
+    elif x_im < 2:
+        b1 = -x_im
+    else:
+        b1 = -2*(x_im - 1)
+
+    b2 = -abs(x_re) + 2
+    b0 = +abs(x_im) - 2
+
+    return [b3, b2, b1, b0]
+sd_qam_16_3 = sd_qam_16_0x3_0_1_2_3
+
+def sd_qam_16_0x0_1_0_2_3(x, Es=1):
+    '''
+    Soft bit LUT generator for constellation:
+
+    0001  0101 | 1101  1001
+
+    0011  0111 | 1111  1011
+    -----------------------
+    0010  0110 | 1110  1010
+
+    0000  0100 | 1100  1000
+    '''
+    x_re = 3*x.real
+    x_im = 3*x.imag
+
+    if x_re < -2:
+        b3 = 2*(x_re + 1)
+    elif x_re < 2:
+        b3 = x_re
+    else:
+        b3 = 2*(x_re - 1)
+
+    if x_im < -2:
+        b0 = 2*(x_im + 1)
+    elif x_im < 2:
+        b0 = x_im
+    else:
+        b0 = 2*(x_im - 1)
+
+    b2 = -abs(x_re) + 2
+    b1 = -abs(x_im) + 2
+
+    return [b3, b2, b1, b0]
+sd_qam_16_4 = sd_qam_16_0x0_1_0_2_3
+
+def sd_qam_16_0x1_1_0_2_3(x, Es=1):
+    '''
+    Soft bit LUT generator for constellation:
+
+    0000  0100 | 1100  1000
+
+    0010  0110 | 1110  1010
+    -----------------------
+    0011  0111 | 1111  1011
+
+    0001  0101 | 1101  1001
+    '''
+    x_re = 3*x.real
+    x_im = 3*x.imag
+
+    if x_re < -2:
+        b3 = 2*(x_re + 1)
+    elif x_re < 2:
+        b3 = x_re
+    else:
+        b3 = 2*(x_re - 1)
+
+    if x_im < -2:
+        b0 = -2*(x_im + 1)
+    elif x_im < 2:
+        b0 = -x_im
+    else:
+        b0 = -2*(x_im - 1)
+
+    b2 = -abs(x_re) + 2
+    b1 = -abs(x_im) + 2
+
+    return [b3, b2, b1, b0]
+sd_qam_16_5 = sd_qam_16_0x1_1_0_2_3
+
+def sd_qam_16_0x2_1_0_2_3(x, Es=1):
+    '''
+    Soft bit LUT generator for constellation:
+
+    0011  0111 | 1111  1011
+
+    0001  0101 | 1101  1001
+    -----------------------
+    0000  0100 | 1100  1000
+
+    0010  0110 | 1110  1010
+    '''
+    x_re = 3*x.real
+    x_im = 3*x.imag
+
+    if x_re < -2:
+        b3 = 2*(x_re + 1)
+    elif x_re < 2:
+        b3 = x_re
+    else:
+        b3 = 2*(x_re - 1)
+
+    if x_im < -2:
+        b0 = 2*(x_im + 1)
+    elif x_im < 2:
+        b0 = x_im
+    else:
+        b0 = 2*(x_im - 1)
+
+    b2 = -abs(x_re) + 2
+    b1 = +abs(x_im) - 2
+
+    return [b3, b2, b1, b0]
+sd_qam_16_6 = sd_qam_16_0x2_1_0_2_3
+
+def sd_qam_16_0x3_1_0_2_3(x, Es=1):
+    '''
+    Soft bit LUT generator for constellation:
+
+    0010  0110 | 1110  1010
+
+    0000  0100 | 1100  1000
+    -----------------------
+    0001  0101 | 1101  1001
+
+    0011  0111 | 1111  1011
+    '''
+    x_re = 3*x.real
+    x_im = 3*x.imag
+
+    if x_re < -2:
+        b3 = 2*(x_re + 1)
+    elif x_re < 2:
+        b3 = x_re
+    else:
+        b3 = 2*(x_re - 1)
+
+    if x_im < -2:
+        b0 = -2*(x_im + 1)
+    elif x_im < 2:
+        b0 = -x_im
+    else:
+        b0 = -2*(x_im - 1)
+
+    b2 = -abs(x_re) + 2
+    b1 = +abs(x_im) - 2
+
+    return [b3, b2, b1, b0]
+sd_qam_16_7 = sd_qam_16_0x3_1_0_2_3
diff --git a/gr-digital/python/digital/soft_dec_lut_gen.py b/gr-digital/python/digital/soft_dec_lut_gen.py
new file mode 100644
index 0000000000..24d8bbdc2e
--- /dev/null
+++ b/gr-digital/python/digital/soft_dec_lut_gen.py
@@ -0,0 +1,232 @@
+#!/usr/bin/env python
+#
+# Copyright 2013 Free Software Foundation, Inc.
+# 
+# This file is part of GNU Radio
+# 
+# GNU Radio is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3, or (at your option)
+# any later version.
+# 
+# GNU Radio is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with GNU Radio; see the file COPYING.  If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street,
+# Boston, MA 02110-1301, USA.
+# 
+
+import numpy
+
+def soft_dec_table_generator(soft_dec_gen, prec, Es=1):
+    '''
+    Builds a LUT that is a list of tuples. The tuple represents the
+    soft decisions for the constellation/bit mapping at any given
+    point in the complex space, (x,y).
+
+    The table is built to a precision specified by the 'prec'
+    argument. There are (2x2)^prec samples in the sample space, so we
+    get the precision of 2^prec samples in both the real and imaginary
+    axes.
+
+    The space is represented where index 0 is the bottom left corner
+    and the maximum index is the upper left. The table index for a
+    surface space with 4 bits of precision looks like the following:
+
+    240 241 242 243 244 245 246 247 | 248 249 250 251 252 253 254 255
+    224 225 226 227 228 229 230 231 | 232 233 234 235 236 237 238 239
+    208 209 210 211 212 213 214 215 | 216 217 218 219 220 221 222 223
+    192 193 194 195 196 197 198 199 | 200 201 202 203 204 205 206 207
+    176 177 178 179 180 181 182 183 | 184 185 186 187 188 189 190 191
+    160 161 162 163 164 165 166 167 | 168 169 170 171 172 173 174 175
+    144 145 146 147 148 149 150 151 | 152 153 154 155 156 157 158 159
+    128 129 130 131 132 133 134 135 | 136 137 138 139 140 141 142 143
+    -----------------------------------------------------------------
+    112 113 114 115 116 117 118 119 | 120 121 122 123 124 125 126 127
+     96  97  98  99 100 101 102 103 | 104 105 106 107 108 109 110 111
+     80  81  82  83  84  85  86  87 |  88  89  90  91  92  93  94  95
+     64  65  66  67  68  69  70  71 |  72  73  74  75  76  77  78  79
+     48  49  50  51  52  53  54  55 |  56  57  58  59  60  61  62  63
+     32  33  34  35  36  37  38  39 |  40  41  42  43  44  45  46  47
+     16  17  18  19  20  21  22  23 |  24  25  26  27  28  29  30  31
+      0   1   2   3   4   5   6   7 |   8   9  10  11  12  13  14  15
+
+    We then calculate coordinates from -1 to 1 with 2^prec points for
+    both the x and y axes. We then sample starting at (-1, -1) and
+    move left to right on the x-axis and then move up a row on the
+    y-axis. For every point in this sampled space, we calculate the
+    soft decisions for the given constellation/mapping. This is done
+    by passing in the function 'soft_dec_gen' as an argument to this
+    function. This takes in the x/y coordinates and outputs the soft
+    decisions. These soft decisions are stored into the list at the
+    index from the above table as a tuple.
+
+    The function 'calc_from_table' takes in a point and reverses this
+    operation. It converts the point from the coordinates (-1,-1) to
+    (1,1) into an index value in the table and returns the tuple of
+    soft decisions at that index.
+
+    Es is the energy per symbol. This is passed to the function to
+    provide the bounds when calling the generator function since they
+    don't know how the constellation was normalized. Using the
+    (maximum) energy per symbol for constellation allows us to provide
+    any scaling of the constellation (normalized to sum to 1,
+    normalized so the outside points sit on +/-1, etc.) but still
+    calculate the soft decisions as we would given the full
+    constellation.
+    '''
+
+    npts = 2.0**prec
+    maxd = Es*numpy.sqrt(2)/2
+    yrng = numpy.linspace(-maxd, maxd, npts)
+    xrng = numpy.linspace(-maxd, maxd, npts)
+
+    table = []
+    for y in yrng:
+        for x in xrng:
+            pt = complex(x, y)
+            decs = soft_dec_gen(pt, Es)
+            table.append(decs)
+    return table
+
+def soft_dec_table(constel, symbols, prec, npwr=1):
+    '''
+    Similar in nature to soft_dec_table_generator above. Instead, this
+    takes in the constellation and symbol points along with the noise
+    power estimate and uses calc_soft_dec (below) to generate the
+    LUT.
+
+    Instead of assuming that the constellation is normalied (e.g., all
+    points are between -1 and 1), this function calculates the min/max
+    of both the real and imaginary axes and uses those when
+    constructing the LUT. So when using this version of the LUT, the
+    samples and the constellations must be working on the same
+    magnitudes.
+
+    Because this uses the calc_soft_dec function, it can be quite
+    a bit more expensive to generate the LUT, though it should be
+    one-time work.
+    '''
+
+    re_min = min(numpy.array(constel).real)
+    im_min = min(numpy.array(constel).imag)
+    re_max = max(numpy.array(constel).real)
+    im_max = max(numpy.array(constel).imag)
+
+    npts = 2.0**prec
+    yrng = numpy.linspace(im_min, im_max, npts)
+    xrng = numpy.linspace(re_min, re_max, npts)
+
+    table = []
+    for y in yrng:
+        for x in xrng:
+            pt = complex(x, y)
+            decs = calc_soft_dec(pt, constel, symbols, npwr)
+            table.append(decs)
+    return table
+
+def calc_soft_dec_from_table(sample, table, prec, Es=1):
+    '''
+    Takes in a complex sample and converts it from the coordinates
+    (-1,-1) to (1,1) into an index value. The index value points to a
+    location in the provided LUT 'table' and returns the soft
+    decisions tuple at that index.
+
+    sample: the complex sample to calculate the soft decisions
+    from.
+
+    table: the LUT.
+
+    prec: the precision used when generating the LUT.
+
+    Es: the energy per symbol. This is passed to the function to
+    provide the bounds when calling the generator function since they
+    don't know how the constellation was normalized. Using the
+    (maximum) energy per symbol for constellation allows us to provide
+    any scaling of the constellation (normalized to sum to 1,
+    normalized so the outside points sit on +/-1, etc.) but still
+    calculate the soft decisions as we would given the full
+    constellation.
+    '''
+    lut_scale = 2**prec
+    maxd = Es*numpy.sqrt(2)/2
+    step = 2*maxd / lut_scale
+    scale = (lut_scale) / (2*maxd) - step
+    
+    xre = sample.real
+    xim = sample.imag
+    xre = int((maxd + min(maxd, max(-maxd, xre))) * scale)
+    xim = int((maxd + min(maxd, max(-maxd, xim))) * scale)
+    index = int(xre + lut_scale*xim)
+
+    max_index = lut_scale**2
+    if(index > max_index):
+        return table[0];
+
+    if(index < 0):
+        raise RuntimeError("calc_from_table: input sample out of range.")
+
+    return table[index]
+
+def calc_soft_dec(sample, constel, symbols, npwr=1):
+    '''
+    This function takes in any consteallation and symbol symbol set
+    (where symbols[i] is the set of bits at constellation point
+    constel[i] and an estimate of the noise power and produces the
+    soft decisions for the given sample.
+
+    If known, the noise power of the received sample may be passed in
+    to this function as npwr.
+
+    This is an incredibly costly algorthm because it must calculate
+    the Euclidean distance between the sample and all points in the
+    constellation to build up its probability
+    calculations. Conversely, it should work for any given
+    constellation/symbol map.
+
+    The function returns a vector of k soft decisions. Decisions less
+    than 0 are more likely to indicate a '0' bit and decisions greater
+    than 0 are more likely to indicate a '1' bit.
+    '''
+    
+    M = len(constel)
+    k = int(numpy.log2(M))
+    tmp = 2*k*[0,]
+    s = k*[0,]
+
+    # Find a scaling factor for the constellation, however it was normalized.
+    constel = numpy.array(constel)
+    scale = min(min(abs(constel.real)), min(abs(constel.imag)))
+
+    for i in range(M):
+        # Calculate the distance between the sample and the current
+        # constellation point.
+        dist = abs(sample - constel[i])**2
+
+        # Calculate the probability factor from the distance and the
+        # scaled noise power.
+        d = numpy.exp(-dist/(2*npwr*scale**2))
+        
+        for j in range(k):
+            # Get the bit at the jth index
+            mask = 1<<j
+            bit = (symbols[i] & mask) >> j
+
+            # If the bit is a 0, add to the probability of a zero
+            if(bit == 0):
+                tmp[2*j+0] += d
+            # else, add to the probability of a one
+            else:
+                tmp[2*j+1] += d
+                
+    # Calculate the log-likelihood ratio for all bits based on the
+    # probability of ones (tmp[2*i+1]) over the probability of a zero
+    # (tmp[2*i+0]).
+    for i in range(k):
+        s[k-1-i] = (numpy.log(tmp[2*i+1]) - numpy.log(tmp[2*i+0])) * scale**2
+
+    return s
diff --git a/gr-digital/python/digital/test_soft_decisions.py b/gr-digital/python/digital/test_soft_decisions.py
new file mode 100755
index 0000000000..78714100b7
--- /dev/null
+++ b/gr-digital/python/digital/test_soft_decisions.py
@@ -0,0 +1,135 @@
+#!/usr/bin/env python
+#
+# Copyright 2013 Free Software Foundation, Inc.
+# 
+# This file is part of GNU Radio
+# 
+# GNU Radio is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3, or (at your option)
+# any later version.
+# 
+# GNU Radio is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with GNU Radio; see the file COPYING.  If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street,
+# Boston, MA 02110-1301, USA.
+# 
+
+import numpy, pylab, sys
+from gnuradio import digital
+from soft_dec_lut_gen import *
+from psk_constellations import *
+from qam_constellations import *
+
+def test_qpsk(i, sample, prec):
+    qpsk_const_list = [psk_4_0, psk_4_1, psk_4_2, psk_4_3,
+                       psk_4_4, psk_4_5, psk_4_6, psk_4_7]
+    qpsk_lut_gen_list = [sd_psk_4_0, sd_psk_4_1, sd_psk_4_2, sd_psk_4_3,
+                         sd_psk_4_4, sd_psk_4_5, sd_psk_4_6, sd_psk_4_7]
+
+    constel, code = qpsk_const_list[i]()
+    qpsk_lut_gen = qpsk_lut_gen_list[i]
+
+    rot_sym = 1
+    side = 2
+    width = 2
+    c = digital.constellation_rect(constel, code, rot_sym,
+                                   side, side, width, width)
+
+    # Get max energy/symbol in constellation
+    constel = c.points()
+    Es = max([numpy.sqrt(constel_i.real**2 + constel_i.imag**2) for constel_i in constel])
+
+    #table = soft_dec_table_generator(qpsk_lut_gen, prec, Es)
+    table = soft_dec_table(constel, code, prec)
+
+    c.gen_soft_dec_lut(prec)
+    #c.set_soft_dec_lut(table, prec)
+
+    y_python_gen_calc = qpsk_lut_gen(sample, Es)
+    y_python_table = calc_soft_dec_from_table(sample, table, prec, Es)
+    y_python_raw_calc = calc_soft_dec(sample, constel, code)
+    y_cpp_table = c.soft_decision_maker(sample)
+    y_cpp_raw_calc = c.calc_soft_dec(sample)
+
+    return (y_python_gen_calc, y_python_table, y_python_raw_calc,
+            y_cpp_table, y_cpp_raw_calc, constel, code, c)
+
+def test_qam16(i, sample, prec):
+    sample = sample/1
+    qam_const_list = [qam_16_0, ]
+    qam_lut_gen_list = [sd_qam_16_0, ]
+
+    constel, code = qam_const_list[i]()
+    qam_lut_gen = qam_lut_gen_list[i]
+
+    rot_sym = 4
+    side = 2
+    width = 2
+    c = digital.constellation_rect(constel, code, rot_sym,
+                                   side, side, width, width)
+
+    # Get max energy/symbol in constellation
+    constel = c.points()
+    Es = max([abs(constel_i) for constel_i in constel])
+
+    #table = soft_dec_table_generator(qam_lut_gen, prec, Es)
+    table = soft_dec_table(constel, code, prec, 1)
+
+    #c.gen_soft_dec_lut(prec)
+    c.set_soft_dec_lut(table, prec)
+
+    y_python_gen_calc = qam_lut_gen(sample, Es)
+    y_python_table = calc_soft_dec_from_table(sample, table, prec, Es)
+    y_python_raw_calc = calc_soft_dec(sample, constel, code, 1)
+    y_cpp_table = c.soft_decision_maker(sample)
+    y_cpp_raw_calc = c.calc_soft_dec(sample)
+
+    return (y_python_gen_calc, y_python_table, y_python_raw_calc,
+            y_cpp_table, y_cpp_raw_calc, constel, code, c)
+
+if __name__ == "__main__":
+
+    index = 0
+    prec = 8
+
+    x_re = 2*numpy.random.random()-1
+    x_im = 2*numpy.random.random()-1
+    x = x_re + x_im*1j
+    #x = -1 + -0.j
+
+    if 1:
+        y_python_gen_calc, y_python_table, y_python_raw_calc, \
+            y_cpp_table, y_cpp_raw_calc, constel, code, c \
+            = test_qpsk(index, x, prec)
+    else:
+        y_python_gen_calc, y_python_table, y_python_raw_calc, \
+            y_cpp_table, y_cpp_raw_calc, constel, code, c \
+            = test_qam16(index, x, prec)
+
+    k = numpy.log2(len(constel))
+
+    print "Sample: ", x
+    print "Python Generator Calculated: ", (y_python_gen_calc)
+    print "Python Generator Table:      ", (y_python_table)
+    print "Python Raw calc:             ", (y_python_raw_calc)
+    print "C++ Table calc:              ", (y_cpp_table)
+    print "C++ Raw calc:                ", (y_cpp_raw_calc)
+
+    fig = pylab.figure(1)
+    sp1 = fig.add_subplot(1,1,1)
+    sp1.plot([c.real for c in constel],
+             [c.imag for c in constel], 'bo')
+    sp1.plot(x.real, x.imag, 'ro')
+    sp1.set_xlim([-1.5, 1.5])
+    sp1.set_ylim([-1.5, 1.5])
+    fill = int(numpy.log2(len(constel)))
+    for i,c in enumerate(constel):
+        sp1.text(1.2*c.real, 1.2*c.imag, bin(code[i])[2:].zfill(fill),
+                 ha='center', va='center', size=18)
+    pylab.show()
author	Tom Rondeau <tom@trondeau.com>	2013-09-04 15:44:20 -0400
committer	Tom Rondeau <tom@trondeau.com>	2013-09-04 15:44:20 -0400
commit	02fd90029a1ec465f7c8da4a0d97dcdb8b3f438f (patch)
tree	9bc144a9cb5b0e877108302d13181e50bd885528
parent	edf2ac24c9c2f723eb80918de7addf6d299958c2 (diff)