diff options
| -rw-r--r-- | gr-digital/include/gnuradio/digital/constellation.h | 45 | ||||
| -rw-r--r-- | gr-digital/lib/constellation.cc | 145 | ||||
| -rw-r--r-- | gr-digital/lib/constellation_soft_decoder_cf_impl.cc | 8 | ||||
| -rw-r--r-- | gr-digital/python/digital/qa_constellation.py | 107 | ||||
| -rw-r--r-- | gr-digital/python/digital/qa_constellation_soft_decoder_cf.py | 153 | ||||
| -rwxr-xr-x | gr-digital/python/digital/qam_constellations.py | 2 | ||||
| -rw-r--r-- | gr-digital/python/digital/soft_dec_lut_gen.py | 96 |
7 files changed, 350 insertions, 206 deletions
diff --git a/gr-digital/include/gnuradio/digital/constellation.h b/gr-digital/include/gnuradio/digital/constellation.h index f2b796ea7f..0c29eff13d 100644 --- a/gr-digital/include/gnuradio/digital/constellation.h +++ b/gr-digital/include/gnuradio/digital/constellation.h @@ -39,7 +39,7 @@ namespace gr { /* */ /* Base class defining interface. */ /************************************************************/ - + class constellation; typedef boost::shared_ptr<constellation> constellation_sptr; @@ -82,13 +82,13 @@ namespace gr { unsigned int decision_maker_pe(const gr_complex *sample, float *phase_error); //! Calculates distance. //unsigned int decision_maker_e(const gr_complex *sample, float *error); - + //! Calculates metrics for all points in the constellation. //! For use with the viterbi algorithm. virtual void calc_metric(const gr_complex *sample, float *metric, gr::digital::trellis_metric_type_t type); virtual void calc_euclidean_metric(const gr_complex *sample, float *metric); virtual void calc_hard_symbol_metric(const gr_complex *sample, float *metric); - + //! Returns the set of points in this constellation. std::vector<gr_complex> points() { return d_constellation;} //! Returns the vector of points in this constellation. @@ -111,7 +111,7 @@ namespace gr { { return floor(log(double(d_constellation.size()))/d_dimensionality/log(2.0)); } - + unsigned int arity() { return d_arity; @@ -120,7 +120,7 @@ namespace gr { constellation_sptr base() { return shared_from_this(); - } + } pmt::pmt_t as_pmt() { @@ -185,7 +185,12 @@ namespace gr { //! Returns True if the soft decision LUT has been defined, False otherwise. bool has_soft_dec_lut(); - + + + std::vector< std::vector<float> > soft_dec_lut(); + + + /*! \brief Returns the soft decisions for the given \p sample. * * \details Returns the soft decisions for the given \p @@ -200,7 +205,7 @@ namespace gr { protected: std::vector<gr_complex> d_constellation; - std::vector<int> d_pre_diff_code; + std::vector<int> d_pre_diff_code; bool d_apply_pre_diff_code; unsigned int d_rotational_symmetry; unsigned int d_dimensionality; @@ -244,7 +249,7 @@ namespace gr { * Make a general constellation object that calculates the Euclidean distance for hard decisions. * * \param constell List of constellation points (order of list matches pre_diff_code) - * \param pre_diff_code List of alphabet symbols (before applying any differential + * \param pre_diff_code List of alphabet symbols (before applying any differential * coding) (order of list matches constell) * \param rotational_symmetry Number of rotations around unit circle that have the same representation. * \param dimensionality Number of dimensions to the constellation. @@ -287,7 +292,7 @@ namespace gr { * Make a sectorized constellation object. * * \param constell List of constellation points (order of list matches pre_diff_code) - * \param pre_diff_code List of alphabet symbols (before applying any differential + * \param pre_diff_code List of alphabet symbols (before applying any differential * coding) (order of list matches constell) * \param rotational_symmetry Number of rotations around unit circle that have the same representation. * \param dimensionality Number of z-axis dimensions to the constellation @@ -342,13 +347,13 @@ namespace gr { * Make a rectangular constellation object. * * \param constell List of constellation points (order of list matches pre_diff_code) - * \param pre_diff_code List of alphabet symbols (before applying any differential + * \param pre_diff_code List of alphabet symbols (before applying any differential * coding) (order of list matches constell) * \param rotational_symmetry Number of rotations around unit circle that have the same representation. * \param real_sectors Number of sectors the real axis is split in to. * \param imag_sectors Number of sectors the imag axis is split in to. * \param width_real_sectors width of each real sector to calculate decision boundaries. - * \param width_imag_sectors width of each imag sector to calculate decision boundaries. + * \param width_imag_sectors width of each imag sector to calculate decision boundaries. */ static constellation_rect::sptr make(std::vector<gr_complex> constell, std::vector<int> pre_diff_code, @@ -402,7 +407,7 @@ namespace gr { * the closest constellation point, however sometimes it's nice to * have the flexibility. */ - class DIGITAL_API constellation_expl_rect + class DIGITAL_API constellation_expl_rect : public constellation_rect { public: @@ -440,17 +445,17 @@ namespace gr { /* constellation_psk */ /************************************************************/ - /*! + /*! * \brief constellation_psk * \ingroup digital * * Constellation space is divided into pie slices sectors. * - * Each slice is associated with the nearest constellation point. + * Each slice is associated with the nearest constellation point. * - * Works well for PSK but nothing else. + * Works well for PSK but nothing else. * - * Assumes that there is a constellation point at 1.x + * Assumes that there is a constellation point at 1.x */ class DIGITAL_API constellation_psk : public constellation_sector { @@ -466,7 +471,7 @@ namespace gr { protected: unsigned int get_sector(const gr_complex *sample); - + unsigned int calc_sector_value(unsigned int sector); constellation_psk(std::vector<gr_complex> constell, @@ -482,7 +487,7 @@ namespace gr { /* */ /************************************************************/ - /*! + /*! * \brief Digital constellation for BPSK . * \ingroup digital * @@ -515,7 +520,7 @@ namespace gr { /* */ /************************************************************/ - /*! + /*! * \brief Digital constellation for QPSK * \ingroup digital * @@ -585,7 +590,7 @@ namespace gr { /* */ /************************************************************/ - /*! + /*! * \brief Digital constellation for 8PSK. * \ingroup digital * diff --git a/gr-digital/lib/constellation.cc b/gr-digital/lib/constellation.cc index 7d2c170d34..b17635fec1 100644 --- a/gr-digital/lib/constellation.cc +++ b/gr-digital/lib/constellation.cc @@ -32,6 +32,7 @@ #include <cfloat> #include <stdexcept> #include <boost/format.hpp> +#include <iostream> namespace gr { namespace digital { @@ -129,7 +130,7 @@ namespace gr { unsigned int min_index = 0; float min_euclid_dist; float euclid_dist; - + min_euclid_dist = get_distance(0, sample); min_index = 0; for(unsigned int j = 1; j < d_arity; j++){ @@ -244,21 +245,22 @@ 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) { + d_lut_scale = powf(2.0f, static_cast<float>(precision)); + + // We know we've normalized the constellation, so the min/max + // dimensions in either direction are scaled to +/-1. + float maxd = 1.0f; + float step = (2.0f*maxd) / (d_lut_scale-1); + float y = -maxd; + while(y < maxd+step) { + float x = -maxd; + while(x < maxd+step) { gr_complex pt = gr_complex(x, y); d_soft_dec_lut.push_back(calc_soft_dec(pt, npwr)); - x += xstep; + x += step; } - y += ystep; + y += step; } d_lut_precision = precision; @@ -267,22 +269,19 @@ namespace gr { std::vector<float> constellation::calc_soft_dec(gr_complex sample, float npwr) { + int v; int M = static_cast<int>(d_constellation.size()); int k = static_cast<int>(log(static_cast<double>(M))/log(2.0)); std::vector<float> tmp(2*k, 0); std::vector<float> s(k, 0); - float scale = d_scalefactor*d_scalefactor; - int v; - 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); - + float dist = std::abs(sample - d_constellation[i]); // Calculate the probability factor from the distance and // the scaled noise power. - float d = expf(-dist / (2.0*npwr*scale)); + float d = expf(-dist/npwr); if(d_apply_pre_diff_code) v = d_pre_diff_code[i]; @@ -307,7 +306,7 @@ namespace gr { // 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; + s[k-1-i] = (logf(tmp[2*i+1]) - logf(tmp[2*i+0])); } return s; @@ -330,30 +329,40 @@ namespace gr { return d_soft_dec_lut.size() > 0; } + std::vector< std::vector<float> > + constellation::soft_dec_lut() + { + return d_soft_dec_lut; + } + 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); + // Clip to just below 1 --> at 1, we can overflow the index + // that will put us in the next row of the 2D LUT. + float xre = branchless_clip(sample.real(), 0.99); + float xim = branchless_clip(sample.imag(), 0.99); + + // We normalize the constellation in the ctor, so we know that + // the maximum dimenions go from -1 to +1. We can infer the x + // and y scale directly. + float scale = d_lut_scale / (2.0f); + + // Convert the clipped x and y samples to nearest index offset + xre = floorf((1.0f + xre) * scale); + xim = floorf((1.0f + xim) * scale); 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."); + // Make sure we are in bounds of the index + while(index >= max_index) { + index -= d_lut_scale; + } + while(index < 0) { + index += d_lut_scale; + } return d_soft_dec_lut[index]; } @@ -412,7 +421,7 @@ namespace gr { unsigned int dimensionality) : constellation(constell, pre_diff_code, rotational_symmetry, dimensionality) {} - + // Chooses points base on shortest distance. // Inefficient. unsigned int @@ -435,11 +444,11 @@ namespace gr { n_sectors(n_sectors) { } - + constellation_sector::~constellation_sector() { } - + unsigned int constellation_sector::decision_maker(const gr_complex *sample) { @@ -447,7 +456,7 @@ namespace gr { sector = get_sector(sample); return sector_values[sector]; } - + void constellation_sector::find_sector_values() { @@ -457,11 +466,11 @@ namespace gr { sector_values.push_back(calc_sector_value(i)); } } - - + + /********************************************************************/ - - + + constellation_rect::sptr constellation_rect::make(std::vector<gr_complex> constell, std::vector<int> pre_diff_code, @@ -495,31 +504,31 @@ namespace gr { d_width_imag_sectors *= d_scalefactor; find_sector_values(); } - + constellation_rect::~constellation_rect() { } - + unsigned int constellation_rect::get_sector(const gr_complex *sample) { int real_sector, imag_sector; unsigned int sector; - + real_sector = int(real(*sample)/d_width_real_sectors + n_real_sectors/2.0); if(real_sector < 0) real_sector = 0; if(real_sector >= (int)n_real_sectors) real_sector = n_real_sectors-1; - + imag_sector = int(imag(*sample)/d_width_imag_sectors + n_imag_sectors/2.0); if(imag_sector < 0) imag_sector = 0; if(imag_sector >= (int)n_imag_sectors) imag_sector = n_imag_sectors-1; - + sector = real_sector * n_imag_sectors + imag_sector; return sector; } @@ -536,7 +545,7 @@ namespace gr { (imag_sector + 0.5 - n_imag_sectors/2.0) * d_width_imag_sectors); return sector_center; } - + unsigned int constellation_rect::calc_sector_value(unsigned int sector) { @@ -547,8 +556,8 @@ namespace gr { } /********************************************************************/ - - constellation_expl_rect::sptr + + constellation_expl_rect::sptr constellation_expl_rect::make(std::vector<gr_complex> constellation, std::vector<int> pre_diff_code, unsigned int rotational_symmetry, @@ -565,7 +574,7 @@ namespace gr { width_real_sectors, width_imag_sectors, sector_values)); } - + constellation_expl_rect::constellation_expl_rect( std::vector<gr_complex> constellation, std::vector<int> pre_diff_code, @@ -580,16 +589,16 @@ namespace gr { d_sector_values(sector_values) { } - + constellation_expl_rect::~constellation_expl_rect() { } - + /********************************************************************/ - - - constellation_psk::sptr - constellation_psk::make(std::vector<gr_complex> constell, + + + constellation_psk::sptr + constellation_psk::make(std::vector<gr_complex> constell, std::vector<int> pre_diff_code, unsigned int n_sectors) { @@ -597,7 +606,7 @@ namespace gr { (constell, pre_diff_code, n_sectors)); } - + constellation_psk::constellation_psk(std::vector<gr_complex> constell, std::vector<int> pre_diff_code, unsigned int n_sectors) : @@ -606,11 +615,11 @@ namespace gr { { find_sector_values(); } - + constellation_psk::~constellation_psk() { } - + unsigned int constellation_psk::get_sector(const gr_complex *sample) { @@ -621,7 +630,7 @@ namespace gr { sector += n_sectors; return sector; } - + unsigned int constellation_psk::calc_sector_value(unsigned int sector) { @@ -635,7 +644,7 @@ namespace gr { /********************************************************************/ - constellation_bpsk::sptr + constellation_bpsk::sptr constellation_bpsk::make() { return constellation_bpsk::sptr(new constellation_bpsk()); @@ -665,7 +674,7 @@ namespace gr { /********************************************************************/ - constellation_qpsk::sptr + constellation_qpsk::sptr constellation_qpsk::make() { return constellation_qpsk::sptr(new constellation_qpsk()); @@ -679,7 +688,7 @@ namespace gr { d_constellation[1] = gr_complex(SQRT_TWO, -SQRT_TWO); d_constellation[2] = gr_complex(-SQRT_TWO, SQRT_TWO); d_constellation[3] = gr_complex(SQRT_TWO, SQRT_TWO); - + /* d_constellation[0] = gr_complex(SQRT_TWO, SQRT_TWO); d_constellation[1] = gr_complex(-SQRT_TWO, SQRT_TWO); @@ -701,7 +710,7 @@ namespace gr { constellation_qpsk::~constellation_qpsk() { } - + unsigned int constellation_qpsk::decision_maker(const gr_complex *sample) { @@ -731,7 +740,7 @@ namespace gr { /********************************************************************/ - constellation_dqpsk::sptr + constellation_dqpsk::sptr constellation_dqpsk::make() { return constellation_dqpsk::sptr(new constellation_dqpsk()); @@ -791,7 +800,7 @@ namespace gr { /********************************************************************/ - constellation_8psk::sptr + constellation_8psk::sptr constellation_8psk::make() { return constellation_8psk::sptr(new constellation_8psk()); diff --git a/gr-digital/lib/constellation_soft_decoder_cf_impl.cc b/gr-digital/lib/constellation_soft_decoder_cf_impl.cc index cd0a844119..aeab660746 100644 --- a/gr-digital/lib/constellation_soft_decoder_cf_impl.cc +++ b/gr-digital/lib/constellation_soft_decoder_cf_impl.cc @@ -1,19 +1,19 @@ /* -*- c++ -*- */ /* * 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, diff --git a/gr-digital/python/digital/qa_constellation.py b/gr-digital/python/digital/qa_constellation.py index 00248f4f09..42e49bb059 100644 --- a/gr-digital/python/digital/qa_constellation.py +++ b/gr-digital/python/digital/qa_constellation.py @@ -1,26 +1,26 @@ #!/usr/bin/env python # # Copyright 2011,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 random +import random, math from cmath import exp, pi, log, sqrt from gnuradio import gr, gr_unittest, digital, blocks @@ -40,7 +40,7 @@ tested_mod_codes = (mod_codes.NO_CODE, mod_codes.GRAY_CODE) def twod_constell(): """ - + """ points = ((1+0j), (0+1j), (-1+0j), (0-1j)) @@ -76,12 +76,12 @@ easy_constellation_info = ( False, None), (qam.qam_constellation, {'constellation_points': (4,), - 'mod_code': tested_mod_codes, + 'mod_code': tested_mod_codes, 'large_ampls_to_corners': [False],}, True, None), (qam.qam_constellation, {'constellation_points': (4, 16, 64), - 'mod_code': tested_mod_codes, + 'mod_code': tested_mod_codes, 'differential': (False,)}, False, None), (digital.constellation_bpsk, {}, True, None), @@ -101,7 +101,7 @@ medium_constellation_info = ( True, None), (qam.qam_constellation, {'constellation_points': (16 ,), - 'mod_code': tested_mod_codes, + 'mod_code': tested_mod_codes, 'large_ampls_to_corners': [False, True],}, True, None), (qamlike.qam32_holeinside_constellation, @@ -113,11 +113,20 @@ medium_constellation_info = ( difficult_constellation_info = ( (qam.qam_constellation, {'constellation_points': (64,), - 'mod_code': tested_mod_codes, + 'mod_code': tested_mod_codes, 'large_ampls_to_corners': [False, True],}, - True, None), + True, None), ) +def slicer(x): + ret = [] + for xi in x: + if(xi < 0): + ret.append(0.0) + else: + ret.append(1.0) + return ret + def tested_constellations(easy=True, medium=True, difficult=True): """ Generator to produce (constellation, differential) tuples for testing purposes. @@ -153,7 +162,7 @@ def tested_constellations(easy=True, medium=True, difficult=True): this_poss_arg[2] = 0 if sum([argindex for argname, argvalues, argindex in poss_args]) == 0: break - + class test_constellation(gr_unittest.TestCase): @@ -170,7 +179,7 @@ class test_constellation(gr_unittest.TestCase): for constellation, differential in tested_constellations(): if differential: rs = constellation.rotational_symmetry() - rotations = [exp(i*2*pi*(0+1j)/rs) for i in range(0, rs)] + rotations = [exp(i*2*pi*(0+1j)/rs) for i in range(0, rs)] else: rotations = [None] for rotation in rotations: @@ -216,17 +225,17 @@ class test_constellation(gr_unittest.TestCase): 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_python_raw_calc += slicer(digital.calc_soft_dec(sample, constel, code)) + y_python_gen_calc += slicer(digital.sd_psk_4_0(sample, Es)) + y_python_table += slicer(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) - + self.assertFloatTuplesAlmostEqual(y_python_raw_calc, y_python_gen_calc, 0) + self.assertFloatTuplesAlmostEqual(y_python_gen_calc, y_python_table, 0) + self.assertFloatTuplesAlmostEqual(y_cpp_raw_calc, y_cpp_table, 0) + def test_soft_qpsk_calc(self): prec = 8 constel, code = digital.psk_4_0() @@ -257,14 +266,54 @@ class test_constellation(gr_unittest.TestCase): 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_python_raw_calc += slicer(digital.calc_soft_dec(sample, constel, code)) + y_python_table += slicer(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) + y_cpp_raw_calc += slicer(c.calc_soft_dec(sample)) + y_cpp_table += slicer(c.soft_decision_maker(sample)) + + self.assertEqual(y_python_raw_calc, y_python_table) + self.assertEqual(y_cpp_raw_calc, y_cpp_table) + + + def test_soft_qam16_calc(self): + prec = 8 + constel, code = digital.qam_16_0() - self.assertFloatTuplesAlmostEqual(y_python_raw_calc, y_python_table, 4) - self.assertFloatTuplesAlmostEqual(y_cpp_raw_calc, y_cpp_table, 4) + 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 += slicer(digital.calc_soft_dec(sample, constel, code)) + y_python_table += slicer(digital.calc_soft_dec_from_table(sample, table, prec, Es)) + + y_cpp_raw_calc += slicer(c.calc_soft_dec(sample)) + y_cpp_table += slicer(c.soft_decision_maker(sample)) + + self.assertFloatTuplesAlmostEqual(y_python_raw_calc, y_python_table, 0) + self.assertFloatTuplesAlmostEqual(y_cpp_raw_calc, y_cpp_table, 0) class mod_demod(gr.hier_block2): def __init__(self, constellation, differential, rotation): @@ -317,7 +366,7 @@ class mod_demod(gr.hier_block2): if self.constellation.apply_pre_diff_code(): self.blocks.append(digital.map_bb( mod_codes.invert_code(self.constellation.pre_diff_code()))) - # unpack the k bit vector into a stream of bits + # unpack the k bit vector into a stream of bits self.blocks.append(blocks.unpack_k_bits_bb( self.constellation.bits_per_symbol())) # connect to block output @@ -326,6 +375,6 @@ class mod_demod(gr.hier_block2): self.blocks.append(weakref.proxy(self)) self.connect(*self.blocks) - + if __name__ == '__main__': gr_unittest.run(test_constellation, "test_constellation.xml") diff --git a/gr-digital/python/digital/qa_constellation_soft_decoder_cf.py b/gr-digital/python/digital/qa_constellation_soft_decoder_cf.py index 6cd757537a..448e76e834 100644 --- a/gr-digital/python/digital/qa_constellation_soft_decoder_cf.py +++ b/gr-digital/python/digital/qa_constellation_soft_decoder_cf.py @@ -1,27 +1,28 @@ #!/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. -# +# from gnuradio import gr, gr_unittest, digital, blocks from math import sqrt +from numpy import random, vectorize class test_constellation_soft_decoder(gr_unittest.TestCase): @@ -31,16 +32,15 @@ class test_constellation_soft_decoder(gr_unittest.TestCase): def tearDown(self): self.tb = None - def test_constellation_soft_decoder_cf_bpsk(self): - prec = 8 - src_data = (0.5 + 0.5j, 0.1 - 1.2j, -0.8 - 0.1j, -0.45 + 0.8j, - 0.8 + 1.0j, -0.5 + 0.1j, 0.1 - 1.2j, 1+1j) - lut = digital.soft_dec_table_generator(digital.sd_psk_2_0x0, prec) + def helper_with_lut(self, prec, src_data, const_gen, const_sd_gen): + cnst_pts, code = const_gen() + Es = max([abs(c) for c in cnst_pts]) + lut = digital.soft_dec_table_generator(const_sd_gen, prec, Es) expected_result = list() for s in src_data: - expected_result += digital.calc_soft_dec_from_table(s, lut, prec, Es=2/sqrt(2)) + res = digital.calc_soft_dec_from_table(s, lut, prec, sqrt(2.0)) + expected_result += res - cnst_pts, code = digital.psk_2_0x0() cnst = digital.constellation_calcdist(cnst_pts, code, 2, 1) cnst.set_soft_dec_lut(lut, int(prec)) src = blocks.vector_source_c(src_data) @@ -49,64 +49,123 @@ class test_constellation_soft_decoder(gr_unittest.TestCase): self.tb.connect(src, op) self.tb.connect(op, dst) - self.tb.run() # run the graph and wait for it to finish + self.tb.run() actual_result = dst.data() # fetch the contents of the sink #print "actual result", actual_result #print "expected result", expected_result - self.assertFloatTuplesAlmostEqual(expected_result, actual_result, 4) + self.assertFloatTuplesAlmostEqual(expected_result, actual_result, 5) - def test_constellation_soft_decoder_cf_qpsk(self): - prec = 8 - src_data = (0.5 + 0.5j, 0.1 - 1.2j, -0.8 - 0.1j, -0.45 + 0.8j, - 0.8 + 1.0j, -0.5 + 0.1j, 0.1 - 1.2j, 1+1j) - lut = digital.soft_dec_table_generator(digital.sd_psk_4_0x0_0_1, prec) + def helper_no_lut(self, prec, src_data, const_gen, const_sd_gen): + cnst_pts, code = const_gen() + cnst = digital.constellation_calcdist(cnst_pts, code, 2, 1) expected_result = list() for s in src_data: - expected_result += digital.calc_soft_dec_from_table(s, lut, prec) + res = digital.calc_soft_dec(s, cnst.points(), code) + expected_result += res - cnst_pts,code = digital.psk_4_0x0_0_1() - cnst = digital.constellation_calcdist(cnst_pts, code, 2, 1) - cnst.set_soft_dec_lut(lut, int(prec)) src = blocks.vector_source_c(src_data) op = digital.constellation_soft_decoder_cf(cnst.base()) dst = blocks.vector_sink_f() self.tb.connect(src, op) self.tb.connect(op, dst) - self.tb.run() # run the graph and wait for it to finish + self.tb.run() actual_result = dst.data() # fetch the contents of the sink #print "actual result", actual_result #print "expected result", expected_result - self.assertFloatTuplesAlmostEqual(expected_result, actual_result, 5) - def test_constellation_soft_decoder_cf_qam16(self): + # Double vs. float precision issues between Python and C++, so + # use only 4 decimals in comparisons. + self.assertFloatTuplesAlmostEqual(expected_result, actual_result, 4) + + def test_constellation_soft_decoder_cf_bpsk_3(self): + prec = 3 + src_data = (-1.0 - 1.0j, 1.0 - 1.0j, -1.0 + 1.0j, 1.0 + 1.0j, + -2.0 - 2.0j, 2.0 - 2.0j, -2.0 + 2.0j, 2.0 + 2.0j, + -0.2 - 0.2j, 0.2 - 0.2j, -0.2 + 0.2j, 0.2 + 0.2j, + 0.3 + 0.4j, 0.1 - 1.2j, -0.8 - 0.1j, -0.4 + 0.8j, + 0.8 + 1.0j, -0.5 + 0.1j, 0.1 + 1.2j, -1.7 - 0.9j) + self.helper_with_lut(prec, src_data, digital.psk_2_0x0, digital.sd_psk_2_0x0) + + def test_constellation_soft_decoder_cf_bpsk_8(self): prec = 8 - src_data = (0.5 + 0.5j, 0.1 - 1.2j, -0.8 - 0.1j, -0.45 + 0.8j, - 0.8 + 1.0j, -0.5 + 0.1j, 0.1 - 1.2j, 1+1j) - lut = digital.soft_dec_table_generator(digital.sd_qam_16_0x0_0_1_2_3, prec) - expected_result = list() - for s in src_data: - expected_result += digital.calc_soft_dec_from_table(s, lut, prec) + src_data = (-1.0 - 1.0j, 1.0 - 1.0j, -1.0 + 1.0j, 1.0 + 1.0j, + -2.0 - 2.0j, 2.0 - 2.0j, -2.0 + 2.0j, 2.0 + 2.0j, + -0.2 - 0.2j, 0.2 - 0.2j, -0.2 + 0.2j, 0.2 + 0.2j, + 0.3 + 0.4j, 0.1 - 1.2j, -0.8 - 0.1j, -0.4 + 0.8j, + 0.8 + 1.0j, -0.5 + 0.1j, 0.1 + 1.2j, -1.7 - 0.9j) + self.helper_with_lut(prec, src_data, digital.psk_2_0x0, digital.sd_psk_2_0x0) + + def test_constellation_soft_decoder_cf_bpsk_8_rand(self): + prec = 8 + src_data = vectorize(complex)(2*random.randn(100), 2*random.randn(100)) + self.helper_with_lut(prec, src_data, digital.psk_2_0x0, digital.sd_psk_2_0x0) - cnst_pts = digital.qam_16_0x0_0_1_2_3() - cnst = digital.constellation_calcdist(cnst_pts[0], cnst_pts[1], 2, 1) - cnst.set_soft_dec_lut(lut, int(prec)) - src = blocks.vector_source_c(src_data) - op = digital.constellation_soft_decoder_cf(cnst.base()) - dst = blocks.vector_sink_f() + def test_constellation_soft_decoder_cf_bpsk_8_rand2(self): + prec = 8 + src_data = vectorize(complex)(2*random.randn(100), 2*random.randn(100)) + self.helper_no_lut(prec, src_data, digital.psk_2_0x0, digital.sd_psk_2_0x0) + + def test_constellation_soft_decoder_cf_qpsk_3(self): + prec = 3 + src_data = (-1.0 - 1.0j, 1.0 - 1.0j, -1.0 + 1.0j, 1.0 + 1.0j, + -2.0 - 2.0j, 2.0 - 2.0j, -2.0 + 2.0j, 2.0 + 2.0j, + -0.2 - 0.2j, 0.2 - 0.2j, -0.2 + 0.2j, 0.2 + 0.2j, + 0.3 + 0.4j, 0.1 - 1.2j, -0.8 - 0.1j, -0.4 + 0.8j, + 0.8 + 1.0j, -0.5 + 0.1j, 0.1 + 1.2j, -1.7 - 0.9j) + self.helper_with_lut(prec, src_data, digital.psk_4_0x0_0_1, digital.sd_psk_4_0x0_0_1) + + def test_constellation_soft_decoder_cf_qpsk_8(self): + prec = 8 + src_data = (-1.0 - 1.0j, 1.0 - 1.0j, -1.0 + 1.0j, 1.0 + 1.0j, + -2.0 - 2.0j, 2.0 - 2.0j, -2.0 + 2.0j, 2.0 + 2.0j, + -0.2 - 0.2j, 0.2 - 0.2j, -0.2 + 0.2j, 0.2 + 0.2j, + 0.3 + 0.4j, 0.1 - 1.2j, -0.8 - 0.1j, -0.4 + 0.8j, + 0.8 + 1.0j, -0.5 + 0.1j, 0.1 + 1.2j, -1.7 - 0.9j) + self.helper_with_lut(prec, src_data, digital.psk_4_0x0_0_1, digital.sd_psk_4_0x0_0_1) + + def test_constellation_soft_decoder_cf_qpsk_8_rand(self): + prec = 8 + src_data = vectorize(complex)(2*random.randn(100), 2*random.randn(100)) + self.helper_with_lut(prec, src_data, digital.psk_4_0x0_0_1, digital.sd_psk_4_0x0_0_1) - self.tb.connect(src, op) - self.tb.connect(op, dst) - self.tb.run() # run the graph and wait for it to finish + def test_constellation_soft_decoder_cf_qpsk_8_rand2(self): + prec = 8 + src_data = vectorize(complex)(2*random.randn(100), 2*random.randn(100)) + self.helper_no_lut(prec, src_data, digital.psk_4_0x0_0_1, digital.sd_psk_4_0x0_0_1) + + def test_constellation_soft_decoder_cf_qam16_3(self): + prec = 3 + src_data = (-1.0 - 1.0j, 1.0 - 1.0j, -1.0 + 1.0j, 1.0 + 1.0j, + -2.0 - 2.0j, 2.0 - 2.0j, -2.0 + 2.0j, 2.0 + 2.0j, + -0.2 - 0.2j, 0.2 - 0.2j, -0.2 + 0.2j, 0.2 + 0.2j, + 0.3 + 0.4j, 0.1 - 1.2j, -0.8 - 0.1j, -0.4 + 0.8j, + 0.8 + 1.0j, -0.5 + 0.1j, 0.1 + 1.2j, -1.7 - 0.9j) + self.helper_with_lut(prec, src_data, digital.qam_16_0x0_0_1_2_3, digital.sd_qam_16_0x0_0_1_2_3) + + def test_constellation_soft_decoder_cf_qam16_8(self): + prec = 8 + src_data = (-1.0 - 1.0j, 1.0 - 1.0j, -1.0 + 1.0j, 1.0 + 1.0j, + -2.0 - 2.0j, 2.0 - 2.0j, -2.0 + 2.0j, 2.0 + 2.0j, + -0.2 - 0.2j, 0.2 - 0.2j, -0.2 + 0.2j, 0.2 + 0.2j, + 0.3 + 0.4j, 0.1 - 1.2j, -0.8 - 0.1j, -0.4 + 0.8j, + 0.8 + 1.0j, -0.5 + 0.1j, 0.1 + 1.2j, -1.7 - 0.9j) + self.helper_with_lut(prec, src_data, digital.qam_16_0x0_0_1_2_3, digital.sd_qam_16_0x0_0_1_2_3) + + def test_constellation_soft_decoder_cf_qam16_8_rand(self): + prec = 8 + src_data = vectorize(complex)(2*random.randn(100), 2*random.randn(100)) + self.helper_with_lut(prec, src_data, digital.qam_16_0x0_0_1_2_3, digital.sd_qam_16_0x0_0_1_2_3) - actual_result = dst.data() # fetch the contents of the sink - #print "actual result", actual_result - #print "expected result", expected_result - self.assertFloatTuplesAlmostEqual(expected_result, actual_result, 5) + def test_constellation_soft_decoder_cf_qam16_8_rand2(self): + prec = 8 + #src_data = vectorize(complex)(2*random.randn(100), 2*random.randn(100)) + src_data = vectorize(complex)(2*random.randn(2), 2*random.randn(2)) + self.helper_no_lut(prec, src_data, digital.qam_16_0x0_0_1_2_3, digital.sd_qam_16_0x0_0_1_2_3) if __name__ == '__main__': - gr_unittest.run(test_constellation_soft_decoder, "test_constellation_soft_decoder.xml") - + #gr_unittest.run(test_constellation_soft_decoder, "test_constellation_soft_decoder.xml") + gr_unittest.run(test_constellation_soft_decoder) diff --git a/gr-digital/python/digital/qam_constellations.py b/gr-digital/python/digital/qam_constellations.py index 6f9f6bfab5..ed86d7e1e3 100755 --- a/gr-digital/python/digital/qam_constellations.py +++ b/gr-digital/python/digital/qam_constellations.py @@ -270,7 +270,7 @@ def sd_qam_16_0x0_0_1_2_3(x, Es=1): b2 = -abs(x_re) + boundary b0 = -abs(x_im) + boundary - return [(Es/2)*b3, (Es/2)*b2, (Es/2)*b1, (Es/2)*b0] + return [(Es/2.0)*b3, (Es/2.0)*b2, (Es/2.0)*b1, (Es/2.0)*b0] sd_qam_16 = sd_qam_16_0x0_0_1_2_3 sd_qam_16_0 = sd_qam_16 diff --git a/gr-digital/python/digital/soft_dec_lut_gen.py b/gr-digital/python/digital/soft_dec_lut_gen.py index 24d8bbdc2e..9c184d60cb 100644 --- a/gr-digital/python/digital/soft_dec_lut_gen.py +++ b/gr-digital/python/digital/soft_dec_lut_gen.py @@ -1,30 +1,29 @@ #!/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 + '''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). @@ -70,18 +69,19 @@ def soft_dec_table_generator(soft_dec_gen, prec, Es=1): (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, + Es is the maximum 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 + maxd = Es*numpy.sqrt(2.0)/2.0 yrng = numpy.linspace(-maxd, maxd, npts) xrng = numpy.linspace(-maxd, maxd, npts) @@ -129,7 +129,7 @@ def soft_dec_table(constel, symbols, prec, npwr=1): table.append(decs) return table -def calc_soft_dec_from_table(sample, table, prec, Es=1): +def calc_soft_dec_from_table(sample, table, prec, Es=1.0): ''' 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 @@ -152,25 +152,25 @@ def calc_soft_dec_from_table(sample, table, prec, Es=1): 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 - + lut_scale = 2.0**prec + maxd = Es*numpy.sqrt(2.0)/2.0 + scale = (lut_scale) / (2.0*maxd) + + alpha = 0.99 # to keep index within bounds 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) + xre = ((maxd + min(alpha*maxd, max(-alpha*maxd, xre))) * scale) + xim = ((maxd + min(alpha*maxd, max(-alpha*maxd, xim))) * scale) + index = int(xre) + lut_scale*int(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.") + while(index >= max_index): + index -= lut_scale; + while(index < 0): + index += lut_scale; - return table[index] + return table[int(index)] def calc_soft_dec(sample, constel, symbols, npwr=1): ''' @@ -192,25 +192,21 @@ def calc_soft_dec(sample, constel, symbols, npwr=1): 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 + dist = abs(sample - constel[i]) # Calculate the probability factor from the distance and the # scaled noise power. - d = numpy.exp(-dist/(2*npwr*scale**2)) - + d = numpy.exp(-dist/npwr) + for j in range(k): # Get the bit at the jth index mask = 1<<j @@ -222,11 +218,37 @@ def calc_soft_dec(sample, constel, symbols, npwr=1): # 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 + s[k-1-i] = (numpy.log(tmp[2*i+1]) - numpy.log(tmp[2*i+0])) return s + + +def show_table(table): + prec = int(numpy.sqrt(len(table))) + nbits = len(table[0]) + pp = "" + subi = 1 + subj = 0 + for i in reversed(xrange(prec+1)): + if(i == prec//2): + pp += "-----" + prec*((nbits*8)+3)*"-" + "\n" + subi = 0 + continue + for j in xrange(prec+1): + if(j == prec//2): + pp += "| " + subj = 1 + else: + item = table[prec*(i-subi) + (j-subj)] + pp += "( " + for t in xrange(nbits-1, -1, -1): + pp += "{0: .4f} ".format(item[t]) + pp += ") " + pp += "\n" + subj = 0 + print pp |