[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Commit-gnuradio] [gnuradio] 01/25: polar: Python test code for polar co
|
From: |
git |
|
Subject: |
[Commit-gnuradio] [gnuradio] 01/25: polar: Python test code for polar codes |
|
Date: |
Wed, 23 Sep 2015 14:51:36 +0000 (UTC) |
This is an automated email from the git hooks/post-receive script.
jcorgan pushed a commit to branch master
in repository gnuradio.
commit 8ac5a86f98d5a545cf65dda77c960cf273ee7aca
Author: Johannes Demel <address@hidden>
Date: Thu May 28 12:39:36 2015 +0200
polar: Python test code for polar codes
---
.gitignore | 3 +
README | 139 ---
README.md | 23 +
gr-fec/python/fec/polar/README.md | 4 +
.../polar/bit_reversed_polar_encoding_scheme.svg | 1160 ++++++++++++++++++++
.../python/fec/polar/channel_construction_bec.py | 110 ++
.../python/fec/polar/channel_construction_bsc.py | 206 ++++
gr-fec/python/fec/polar/common.py | 69 ++
gr-fec/python/fec/polar/decoder.py | 278 +++++
gr-fec/python/fec/polar/encoder.py | 135 +++
gr-fec/python/fec/polar/helper_functions.py | 84 ++
gr-fec/python/fec/polar/polar_code.py | 493 +++++++++
gr-fec/python/fec/polar/testbed.py | 143 +++
13 files changed, 2708 insertions(+), 139 deletions(-)
diff --git a/.gitignore b/.gitignore
index a185998..5140c03 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,6 @@
*.pyc
*.pyo
build/
+.idea/
+*.npy
+*.lprof
diff --git a/README b/README
deleted file mode 100644
index f460dd7..0000000
--- a/README
+++ /dev/null
@@ -1,139 +0,0 @@
-#
-# Copyright 2001-2007,2009,2012 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.
-#
-
-Welcome to GNU Radio!
-
-
-Please see http://gnuradio.org for the wiki, bug tracking,
-and source code viewer.
-
-If you've got questions about GNU Radio, please subscribe to the
-discuss-gnuradio mailing list and post your questions there.
-http://gnuradio.org/redmine/projects/gnuradio/wiki/MailingLists
-
-There is also a "Build Guide" in the wiki that contains OS specific
-recommendations:
-http://gnuradio.org/redmine/projects/gnuradio/wiki/BuildGuide
-
-
-The bleeding edge code can be found in our git repository at
-http://gnuradio.org/git/gnuradio.git/. To checkout the latest, use
-this command:
-
- $ git clone git://git.gnuradio.org/gnuradio
-
-For information about using Git, please see:
-http://gnuradio.org/redmine/projects/gnuradio/wiki/DevelopingWithGit
-
-
-How to Build GNU Radio:
-
-For more complete instructions, see the "Building GNU Radio" page in
-the GNU Radio manual (can be built or found online at
-http://gnuradio.org/doc/doxygen/build_guide.html).
-
-See these steps for a quick build guide.
-
- (1) Ensure that you've satisfied the external dependencies. These
- dependencies are listed in the manual's build page and are not
- presented here to reduce duplication errors.
-
- See the wiki at http://gnuradio.org for details.
-
-
- (2) Building from cmake:
-
- $ mkdir $(builddir)
- $ cd $(builddir)
- $ cmake [OPTIONS] $(srcdir)
- $ make
- $ make test
- $ sudo make install
-
-
-That's it!
-
-Options:
-Useful options include setting the install prefix and the build type:
-
- -DCMAKE_INSTALL_PREFIX=<directory to install to>
- -DCMAKE_BUILD_TYPE="<type>"
-
-Currently, GNU Radio has a "Debug" type that builds with '-g -O2'
-useful for debugging the software and a "Release" type that builds
-with '-O3', which is the default.
-
-
--------------------------------------------------------------------------------
-
- KNOWN INCOMPATIBILITIES
-
-
- GNU Radio triggers bugs in g++ 3.3 for X86. DO NOT USE GCC 3.3 on
- the X86 platform. g++ 3.2, 3.4, and the 4.* series are known to work well.
-
--------------------------------------------------------------------------------
-
--------------------------------------------------------------------------------
-
- NOTES
-
--------------------------------------------------------------------------------
-
-To run the examples you may need to set PYTHONPATH. Note that the
-prefix and python version number in the path needs to match your
-installed version of python.
-
- $ export PYTHONPATH=/usr/local/lib/python2.7/dist-packages
-
-You may want to add this to your shell init file (~/.bash_profile if
-you use bash).
-
-
-Another handy trick if for example your fftw includes and libs are
-installed in, say ~/local/include and ~/local/lib, instead of
-/usr/local is this:
-
- $ export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$HOME/local/lib
- $ make CPPFLAGS="-I$HOME/local/include"
-
-
-Sometimes the prerequisites are installed in a location which is not
-included in the default compiler and linker search paths. This
-happens with pkgsrc and NetBSD. To build, tell configure to use these
-locations:
-
- LDFLAGS="-L/usr/pkg/lib -R/usr/pkg/lib" CPPFLAGS="-I/usr/pkg/include"
./configure --prefix=/usr/gnuradio
-
-
--------------------------------------------------------------------------------
-
- Legal Matters
-
--------------------------------------------------------------------------------
-
-Some files have been changed many times throughout the
-years. Copyright notices at the tops of these files list which years
-changes have been made. For some files, changes have occurred in many
-consecutive years. These files may often have the format of a year
-range (e.g., "2006 - 2011"), which indicates that these files have had
-copyrightable changes made during each year in the range, inclusive.
-
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..8d6fc4b
--- /dev/null
+++ b/README.md
@@ -0,0 +1,23 @@
+ESA Summer of Code in Space 2015 POLAR Code implementation for GNU Radio
+================
+
+This is a GNU Radio clone which contains all the files for my POLAR code
project for ESA Summer of Code in Space 2015. My name is Johannes Demel and I
am a Master degree student at *Karlsruhe Institute of Technology (KIT)*.
+
+Basics
+----------
+The project is based on Erdal Arikan's 2009 paper *Channel Polarization: A
Method for Constructing Capacity-Achieving Codes for Symmetric Binary-Input
Memoryless Channels*. An in-depth explanation of my project is given in my
project proposal <https://github.com/jdemel/socis-proposal>.
+
+Structure
+----------
+The aim of the project is to add POLAR code encoder, decoder, channel
construction capabilities to GNU Radio. POLAR codes are used for channel
coding, thus all files shall reside within the *gr-fec* submodule. Also encoder
and decoder shall use the FECAPI. The decision to do all development in-tree is
driven by previous experience that out-of-tree modules are often terribly hard
to merge later on and it may take years to do so.
+
+Project
+-------
+A quick overview of my project timeline goes as follows.
+
+1. create Python testcode
+2. implement encoder and decoder in C++ for FECAPI.
+3. optimize blocks and create VOLK kernels.
+4. implement channel construction algorithms.
+
+This may get continuously updated and I will add more details during project
development.
\ No newline at end of file
diff --git a/gr-fec/python/fec/polar/README.md
b/gr-fec/python/fec/polar/README.md
new file mode 100644
index 0000000..2bd00dc
--- /dev/null
+++ b/gr-fec/python/fec/polar/README.md
@@ -0,0 +1,4 @@
+POLAR Code Python test functions module
+===========
+
+This folder contains all the necessary files for POLAR code testcode. It shall
serve as a reference later on.
\ No newline at end of file
diff --git a/gr-fec/python/fec/polar/bit_reversed_polar_encoding_scheme.svg
b/gr-fec/python/fec/polar/bit_reversed_polar_encoding_scheme.svg
new file mode 100644
index 0000000..84849c9
--- /dev/null
+++ b/gr-fec/python/fec/polar/bit_reversed_polar_encoding_scheme.svg
@@ -0,0 +1,1160 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!-- Created with Inkscape (http://www.inkscape.org/) -->
+
+<svg
+ xmlns:dc="http://purl.org/dc/elements/1.1/";
+ xmlns:cc="http://creativecommons.org/ns#";
+ xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#";
+ xmlns:svg="http://www.w3.org/2000/svg";
+ xmlns="http://www.w3.org/2000/svg";
+ xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd";
+ xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape";
+ width="744.09448819"
+ height="1052.3622047"
+ id="svg2"
+ version="1.1"
+ inkscape:version="0.48.4 r9939"
+ sodipodi:docname="bit_reversed_polar_encoding_scheme.svg">
+ <defs
+ id="defs4" />
+ <sodipodi:namedview
+ id="base"
+ pagecolor="#ffffff"
+ bordercolor="#666666"
+ borderopacity="1.0"
+ inkscape:pageopacity="0.0"
+ inkscape:pageshadow="2"
+ inkscape:zoom="0.94"
+ inkscape:cx="186.70213"
+ inkscape:cy="520"
+ inkscape:document-units="px"
+ inkscape:current-layer="layer1"
+ showgrid="false"
+ inkscape:window-width="1920"
+ inkscape:window-height="1056"
+ inkscape:window-x="1920"
+ inkscape:window-y="24"
+ inkscape:window-maximized="1" />
+ <metadata
+ id="metadata7">
+ <rdf:RDF>
+ <cc:Work
+ rdf:about="">
+ <dc:format>image/svg+xml</dc:format>
+ <dc:type
+ rdf:resource="http://purl.org/dc/dcmitype/StillImage"; />
+ </cc:Work>
+ </rdf:RDF>
+ </metadata>
+ <g
+ inkscape:label="Layer 1"
+ inkscape:groupmode="layer"
+ id="layer1">
+ <flowRoot
+ xml:space="preserve"
+ id="flowRoot3755"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"><flowRegion
+ id="flowRegion3757"><rect
+ id="rect3759"
+ width="190.42554"
+ height="154.25533"
+ x="75.531914"
+ y="201.29836" /></flowRegion><flowPara
+ id="flowPara3761">000sfdgkzhlbx000000</flowPara></flowRoot> <g
+ id="g4231"
+ transform="translate(-11.702127,-12.23403)">
+ <g
+ transform="translate(-77.659574,-114.89362)"
+ id="g3767">
+ <rect
+ id="rect2985"
+ width="42.553192"
+ height="41.489372"
+ x="136.17021"
+ y="258.74515"
+ style="fill:none;stroke:#000000;stroke-opacity:1" />
+ <text
+ xml:space="preserve"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ x="144.68085"
+ y="293.85156"
+ id="text3763"
+ sodipodi:linespacing="125%"><tspan
+ sodipodi:role="line"
+ id="tspan3765"
+ x="144.68085"
+ y="293.85156">0</tspan></text>
+ </g>
+ </g>
+ <g
+ id="g4225"
+ transform="translate(-5.319148,-34.650486)">
+ <g
+ transform="translate(-84.042553,-23.936162)"
+ id="g3767-0">
+ <rect
+ id="rect2985-6"
+ width="42.553192"
+ height="41.489372"
+ x="136.17021"
+ y="258.74515"
+ style="fill:none;stroke:#000000;stroke-opacity:1" />
+ <text
+ xml:space="preserve"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ x="144.68085"
+ y="293.85156"
+ id="text3763-6"
+ sodipodi:linespacing="125%"><tspan
+ sodipodi:role="line"
+ id="tspan3765-8"
+ x="144.68085"
+ y="293.85156">4</tspan></text>
+ </g>
+ </g>
+ <g
+ id="g4219"
+ transform="translate(-10.638297,-24.62003)">
+ <g
+ transform="translate(-78.723404,34.574476)"
+ id="g3767-8">
+ <rect
+ id="rect2985-8"
+ width="42.553192"
+ height="41.489372"
+ x="136.17021"
+ y="258.74515"
+ style="fill:none;stroke:#000000;stroke-opacity:1" />
+ <text
+ xml:space="preserve"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ x="144.68085"
+ y="293.85156"
+ id="text3763-62"
+ sodipodi:linespacing="125%"><tspan
+ sodipodi:role="line"
+ id="tspan3765-7"
+ x="144.68085"
+ y="293.85156">2</tspan></text>
+ </g>
+ </g>
+ <g
+ id="g4213"
+ transform="translate(-7.446808,-9.270503)">
+ <g
+ transform="translate(-81.914893,87.765966)"
+ id="g3767-5">
+ <rect
+ id="rect2985-3"
+ width="42.553192"
+ height="41.489372"
+ x="136.17021"
+ y="258.74515"
+ style="fill:none;stroke:#000000;stroke-opacity:1" />
+ <text
+ xml:space="preserve"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ x="144.68085"
+ y="293.85156"
+ id="text3763-7"
+ sodipodi:linespacing="125%"><tspan
+ sodipodi:role="line"
+ id="tspan3765-81"
+ x="144.68085"
+ y="293.85156">6</tspan></text>
+ </g>
+ </g>
+ <g
+ id="g4207"
+ transform="translate(-5.319148,-12.00608)">
+ <g
+ transform="translate(-84.042553,159.04256)"
+ id="g3767-1">
+ <rect
+ id="rect2985-5"
+ width="42.553192"
+ height="41.489372"
+ x="136.17021"
+ y="258.74515"
+ style="fill:none;stroke:#000000;stroke-opacity:1" />
+ <text
+ xml:space="preserve"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ x="144.68085"
+ y="293.85156"
+ id="text3763-3"
+ sodipodi:linespacing="125%"><tspan
+ sodipodi:role="line"
+ id="tspan3765-87"
+ x="144.68085"
+ y="293.85156">1</tspan></text>
+ </g>
+ </g>
+ <g
+ id="g4201"
+ transform="translate(-11.702127,0.15198)">
+ <g
+ transform="translate(-77.659574,215.42554)"
+ id="g3767-54">
+ <rect
+ id="rect2985-4"
+ width="42.553192"
+ height="41.489372"
+ x="136.17021"
+ y="258.74515"
+ style="fill:none;stroke:#000000;stroke-opacity:1" />
+ <text
+ xml:space="preserve"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ x="144.68085"
+ y="293.85156"
+ id="text3763-4"
+ sodipodi:linespacing="125%"><tspan
+ sodipodi:role="line"
+ id="tspan3765-6"
+ x="144.68085"
+ y="293.85156">5</tspan></text>
+ </g>
+ </g>
+ <g
+ id="g4195"
+ transform="translate(-3.191489,-1.51975)">
+ <g
+ transform="translate(-86.170212,285.63831)"
+ id="g3767-9">
+ <rect
+ id="rect2985-7"
+ width="42.553192"
+ height="41.489372"
+ x="136.17021"
+ y="258.74515"
+ style="fill:none;stroke:#000000;stroke-opacity:1" />
+ <text
+ xml:space="preserve"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ x="144.68085"
+ y="293.85156"
+ id="text3763-8"
+ sodipodi:linespacing="125%"><tspan
+ sodipodi:role="line"
+ id="tspan3765-0"
+ x="144.68085"
+ y="293.85156">3</tspan></text>
+ </g>
+ </g>
+ <g
+ id="g4189"
+ transform="translate(0,-2.12766)">
+ <g
+ transform="translate(-89.361701,354.78724)"
+ id="g3767-6">
+ <rect
+ id="rect2985-73"
+ width="42.553192"
+ height="41.489372"
+ x="136.17021"
+ y="258.74515"
+ style="fill:none;stroke:#000000;stroke-opacity:1" />
+ <text
+ xml:space="preserve"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ x="144.68085"
+ y="293.85156"
+ id="text3763-64"
+ sodipodi:linespacing="125%"><tspan
+ sodipodi:role="line"
+ id="tspan3765-07"
+ x="144.68085"
+ y="293.85156">7</tspan></text>
+ </g>
+ </g>
+ <g
+ id="g3767-4"
+ transform="translate(188.08511,-127.12765)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-0" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-38"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-9"
+ sodipodi:role="line">0</tspan></text>
+ </g>
+ <g
+ id="g3767-01"
+ transform="translate(188.08511,-58.586648)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-68" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-40"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-69"
+ sodipodi:role="line">1</tspan></text>
+ </g>
+ <g
+ id="g3767-83"
+ transform="translate(188.08511,9.9544464)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-44" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-9"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-1"
+ sodipodi:role="line">2</tspan></text>
+ </g>
+ <g
+ id="g3767-010"
+ transform="translate(188.08511,78.495463)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-9" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-1"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-68"
+ sodipodi:role="line">3</tspan></text>
+ </g>
+ <g
+ id="g3767-7"
+ transform="translate(188.08511,147.03648)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-42" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-93"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-89"
+ sodipodi:role="line">4</tspan></text>
+ </g>
+ <g
+ id="g3767-87"
+ transform="translate(188.08511,215.57752)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-1" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-90"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-97"
+ sodipodi:role="line">5</tspan></text>
+ </g>
+ <g
+ id="g3767-79"
+ transform="translate(188.08511,284.11856)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-36" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-99"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-3"
+ sodipodi:role="line">6</tspan></text>
+ </g>
+ <g
+ id="g3767-40"
+ transform="translate(188.08511,352.65958)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-66" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-16"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-5"
+ sodipodi:role="line">7</tspan></text>
+ </g>
+ <g
+ id="g3767-544"
+ transform="translate(344.97872,-127.12765)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-2" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-88"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-32"
+ sodipodi:role="line">0</tspan></text>
+ </g>
+ <g
+ id="g3767-94"
+ transform="translate(344.97872,-58.586648)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-08" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-5"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-10"
+ sodipodi:role="line">1</tspan></text>
+ </g>
+ <g
+ id="g3767-59"
+ transform="translate(344.97872,9.9544464)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-76" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-46"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-65"
+ sodipodi:role="line">2</tspan></text>
+ </g>
+ <g
+ id="g3767-12"
+ transform="translate(344.97872,78.495463)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-67" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-87"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-56"
+ sodipodi:role="line">3</tspan></text>
+ </g>
+ <g
+ id="g3767-2"
+ transform="translate(344.97872,147.03648)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-92" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-47"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-09"
+ sodipodi:role="line">4</tspan></text>
+ </g>
+ <g
+ id="g3767-11"
+ transform="translate(344.97872,215.57752)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-449" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-95"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-06"
+ sodipodi:role="line">5</tspan></text>
+ </g>
+ <g
+ id="g3767-69"
+ transform="translate(344.97872,284.11856)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-31" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-51"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-683"
+ sodipodi:role="line">6</tspan></text>
+ </g>
+ <g
+ id="g3767-47"
+ transform="translate(344.97872,352.65958)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-32" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-2"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-94"
+ sodipodi:role="line">7</tspan></text>
+ </g>
+ <g
+ id="g3767-111"
+ transform="translate(502,-127.12765)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-14" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-0"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-4"
+ sodipodi:role="line">0</tspan></text>
+ </g>
+ <g
+ id="g3767-56"
+ transform="translate(502,-58.586648)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-85" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-66"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-74"
+ sodipodi:role="line">1</tspan></text>
+ </g>
+ <g
+ id="g3767-46"
+ transform="translate(502,9.9544464)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-75" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-48"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-42"
+ sodipodi:role="line">2</tspan></text>
+ </g>
+ <g
+ id="g3767-18"
+ transform="translate(502,78.495463)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-27" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-24"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-67"
+ sodipodi:role="line">3</tspan></text>
+ </g>
+ <g
+ id="g3767-80"
+ transform="translate(502,147.03648)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-89" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-405"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-17"
+ sodipodi:role="line">4</tspan></text>
+ </g>
+ <g
+ id="g3767-66"
+ transform="translate(502,215.57752)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-52" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-39"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-82"
+ sodipodi:role="line">5</tspan></text>
+ </g>
+ <g
+ id="g3767-946"
+ transform="translate(502,284.11856)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-70" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-01"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-02"
+ sodipodi:role="line">6</tspan></text>
+ </g>
+ <g
+ id="g3767-837"
+ transform="translate(502,352.65958)">
+ <rect
+ style="fill:none;stroke:#000000;stroke-opacity:1"
+ y="258.74515"
+ x="136.17021"
+ height="41.489372"
+ width="42.553192"
+ id="rect2985-60" />
+ <text
+ sodipodi:linespacing="125%"
+ id="text3763-76"
+ y="293.85156"
+ x="144.68085"
+
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
+ xml:space="preserve"><tspan
+ y="293.85156"
+ x="144.68085"
+ id="tspan3765-08"
+ sodipodi:role="line">7</tspan></text>
+ </g>
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 89.361703,152.36218 234.893617,0"
+ id="path4237"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g4231"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-4"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 366.80851,152.36218 114.34042,0"
+ id="path4239"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-4"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-544"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 523.70212,152.36218 114.46809,0"
+ id="path4241"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-544"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-111"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 89.361703,220.90319 234.893617,0"
+ id="path4243"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g4225"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-01"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 366.80851,220.90319 114.34042,0"
+ id="path4245"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-01"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-94"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 523.70212,220.90319 114.46809,0"
+ id="path4247"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-94"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-56"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 89.361703,289.44428 234.893617,0"
+ id="path4249"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g4219"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-83"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 366.80851,289.44428 114.34042,0"
+ id="path4251"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-83"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-59"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 523.70212,289.44428 114.46809,0"
+ id="path4253"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-59"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-46"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 89.361703,357.9853 234.893617,0"
+ id="path4255"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g4213"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-010"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 366.80851,357.9853 114.34042,0"
+ id="path4257"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-010"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-12"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 523.70212,357.9853 114.46809,0"
+ id="path4259"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-12"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-18"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 89.361703,426.52631 234.893617,0"
+ id="path4261"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g4207"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-7"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 366.80851,426.52631 114.34042,0"
+ id="path4263"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-7"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-2"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 523.70212,426.52631 114.46809,0"
+ id="path4265"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-2"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-80"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 89.361703,495.06735 234.893617,0"
+ id="path4267"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g4201"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-87"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 366.80851,495.06735 114.34042,0"
+ id="path4269"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-87"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-11"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 523.70212,495.06735 114.46809,0"
+ id="path4271"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-11"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-66"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 366.80851,563.60839 114.34042,0"
+ id="path4275"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-79"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-69"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 523.70212,563.60839 114.46809,0"
+ id="path4277"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-69"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-946"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 366.80851,632.14941 114.34042,0"
+ id="path4285"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-40"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-47"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 523.70212,632.14941 114.46809,0"
+ id="path4287"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-47"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-837"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 89.361703,563.60839 234.893617,0"
+ id="path4289"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g4195"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-79"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m 89.361703,632.14941 234.893617,0"
+ id="path4291"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g4189"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-40"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="M 523.70212,211.61579 638.17021,161.64958"
+ id="path4293"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-94"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-111"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="M 523.70212,348.6979 638.17021,298.73168"
+ id="path4295"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-12"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-46"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="M 523.70212,485.77995 638.17021,435.81372"
+ id="path4297"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-11"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-80"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="M 523.70212,622.86201 638.17021,572.89579"
+ id="path4299"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-47"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-946"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="M 366.80851,545.01847 481.14893,445.11624"
+ id="path4301"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-79"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-2"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="M 366.80851,613.55949 481.14893,513.65727"
+ id="path4303"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-40"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-11"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="M 366.80851,270.85436 481.14893,170.95211"
+ id="path4305"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-83"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-544"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="M 366.80851,339.39537 481.14893,239.49311"
+ id="path4307"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g3767-010"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-94"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="M 89.078178,405.78163 324.53885,173.10687"
+ id="path4309"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g4207"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-4"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="M 89.078175,474.32267 324.53885,241.64787"
+ id="path4311"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g4201"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-01"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="M 89.078179,542.86371 324.53885,310.18897"
+ id="path4313"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g4195"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-83"
+ inkscape:connection-end-point="d4" />
+ <path
+
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="M 89.078179,611.40473 324.53885,378.72998"
+ id="path4315"
+ inkscape:connector-type="polyline"
+ inkscape:connector-curvature="0"
+ inkscape:connection-start="#g4189"
+ inkscape:connection-start-point="d4"
+ inkscape:connection-end="#g3767-010"
+ inkscape:connection-end-point="d4" />
+ </g>
+</svg>
diff --git a/gr-fec/python/fec/polar/channel_construction_bec.py
b/gr-fec/python/fec/polar/channel_construction_bec.py
new file mode 100644
index 0000000..cd61f46
--- /dev/null
+++ b/gr-fec/python/fec/polar/channel_construction_bec.py
@@ -0,0 +1,110 @@
+#!/usr/bin/env python
+#
+# Copyright 2015 Free Software Foundation, Inc.
+#
+# 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 as np
+import matplotlib.pyplot as plt
+
+def bsc_channel(p):
+ '''
+ binary symmetric channel (BSC)
+ output alphabet Y = {0, 1} and
+ W(0|0) = W(1|1) and W(1|0) = W(0|1)
+
+ this function returns a prob's vector for a BSC
+ p denotes an erroneous transistion
+ '''
+ if not (p >= 0.0 and p <= 1.0):
+ print "given p is out of range!"
+ return np.array([], dtype=float)
+
+ # 0 -> 0, 0 -> 1, 1 -> 0, 1 -> 1
+ W = np.array([[1 - p, p], [p, 1 - p]], dtype=float)
+ return W
+
+
+def bec_channel(eta):
+ '''
+ binary erasure channel (BEC)
+ for each y e Y
+ W(y|0) * W(y|1) = 0 or W(y|0) = W(y|1)
+ transistions are 1 -> 1 or 0 -> 0 or {0, 1} -> ? (erased symbol)
+ '''
+
+ # looks like BSC but should be interpreted differently.
+ W = np.array((1 - eta, eta, 1 - eta), dtype=float)
+ return W
+
+
+def odd_rec(iwn):
+ return iwn ** 2
+
+
+def even_rec(iwn):
+ return 2 * iwn - iwn ** 2
+
+
+def calc_one_recursion(iw0):
+ iw1 = np.zeros(2 * len(iw0)) # double values
+ for i in range(len(iw0)):
+ # careful indices screw you because paper is '1' based :(
+ iw1[2 * i] = odd_rec(iw0[i])
+ iw1[2 * i + 1] = even_rec(iw0[i])
+ return iw1
+
+
+def calculate_bec_channel_capacities(eta, n):
+ iw = 1 - eta # holds for BEC as stated in paper
+ iw = np.array([iw, ], dtype=float)
+ lw = int(np.log2(n))
+ for i in range(lw):
+ iw = calc_one_recursion(iw)
+ return iw
+
+
+def get_frozen_bit_indices_from_capacities(chan_caps, nfrozen):
+ indexes = np.array([], dtype=int)
+ while indexes.size < nfrozen:
+ index = np.argmin(chan_caps)
+ indexes = np.append(indexes, index)
+ chan_caps[index] = 1.0
+ return np.sort(indexes)
+
+
+def bec_channel_contruction_tests():
+ n = 2 ** 10
+ k = n // 2
+ eta = 0.3
+ bec_capacities = calculate_bec_channel_capacities(eta, n)
+ print(bec_capacities)
+
+ frozenbits_position =
get_frozen_bit_indices_from_capacities(bec_capacities, k)
+
+ print('frozenbit_positions:')
+ print(frozenbits_position)
+ print('frozenbit_num:', frozenbits_position.size)
+
+
+def main():
+ print 'channel construction main'
+ bec_channel_contruction_tests()
+
+
+if __name__ == '__main__':
+ main()
diff --git a/gr-fec/python/fec/polar/channel_construction_bsc.py
b/gr-fec/python/fec/polar/channel_construction_bsc.py
new file mode 100755
index 0000000..8cee7ce
--- /dev/null
+++ b/gr-fec/python/fec/polar/channel_construction_bsc.py
@@ -0,0 +1,206 @@
+#!/usr/bin/env python
+#
+# Copyright 2015 Free Software Foundation, Inc.
+#
+# 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 as np
+from helper_functions import *
+import matplotlib.pyplot as plt
+
+
+# def bit_reverse(value, n):
+# # is this really missing in NumPy???
+# bits = np.zeros(n, type(value))
+# for index in range(n):
+# mask = 1
+# mask = np.left_shift(mask, index)
+# bit = np.bitwise_and(value, mask)
+# bit = np.right_shift(bit, index)
+# bits[index] = bit
+# bits = bits[::-1]
+# result = 0
+# for index, bit in enumerate(bits):
+# bit = np.left_shift(bit, index)
+# result += bit
+# return result
+
+
+def get_Bn(n):
+ # this is a bit reversal matrix.
+ lw = int(np.log2(n)) # number of used bits
+ indexes = [bit_reverse(i, lw) for i in range(n)]
+ Bn = np.zeros((n, n), type(n))
+ for i, index in enumerate(indexes):
+ Bn[i][index] = 1
+ return Bn
+
+
+def get_Fn(n):
+ # this matrix defines the actual channel combining.
+ if n == 1:
+ return np.array([1, ])
+ F2 = np.array([[1, 0], [1, 1]], np.int)
+ nump = int(np.log2(n)) - 1 # number of Kronecker products to calculate
+ Fn = F2
+ for i in range(nump):
+ Fn = np.kron(Fn, F2)
+ return Fn
+
+def get_Gn(n):
+ # this matrix is called generator matrix
+ if not is_power_of_two(n):
+ print "invalid input"
+ return None
+ if n == 1:
+ return np.array([1, ])
+ Bn = get_Bn(n)
+ Fn = get_Fn(n)
+ Gn = np.dot(Bn, Fn)
+ return Gn
+
+
+def mutual_information(w):
+ '''
+ calculate mutual information I(W)
+ I(W) = sum over y e Y ( sum over x e X ( ... ) )
+ .5 W(y|x) log frac { W(y|x) }{ .5 W(y|0) + .5 W(y|1) }
+ '''
+ ydim, xdim = np.shape(w)
+ i = 0.0
+ for y in range(ydim):
+ for x in range(xdim):
+ v = w[y][x] * np.log2(w[y][x] / (0.5 * w[y][0] + 0.5 * w[y][1]))
+ i += v
+ i /= 2.0
+ return i
+
+
+def bhattacharyya_parameter(w):
+ '''bhattacharyya parameter is a measure of similarity between two prob.
distributions'''
+ # sum over all y e Y for sqrt( W(y|0) * W(y|1) )
+ dim = np.shape(w)
+ ydim = dim[0]
+ xdim = dim[1]
+ z = 0.0
+ for y in range(ydim):
+ z += np.sqrt(w[y][0] * w[y][1])
+ # need all
+ return z
+
+
+def w_transition_prob(y, u, p):
+ return p[y == u]
+ # return 1 - p if y == u else p
+
+# @profile
+def calculate_joint_transition_probability(N, y, x, transition_probs):
+ single_ws = np.empty(N)
+ single_ws[y == x] = transition_probs[True]
+ single_ws[y != x] = transition_probs[False]
+ return np.prod(single_ws)
+
+
+# @profile
+def w_split_prob(y, u, G, transition_probs):
+ ''' Calculate channel splitting probabilities. '''
+ N = len(y) # number of combined channels
+ df = N - len(u) # degrees of freedom.
+ prob = 0.0
+ for uf in range(2 ** df):
+ utemp = unpack_byte(np.array([uf, ]), df)
+ ub = np.concatenate((u, utemp))
+ x = np.dot(ub, G) % 2
+ true_num = np.sum(y == x)
+ false_num = N - true_num
+ w = (transition_probs[True] ** true_num) * (transition_probs[False] **
false_num)
+ # w = calculate_joint_transition_probability(N, y, x, transition_probs)
+ prob += w
+ divider = 1.0 / (2 ** (N - 1))
+ return divider * prob
+
+# @profile
+def wn_split_channel(N, p):
+ G = get_Gn(N)
+ y = np.zeros((N, ), dtype=np.uint8)
+ n = 1
+ u = np.zeros((n + 1, ), dtype=np.uint8)
+ transition_probs = np.array([p, 1 - p], dtype=float)
+
+ z_params = []
+ c_params = []
+ for w in range(N):
+ nentries = (2 ** N) * (2 ** w)
+ print "for w=", w, " nentries=", nentries
+ w_probs = np.zeros((nentries, 2), dtype=float)
+ for y in range(2 ** N):
+ yb = unpack_byte(np.array([y, ]), N)
+ for u in range(2 ** (w + 1)):
+ ub = unpack_byte(np.array([u, ]), w + 1)
+ wp = w_split_prob(yb, ub, G, transition_probs)
+ ufixed = ub[0:-1]
+ yindex = y * (2 ** w) + pack_byte(ufixed)
+ uindex = ub[-1]
+ w_probs[yindex][uindex] = wp
+
+ z = bhattacharyya_parameter(w_probs)
+ z_params.append(z)
+ c = mutual_information(w_probs)
+ c_params.append(c)
+ print "W=", w, "Z=", z, "capacity=", c
+
+ return z_params, c_params
+
+
+def calculate_z_param(x):
+ # variables etc taken from paper bei Ido Tal et al.
+ # name there is f(x)
+ # x is the cross over probability of a BSC.
+ return 2 * np.sqrt(x * (1 - x))
+
+
+def calculate_capacity(x):
+ # in paper it is called g(x)
+ return -1. * x * np.log(x) - (1 - x) * np.log(1 - x)
+
+
+def splitting_masses_algorithm(n, k):
+ m = 2 ** n
+ p0 = 1.0 / m
+ mass_vec
+
+
+def main():
+ print 'channel construction BSC main'
+ n = 3
+ m = 2 ** n
+ k = m // 2
+ eta = 0.3
+ p = 0.1
+ # z, c = wn_split_channel(m, p)
+ # print(z)
+ # print(c)
+
+ u = np.zeros(m, dtype=bool)
+ G = get_Gn(m).astype(dtype=bool)
+ print G
+ print np.dot(u, G)
+
+
+
+if __name__ == '__main__':
+ main()
diff --git a/gr-fec/python/fec/polar/common.py
b/gr-fec/python/fec/polar/common.py
new file mode 100644
index 0000000..b4b152d
--- /dev/null
+++ b/gr-fec/python/fec/polar/common.py
@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+#
+# Copyright 2015 Free Software Foundation, Inc.
+#
+# 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 as np
+from helper_functions import *
+
+'''
+PolarCommon holds value checks and common initializer code for both Encoder
and Decoder.
+'''
+
+
+class PolarCommon:
+ def __init__(self, n, k, frozen_bit_position, frozenbits=None):
+ if not is_power_of_two(n):
+ raise ValueError("n={0} is not a power of 2!".format(n))
+ if frozenbits is None:
+ frozenbits = np.zeros(n - k, dtype=np.int)
+ if not len(frozenbits) == n - k:
+ raise ValueError("len(frozenbits)={0} is not equal to
n-k={1}!".format(len(frozenbits), n - k))
+ if not frozenbits.dtype == np.int:
+ frozenbits = frozenbits.astype(dtype=int)
+ if not len(frozen_bit_position) == (n - k):
+ raise ValueError("len(frozen_bit_position)={0} is not equal to
n-k={1}!".format(len(frozen_bit_position), n - k))
+ if not frozen_bit_position.dtype == np.int:
+ frozen_bit_position = frozen_bit_position.astype(dtype=int)
+
+ self.bit_reverse_positions = self._vector_bit_reversed(np.arange(n,
dtype=int), int(np.log2(n)))
+ self.N = n
+ self.power = int(np.log2(self.N))
+ self.K = k
+ self.frozenbits = frozenbits
+ self.frozen_bit_position = frozen_bit_position
+ self.info_bit_position = np.delete(np.arange(self.N),
self.frozen_bit_position)
+
+ def _insert_frozen_bits(self, u):
+ prototype = np.empty(self.N, dtype=int)
+ prototype[self.frozen_bit_position] = self.frozenbits
+ prototype[self.info_bit_position] = u
+ return prototype
+
+ def _extract_info_bits(self, y):
+ return y[self.info_bit_position]
+
+ def _reverse_bits(self, vec):
+ return vec[self.bit_reverse_positions]
+
+ def _vector_bit_reversed(self, vec, n):
+ return bit_reverse_vector(vec, n)
+
+ def info_print(self):
+ print "POLAR code ({0}, {1})".format(self.N, self.K)
diff --git a/gr-fec/python/fec/polar/decoder.py
b/gr-fec/python/fec/polar/decoder.py
new file mode 100644
index 0000000..d74f1f9
--- /dev/null
+++ b/gr-fec/python/fec/polar/decoder.py
@@ -0,0 +1,278 @@
+#!/usr/bin/env python
+#
+# Copyright 2015 Free Software Foundation, Inc.
+#
+# 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 as np
+from common import PolarCommon
+
+# for dev
+from encoder import PolarEncoder
+from matplotlib import pyplot as plt
+
+
+class PolarDecoder(PolarCommon):
+ def __init__(self, n, k, frozen_bit_position, frozenbits=None):
+ PolarCommon.__init__(self, n, k, frozen_bit_position, frozenbits)
+
+ self.error_probability = 0.1 # this is kind of a dummy value. usually
chosen individually.
+ self.bsc_lr = ((1 - self.error_probability) / self.error_probability,
self.error_probability / (1 - self.error_probability))
+ self.bsc_llrs = np.log(self.bsc_lr)
+
+ def _llr_bit(self, bit):
+ return self.bsc_llrs[bit]
+
+ def _llr_odd(self, la, lb):
+ # this functions uses the min-sum approximation
+ # exact formula: np.log((np.exp(la + lb) + 1) / (np.exp(la) +
np.exp(lb)))
+ return np.sign(la) * np.sign(lb) * np.minimum(np.abs(la), np.abs(lb))
+
+ _f_vals = np.array((1.0, -1.0), dtype=float)
+
+ def _llr_even(self, la, lb, f):
+ return (la * self._f_vals[f]) + lb
+
+ def _llr_bit_decision(self, llr):
+ if llr < 0.0:
+ ui = int(1)
+ else:
+ ui = int(0)
+ return ui
+
+ def _retrieve_bit_from_llr(self, lr, pos):
+ f_index = np.where(self.frozen_bit_position == pos)[0]
+ if not f_index.size == 0:
+ ui = self.frozenbits[f_index][0]
+ else:
+ ui = self._llr_bit_decision(lr)
+ return ui
+
+ def _lr_bit(self, bit):
+ return self.bsc_lr[bit]
+
+ def _lr_odd(self, la, lb):
+ # la is upper branch and lb is lower branch
+ return (la * lb + 1) / (la + lb)
+
+ def _lr_even(self, la, lb, f):
+ # la is upper branch and lb is lower branch, f is last decoded bit.
+ return (la ** (1 - (2 * f))) * lb
+
+ def _lr_bit_decision(self, lr):
+ if lr < 1:
+ return int(1)
+ return int(0)
+
+ def _get_even_indices_values(self, u_hat):
+ # looks like overkill for some indexing, but zero and one based
indexing mix-up gives you haedaches.
+ return u_hat[1::2]
+
+ def _get_odd_indices_values(self, u_hat):
+ return u_hat[0::2]
+
+ def _calculate_lrs(self, y, u):
+ ue = self._get_even_indices_values(u)
+ uo = self._get_odd_indices_values(u)
+ ya = y[0:y.size//2]
+ yb = y[(y.size//2):]
+ la = self._lr_decision_element(ya, (ue + uo) % 2)
+ lb = self._lr_decision_element(yb, ue)
+ return la, lb
+
+ def _lr_decision_element(self, y, u):
+ if y.size == 1:
+ return self._llr_bit(y[0])
+ if u.size % 2 == 0: # use odd branch formula
+ la, lb = self._calculate_lrs(y, u)
+ return self._llr_odd(la, lb)
+ else:
+ ui = u[-1]
+ la, lb = self._calculate_lrs(y, u[0:-1])
+ return self._llr_even(la, lb, ui)
+
+ def _retrieve_bit_from_lr(self, lr, pos):
+ f_index = np.where(self.frozen_bit_position == pos)[0]
+ if not f_index.size == 0:
+ ui = self.frozenbits[f_index][0]
+ else:
+ ui = self._lr_bit_decision(lr)
+ return ui
+
+ def _lr_sc_decoder(self, y):
+ # this is the standard SC decoder as derived from the formulas. It
sticks to natural bit order.
+ u = np.array([], dtype=int)
+ for i in range(y.size):
+ lr = self._lr_decision_element(y, u)
+ ui = self._retrieve_bit_from_llr(lr, i)
+ u = np.append(u, ui)
+ return u
+
+ def _butterfly_decode_bits(self, pos, graph, u):
+ llr = graph[pos][0]
+ ui = self._llr_bit_decision(llr)
+ u = np.append(u, ui)
+ lower_right = pos + (self.N // 2)
+ la = graph[pos][1]
+ lb = graph[lower_right][1]
+ graph[lower_right][0] = self._llr_even(la, lb, ui)
+ llr = graph[lower_right][0]
+ ui = self._llr_bit_decision(llr)
+ u = np.append(u, ui)
+ return graph, u
+
+ def _lr_sc_decoder_efficient(self, y):
+ graph = np.full((self.N, self.power + 1), np.NaN, dtype=float)
+ for i in range(self.N):
+ graph[i][self.power] = self._llr_bit(y[i])
+ decode_order = self._vector_bit_reversed(np.arange(self.N), self.power)
+ decode_order = np.delete(decode_order, np.where(decode_order >= self.N
// 2))
+ u = np.array([], dtype=int)
+ for pos in decode_order:
+ graph = self._butterfly(pos, 0, graph, u)
+ graph, u = self._butterfly_decode_bits(pos, graph, u)
+ return u
+
+ def _stop_propagation(self, bf_entry_row, stage):
+ # calculate break condition
+ modulus = 2 ** (self.power - stage)
+ # stage_size = self.N // (2 ** stage)
+ # half_stage_size = stage_size // 2
+ half_stage_size = self.N // (2 ** (stage + 1))
+ stage_pos = bf_entry_row % modulus
+ return stage_pos >= half_stage_size
+
+ def _butterfly(self, bf_entry_row, stage, graph, u):
+ if not self.power > stage:
+ return graph
+
+ if self._stop_propagation(bf_entry_row, stage):
+ upper_right = bf_entry_row - self.N // (2 ** (stage + 1))
+ la = graph[upper_right][stage + 1]
+ lb = graph[bf_entry_row][stage + 1]
+ ui = u[-1]
+ graph[bf_entry_row][stage] = self._llr_even(la, lb, ui)
+ return graph
+
+ # activate right side butterflies
+ u_even = self._get_even_indices_values(u)
+ u_odd = self._get_odd_indices_values(u)
+ graph = self._butterfly(bf_entry_row, stage + 1, graph, (u_even +
u_odd) % 2)
+ lower_right = bf_entry_row + self.N // (2 ** (stage + 1))
+ graph = self._butterfly(lower_right, stage + 1, graph, u_even)
+
+ la = graph[bf_entry_row][stage + 1]
+ lb = graph[lower_right][stage + 1]
+ graph[bf_entry_row][stage] = self._llr_odd(la, lb)
+ return graph
+
+ def decode(self, data, is_packed=False):
+ if not len(data) == self.N:
+ raise ValueError("len(data)={0} is not equal to
n={1}!".format(len(data), self.N))
+ if is_packed:
+ data = np.unpackbits(data)
+ data = self._lr_sc_decoder_efficient(data)
+ data = self._extract_info_bits(data)
+ if is_packed:
+ data = np.packbits(data)
+ return data
+
+
+def test_reverse_enc_dec():
+ n = 16
+ k = 8
+ frozenbits = np.zeros(n - k)
+ frozenbitposition = np.array((0, 1, 2, 3, 4, 5, 8, 9), dtype=int)
+ bits = np.random.randint(2, size=k)
+ encoder = PolarEncoder(n, k, frozenbitposition, frozenbits)
+ decoder = PolarDecoder(n, k, frozenbitposition, frozenbits)
+ encoded = encoder.encode(bits)
+ print 'encoded:', encoded
+ rx = decoder.decode(encoded)
+ print 'bits:', bits
+ print 'rx :', rx
+ print (bits == rx).all()
+
+
+def compare_decoder_impls():
+ print '\nthis is decoder test'
+ n = 8
+ k = 4
+ frozenbits = np.zeros(n - k)
+ # frozenbitposition = np.array((0, 1, 2, 3, 4, 5, 8, 9), dtype=int)
+ frozenbitposition = np.array((0, 1, 2, 4), dtype=int)
+ # bits = np.ones(k, dtype=int)
+ bits = np.random.randint(2, size=k)
+ # bits = np.array([0, 1, 1, 1])
+ # bits = np.array([0, 1, 1, 0])
+ # bits = np.array([1, 0, 1, 0])
+ print 'bits:', bits
+ encoder = PolarEncoder(n, k, frozenbitposition, frozenbits)
+ decoder = PolarDecoder(n, k, frozenbitposition, frozenbits)
+ encoded = encoder.encode(bits)
+ print 'encoded:', encoded
+ rx_st = decoder._lr_sc_decoder(encoded)
+ rx_eff = decoder._lr_sc_decoder_efficient(encoded)
+ print 'standard :', rx_st
+ print 'efficient:', rx_eff
+ print (rx_st == rx_eff).all()
+
+
+
+def main():
+ power = 3
+ n = 2 ** power
+ k = 4
+ frozenbits = np.zeros(n - k, dtype=int)
+ frozenbitposition = np.array((0, 1, 2, 4), dtype=int)
+ frozenbitposition4 = np.array((0, 1), dtype=int)
+
+
+ encoder = PolarEncoder(n, k, frozenbitposition, frozenbits)
+ decoder = PolarDecoder(n, k, frozenbitposition, frozenbits)
+
+ bits = np.ones(k, dtype=int)
+ # bits = np.array([1, 0, 1, 0], dtype=int)
+ print "bits: ", bits
+ evec = encoder.encode(bits)
+ print "froz: ", encoder._insert_frozen_bits(bits)
+ print "evec: ", evec
+ # dvec = decoder.decode(evec)
+ # print "dec: ", dvec
+
+ # llr = decoder._llr(4, evec, np.array([0, 0, 0]))
+ # print "llr=", llr
+ evec[1] = 0
+ deced = decoder._lr_sc_decoder(evec)
+ print 'SC decoded:', deced
+
+
+
+
+ # test_reverse_enc_dec()
+ compare_decoder_impls()
+
+ # graph_decode()
+
+
+
+
+
+
+
+if __name__ == '__main__':
+ main()
\ No newline at end of file
diff --git a/gr-fec/python/fec/polar/encoder.py
b/gr-fec/python/fec/polar/encoder.py
new file mode 100644
index 0000000..6f87a22
--- /dev/null
+++ b/gr-fec/python/fec/polar/encoder.py
@@ -0,0 +1,135 @@
+#!/usr/bin/env python
+#
+# Copyright 2015 Free Software Foundation, Inc.
+#
+# 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 as np
+from common import PolarCommon
+
+
+class PolarEncoder(PolarCommon):
+ def __init__(self, n, k, frozen_bit_position, frozenbits=None):
+ PolarCommon.__init__(self, n, k, frozen_bit_position, frozenbits)
+ self.G = self._gn(n)
+
+ def _gn(self, n):
+ # this matrix is called generator matrix
+ if n == 1:
+ return np.array([1, ])
+ f = self._fn(n)
+ return f
+
+ def _fn(self, n):
+ # this matrix defines the actual channel combining.
+ if n == 1:
+ return np.array([1, ])
+ f2 = np.array([[1, 0], [1, 1]], np.int)
+ nump = int(np.log2(n)) - 1 # number of Kronecker products to calculate
+ fn = f2
+ for i in range(nump):
+ fn = np.kron(fn, f2)
+ return fn
+
+ def get_gn(self):
+ return self.G
+
+ def _prepare_input_data(self, vec):
+ vec = self._insert_frozen_bits(vec)
+ vec = self._reverse_bits(vec)
+ return vec
+
+ def _encode_matrix(self, data):
+ data = np.dot(data, self.G) % 2
+ data = data.astype(dtype=int)
+ return data
+
+ def _encode_efficient(self, vec):
+ nstages = int(np.log2(self.N))
+ pos = np.arange(self.N, dtype=int)
+ for i in range(nstages):
+ splitted = np.reshape(pos, (2 ** (i + 1), -1))
+ upper_branch = splitted[0::2].flatten()
+ lower_branch = splitted[1::2].flatten()
+ vec[upper_branch] = (vec[upper_branch] + vec[lower_branch]) % 2
+ return vec
+
+ def encode(self, data, is_packed=False):
+ if not len(data) == self.K:
+ raise ValueError("len(data)={0} is not equal to
k={1}!".format(len(data), self.K))
+ if is_packed:
+ data = np.unpackbits(data)
+ if np.max(data) > 1 or np.min(data) < 0:
+ raise ValueError("can only encode bits!")
+ data = self._prepare_input_data(data)
+ data = self._encode_efficient(data)
+ if is_packed:
+ data = np.packbits(data)
+ return data
+
+
+def compare_results(encoder, ntests, k):
+ for n in range(ntests):
+ bits = np.random.randint(2, size=k)
+ preped = encoder._prepare_input_data(bits)
+ menc = encoder._encode_matrix(preped)
+ fenc = encoder._encode_efficient(preped)
+ if (menc == fenc).all() == False:
+ return False
+ return True
+
+
+def test_pseudo_rate_1_encoder(encoder, ntests, k):
+ for n in range(ntests):
+ bits = np.random.randint(2, size=k)
+ u = encoder._prepare_input_data(bits)
+ fenc = encoder._encode_efficient(u)
+ u_hat = encoder._encode_efficient(fenc)
+ if not (u_hat == u).all():
+ print('rate-1 encoder/decoder failed')
+ print u
+ print u_hat
+ return False
+ return True
+
+
+def test_encoder_impls():
+ print('comparing encoder implementations, matrix vs. efficient')
+ ntests = 1000
+ n = 16
+ k = 8
+ frozenbits = np.zeros(n - k)
+ # frozenbitposition8 = np.array((0, 1, 2, 4), dtype=int) # keep it!
+ frozenbitposition = np.array((0, 1, 2, 3, 4, 5, 8, 9), dtype=int)
+ encoder = PolarEncoder(n, k, frozenbitposition, frozenbits)
+ print 'result:', compare_results(encoder, ntests, k)
+
+ print('Test rate-1 encoder/decoder chain results')
+ r1_test = test_pseudo_rate_1_encoder(encoder, ntests, k)
+ print 'test rate-1 encoder/decoder:', r1_test
+
+
+
+def main():
+ print "main in encoder"
+ test_encoder_impls()
+
+
+
+
+if __name__ == '__main__':
+ main()
\ No newline at end of file
diff --git a/gr-fec/python/fec/polar/helper_functions.py
b/gr-fec/python/fec/polar/helper_functions.py
new file mode 100644
index 0000000..1d82457
--- /dev/null
+++ b/gr-fec/python/fec/polar/helper_functions.py
@@ -0,0 +1,84 @@
+#!/usr/bin/env python
+#
+# Copyright 2015 Free Software Foundation, Inc.
+#
+# 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 as np
+
+
+def is_power_of_two(num):
+ if type(num) != int:
+ return False # make sure we only compute integers.
+ return num != 0 and ((num & (num - 1)) == 0)
+
+
+def bit_reverse(value, n):
+ # is this really missing in NumPy???
+ seq = np.int(value)
+ rev = np.int(0)
+ rmask = np.int(1)
+ lmask = np.int(2 ** (n - 1))
+ for i in range(n // 2):
+ shiftval = n - 1 - (i * 2)
+ rshift = np.left_shift(np.bitwise_and(seq, rmask), shiftval)
+ lshift = np.right_shift(np.bitwise_and(seq, lmask), shiftval)
+ rev = np.bitwise_or(rev, rshift)
+ rev = np.bitwise_or(rev, lshift)
+ rmask = np.left_shift(rmask, 1)
+ lmask = np.right_shift(lmask, 1)
+ if not n % 2 == 0:
+ rev = np.bitwise_or(rev, np.bitwise_and(seq, rmask))
+ return rev
+
+
+def bit_reverse_vector(vec, n):
+ return np.array([bit_reverse(e, n) for e in vec], dtype=vec.dtype)
+
+
+def unpack_byte(byte, nactive):
+ if np.amin(byte) < 0 or np.amax(byte) > 255:
+ return None
+ if not byte.dtype == np.uint8:
+ byte = byte.astype(np.uint8)
+ if nactive == 0:
+ return np.array([], dtype=np.uint8)
+ return np.unpackbits(byte)[-nactive:]
+
+
+def pack_byte(bits):
+ if len(bits) == 0:
+ return 0
+ if np.amin(bits) < 0 or np.amax(bits) > 1: # only '1' and '0' in bits
array allowed!
+ return None
+ bits = np.concatenate((np.zeros(8 - len(bits), dtype=np.uint8), bits))
+ res = np.packbits(bits)[0]
+ return res
+
+
+def main():
+ print 'helper functions'
+
+ for i in range(8):
+ print(i, 'is power of 2: ', is_power_of_two(i))
+ n = 2 ** 6
+ k = n // 2
+ eta = 0.3
+
+
+if __name__ == '__main__':
+ main()
diff --git a/gr-fec/python/fec/polar/polar_code.py
b/gr-fec/python/fec/polar/polar_code.py
new file mode 100644
index 0000000..f263c6d
--- /dev/null
+++ b/gr-fec/python/fec/polar/polar_code.py
@@ -0,0 +1,493 @@
+#!/usr/bin/env python
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from encoder import PolarEncoder as enc
+
+# input alphabet X is always {0,1}
+# output alphabet is arbitrary
+# transition probabilities are arbitrary W(y|x)
+
+
+def bit_reverse(value, n):
+ # is this really missing in NumPy???
+ bits = np.zeros(n, type(value))
+ for index in range(n):
+ mask = 1
+ mask = np.left_shift(mask, index)
+ bit = np.bitwise_and(value, mask)
+ bit = np.right_shift(bit, index)
+ bits[index] = bit
+ bits = bits[::-1]
+ result = 0
+ for index, bit in enumerate(bits):
+ bit = np.left_shift(bit, index)
+ result += bit
+ return result
+
+
+def get_Bn(n):
+ # this is a bit reversal matrix.
+ lw = int(np.log2(n)) # number of used bits
+ indexes = [bit_reverse(i, lw) for i in range(n)]
+ Bn = np.zeros((n, n), type(n))
+ for i, index in enumerate(indexes):
+ Bn[i][index] = 1
+ return Bn
+
+
+def get_Fn(n):
+ # this matrix defines the actual channel combining.
+ if n == 1:
+ return np.array([1, ])
+ F2 = np.array([[1, 0], [1, 1]], np.int)
+ nump = int(np.log2(n)) - 1 # number of Kronecker products to calculate
+ Fn = F2
+ for i in range(nump):
+ Fn = np.kron(Fn, F2)
+ return Fn
+
+def get_Gn(n):
+ # this matrix is called generator matrix
+ if not is_power_of_2(n):
+ print "invalid input"
+ return None
+ if n == 1:
+ return np.array([1, ])
+ Bn = get_Bn(n)
+ Fn = get_Fn(n)
+ Gn = np.dot(Bn, Fn)
+ return Gn
+
+
+def odd_rec(iwn):
+ return iwn ** 2
+
+
+def even_rec(iwn):
+ return 2 * iwn - iwn ** 2
+
+
+def calc_one_recursion(iw0):
+ iw1 = np.zeros(2 * len(iw0)) # double values
+ for i in range(len(iw0)):
+ # careful indices screw you because paper is '1' based :(
+ iw1[2 * i] = odd_rec(iw0[i])
+ iw1[2 * i + 1] = even_rec(iw0[i])
+ return iw1
+
+
+def calculate_bec_channel_capacities(eta, n):
+ iw = 1 - eta # holds for BEC as stated in paper
+ lw = int(np.log2(n))
+ iw = [iw, ]
+ for i in range(lw):
+ iw = calc_one_recursion(iw)
+ return iw
+
+
+def bsc_channel(p):
+ '''
+ binary symmetric channel (BSC)
+ output alphabet Y = {0, 1} and
+ W(0|0) = W(1|1) and W(1|0) = W(0|1)
+
+ this function returns a prob's vector for a BSC
+ p denotes an erroneous transistion
+ '''
+ if not (p >= 0.0 and p <= 1.0):
+ print "given p is out of range!"
+ return ()
+
+ # 0 -> 0, 0 -> 1, 1 -> 0, 1 -> 1
+ W = np.array([[1 - p, p], [p, 1 - p]])
+ return W
+
+
+def bec_channel(eta):
+ '''
+ binary erasure channel (BEC)
+ for each y e Y
+ W(y|0) * W(y|1) = 0 or W(y|0) = W(y|1)
+ transistions are 1 -> 1 or 0 -> 0 or {0, 1} -> ? (erased symbol)
+ '''
+ print eta
+
+ # looks like BSC but should be interpreted differently.
+ W = (1 - eta, eta, 1 - eta)
+ return W
+
+
+def is_power_of_2(num):
+ if type(num) != int:
+ return False # make sure we only compute integers.
+ return num != 0 and ((num & (num - 1)) == 0)
+
+
+def combine_W2(u):
+ # This function applies the channel combination for W2
+ x1 = (u[0] + u[1]) % 2
+ x2 = (u[1])
+ return np.array([x1, x2], np.int)
+
+
+def combine_W4(u):
+ im = np.concatenate((combine_W2(u[0:2]), combine_W2(u[2:4])))
+ print "combine_W4.im = ", im
+ swappy = im[1]
+ im[1] = im[2]
+ im[2] = swappy
+ print "combine_W4.reverse_shuffled = ", im
+
+ return np.concatenate((combine_W2(im[0:2]), combine_W2(im[2:4])))
+
+
+def combine_Wn(u):
+ '''Combine vector u for encoding. It's described in the "Channel
combining" section'''
+ # will throw error if len(u) isn't a power of 2!
+ n = len(u)
+ G = get_Gn(n)
+ return np.dot(u, G) % 2
+
+
+def unpack_byte(byte, nactive):
+ if np.amin(byte) < 0 or np.amax(byte) > 255:
+ return None
+ if not byte.dtype == np.uint8:
+ byte = byte.astype(np.uint8)
+ if nactive == 0:
+ return np.array([], dtype=np.uint8)
+ return np.unpackbits(byte)[-nactive:]
+
+
+def pack_byte(bits):
+ if len(bits) == 0:
+ return 0
+ if np.amin(bits) < 0 or np.amax(bits) > 1: # only '1' and '0' in bits
array allowed!
+ return None
+ bits = np.concatenate((np.zeros(8 - len(bits), dtype=np.uint8), bits))
+ res = np.packbits(bits)[0]
+ return res
+
+
+def calculate_conditional_prob(y, u, channel):
+ # y and u same size and 1d!
+ # channel 2 x 2 matrix!
+ x1 = (u[0] + u[1]) % 2
+ s = 0 if y[0] == x1 else 1
+ # print "W(", y[0], "|", u[0], "+", u[1], "=", s, ") =", channel[y[0]][x1]
+ w112 = channel[y[0]][x1]
+ w22 = channel[y[1]][u[1]]
+ return w112 * w22
+
+
+def calculate_w2_probs():
+ w0 = np.array([[0.9, 0.1], [0.1, 0.9]])
+ print w0
+
+ w2 = np.zeros((4, 4), dtype=float)
+
+ print "W(y1|u1+u2)"
+ for y in range(4):
+ ybits = unpack_byte(np.array([y, ]), 2)
+ for u in range(4):
+ ubits = unpack_byte(np.array([u, ]), 2)
+ prob = calculate_conditional_prob(ybits, ubits, w0)
+ w2[y][u] = prob
+ # print "W(y1,y2|u1,u2) =", prob
+ return w2
+
+
+def get_prob(y, u, channel):
+ return channel[y][u]
+
+
+def get_wn_prob(y, u, channel):
+ x = combine_Wn(u)
+ result = 1
+ for i in range(len(x)):
+ result *= get_prob(y[i], x[i], channel)
+ return result
+
+
+def calculate_Wn_probs(n, channel):
+ # print "calculate_wn_probs"
+ ncomb = 2 ** n
+ wn = np.ones((ncomb, ncomb), dtype=float)
+ for y in range(ncomb):
+ ybits = unpack_byte(np.array([y, ]), n)
+ for u in range(ncomb):
+ ubits = unpack_byte(np.array([u, ]), n)
+ xbits = combine_Wn(ubits)
+ wn[y][u] = get_wn_prob(ybits, ubits, channel)
+ return wn
+
+
+def calculate_channel_splitting_probs(wn, n):
+ print wn
+
+ wi = np.zeros((n, 2 ** n, 2 ** n), dtype=float)
+ divider = (1.0 / (2 ** (n - 1)))
+ for i in range(n):
+ for y in range(2 ** n):
+ ybits = unpack_byte(np.array([y, ]), n)
+ print
+ for u in range(2 ** n):
+ ubits = unpack_byte(np.array([u, ]), n)
+ usc = ubits[0:i]
+ u = ubits[i]
+ usum = ubits[i+1:]
+ fixed_pu = np.append(usc, u)
+
+ my_iter = []
+ if len(usum) == 0:
+ my_iter.append(np.append(np.zeros(8 - len(fixed_pu),
dtype=np.uint8), fixed_pu))
+ else:
+ uiter = np.arange(2 ** len(usum), dtype=np.uint8)
+ for ui in uiter:
+ element = unpack_byte(ui, len(usum))
+ item = np.append(fixed_pu, element)
+ item = np.append(np.zeros(8 - len(item),
dtype=np.uint8), item)
+ my_iter.append(item)
+ my_iter = np.array(my_iter)
+
+ prob_sum = 0
+ for item in my_iter:
+ upos = np.packbits(item)
+ # print "y=", np.binary_repr(y, n), "item=", item, "u=",
np.binary_repr(upos, n)
+ wni = wn[y][upos]
+ prob_sum += wni
+ prob_sum *= divider
+
+ print "W[{5}]({0},{1},{2}|{3}) = {4}".format(ybits[0],
ybits[1], usc, u, prob_sum, i)
+
+ # print "i=", i, "y=", np.binary_repr(y, n), " fixed=",
fixed_pu, "usum=", usum, " WN(i) =", prob_sum
+ wi[i][y][u] = prob_sum
+
+ for i in range(n):
+ print
+ for y in range(2 ** n):
+ ybits = unpack_byte(np.array([y, ]), n)
+ for u in range(2 ** n):
+
+ print "W[{}] ({},{}|{})".format(i, ybits[0], ybits[1], u)
+
+
+
+ return wi
+
+
+
+def mutual_information(w):
+ '''
+ calculate mutual information I(W)
+ I(W) = sum over y e Y ( sum over x e X ( ... ) )
+ .5 W(y|x) log frac { W(y|x) }{ .5 W(y|0) + .5 W(y|1) }
+ '''
+ ydim, xdim = np.shape(w)
+ i = 0.0
+ for y in range(ydim):
+ for x in range(xdim):
+ v = w[y][x] * np.log2(w[y][x] / (0.5 * w[y][0] + 0.5 * w[y][1]))
+ i += v
+ i /= 2.0
+ return i
+
+
+def capacity(w):
+ '''
+ find supremum I(W) of a channel.
+ '''
+ return w
+
+
+def bhattacharyya_parameter(w):
+ '''bhattacharyya parameter is a measure of similarity between two prob.
distributions'''
+ # sum over all y e Y for sqrt( W(y|0) * W(y|1) )
+ dim = np.shape(w)
+ ydim = dim[0]
+ xdim = dim[1]
+ z = 0.0
+ for y in range(ydim):
+ z += np.sqrt(w[y][0] * w[y][1])
+ # need all
+ return z
+
+
+def w_transition_prob(y, u, p):
+ return 1 - p if y == u else p
+
+
+def w_split_prob(y, u, G, p):
+ ''' Calculate channel splitting probabilities. '''
+ N = len(y) # number of combined channels
+ df = N - len(u) # degrees of freedom.
+ prob = 0.0
+ for uf in range(2 ** df):
+ utemp = unpack_byte(np.array([uf, ]), df)
+ ub = np.concatenate((u, utemp))
+ # print "y=", y, "\tu=", ub
+ x = np.dot(ub, G) % 2
+ # print "x=", x
+ w = 1.0
+ for i in range(N):
+ w *= w_transition_prob(y[i], x[i], p)
+ # print "for u1=", uf, "prob=", w
+ prob += w
+ divider = 1.0 / (2 ** (N - 1))
+ return divider * prob
+
+
+def wn_split_channel(N, p):
+ G = get_Gn(N)
+ y = np.zeros((N, ), dtype=np.uint8)
+ n = 1
+ u = np.zeros((n + 1, ), dtype=np.uint8)
+
+ pn = w_split_prob(y, u, G, p)
+ # print "pn=", pn
+ z_params = []
+ c_params = []
+ for w in range(N):
+ nentries = (2 ** N) * (2 ** w)
+ print "for w=", w, " nentries=", nentries
+ w_probs = np.zeros((nentries, 2), dtype=float)
+ for y in range(2 ** N):
+ yb = unpack_byte(np.array([y, ]), N)
+ # print "\n\nyb", yb
+ for u in range(2 ** (w + 1)):
+ ub = unpack_byte(np.array([u, ]), w + 1)
+ # print "ub", ub
+ wp = w_split_prob(yb, ub, G, p)
+ ufixed = ub[0:-1]
+ yindex = y * (2 ** w) + pack_byte(ufixed)
+ uindex = ub[-1]
+ # print "y=", y, "ub", u, " prob=", wp, "index=[", yindex, ",
", uindex, "]"
+ w_probs[yindex][uindex] = wp
+ print "\n", np.sum(w_probs, axis=0)
+ z = bhattacharyya_parameter(w_probs)
+ z_params.append(z)
+ c = mutual_information(w_probs)
+ c_params.append(c)
+ print "W=", w, "Z=", z, "capacity=", c
+
+ return z_params, c_params
+
+
+def wminus(y0, y1, u0, p):
+ prob = 0.0
+ for i in range(2):
+ # print y0, y1, u0, i
+ u0u1 = (i + u0) % 2
+ w0 = w_transition_prob(y0, u0u1, p)
+ w1 = w_transition_prob(y1, i, p)
+ # print "w0=", w0, "\tw1=", w1
+ prob += w0 * w1
+ prob *= 0.5
+ return prob
+
+
+def wplus(y0, y1, u0, u1, p):
+ u0u1 = (u0 + u1) % 2
+ w0 = w_transition_prob(y0, u0u1, p)
+ w1 = w_transition_prob(y1, u1, p)
+ return 0.5 * w0 * w1
+
+
+def w2_split_channel(p):
+
+ wm_probs = np.zeros((4, 2), dtype=float)
+ for y0 in range(2):
+ for y1 in range(2):
+ for u0 in range(2):
+ wm = wminus(y0, y1, u0, p)
+ wm_probs[y0 * 2 + y1][u0] = wm
+
+ print wm_probs
+ print "W- Z parameter=", bhattacharyya_parameter(wm_probs)
+ print "I(W-)=", mutual_information(wm_probs)
+
+ chan_prob = 0.0
+ wp_probs = np.zeros((8, 2), dtype=float)
+ for y0 in range(2):
+ for y1 in range(2):
+ for u0 in range(2):
+ for u1 in range(2):
+ wp = wplus(y0, y1, u0, u1, p)
+ wp_probs[4 * y0 + 2 * y1 + u0][u1] = wp
+
+ print wp_probs
+ print "W+ Z parameter=", bhattacharyya_parameter(wp_probs)
+ print "I(W+)=", mutual_information(wp_probs)
+
+
+def plot_channel(probs, p, name):
+ nc = len(probs)
+ plt.plot(probs)
+ plt.grid()
+ plt.xlim((0, nc - 1))
+ plt.ylim((0.0, 1.0))
+ plt.xlabel('channel')
+ plt.ylabel('capacity')
+
+ fname = name, '_', p, '_', nc
+ print fname
+ fname = ''.join(str(n) for n in fname)
+ print fname
+ fname = fname.replace('.', '_')
+ print fname
+ fname = ''.join((fname, '.pdf'))
+ print fname
+ plt.savefig(fname)
+ # plt.show()
+
+
+def main():
+ np.set_printoptions(precision=4, linewidth=150)
+ print "POLAR code!"
+ w2 = calculate_w2_probs()
+ print w2
+
+ p = 0.1
+ n = 2
+ wn = calculate_Wn_probs(n, bsc_channel(p))
+
+ wi = calculate_channel_splitting_probs(wn, n)
+
+ w0 = bsc_channel(p)
+ print w0
+ iw = mutual_information(w0)
+ print "I(W)=", iw
+ print 1 - (p * np.log2(1.0 / p) + (1 - p) * np.log2(1.0 / (1 - p)))
+
+ print "Z param call"
+ bhattacharyya_parameter(w0)
+
+ print "old vs new encoder"
+
+ p = 0.1
+ w2_split_channel(p)
+ z, c = wn_split_channel(4, p)
+
+ eta = 0.3
+ nc = 1024
+ probs = calculate_bec_channel_capacities(eta, nc)
+ # plot_channel(probs, eta, 'bec_capacity')
+ # plot_channel(c, p, 'bsc_capacity')
+
+ nt = 2 ** 5
+
+
+ fG = get_Gn(nt)
+ # print fG
+ encoder = enc(nt, 4, np.arange(nt - 4, dtype=int))
+ # print encoder.get_gn()
+ comp = np.equal(fG, encoder.get_gn())
+
+ print "is equal?", np.all(comp)
+
+
+
+if __name__ == '__main__':
+ main()
diff --git a/gr-fec/python/fec/polar/testbed.py
b/gr-fec/python/fec/polar/testbed.py
new file mode 100755
index 0000000..f34f8a6
--- /dev/null
+++ b/gr-fec/python/fec/polar/testbed.py
@@ -0,0 +1,143 @@
+#!/usr/bin/env python
+#
+# Copyright 2015 Free Software Foundation, Inc.
+#
+# 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 as np
+from encoder import PolarEncoder
+from decoder import PolarDecoder
+
+import matplotlib.pyplot as plt
+
+
+def get_frozen_bit_position():
+ # frozenbitposition = np.array((0, 1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 16, 17,
18, 20, 24), dtype=int)
+ # frozenbitposition = np.array((0, 1, 2, 3, 4, 5, 8, 9), dtype=int)
+ frozenbitposition =
np.load('frozen_bit_positions_n256_k128_p0.11.npy').flatten()
+ print(frozenbitposition)
+ return frozenbitposition
+
+
+def test_enc_dec_chain():
+ ntests = 100
+ n = 256
+ k = n // 2
+ frozenbits = np.zeros(n - k)
+ frozenbitposition = get_frozen_bit_position()
+ for i in range(ntests):
+ bits = np.random.randint(2, size=k)
+ encoder = PolarEncoder(n, k, frozenbitposition, frozenbits)
+ decoder = PolarDecoder(n, k, frozenbitposition, frozenbits)
+ encoded = encoder.encode(bits)
+ rx = decoder.decode(encoded)
+ if not is_equal(bits, rx):
+ raise ValueError('Test #', i, 'failed, input and output differ',
bits, '!=', rx)
+ return
+
+
+def is_equal(first, second):
+ if not (first == second).all():
+ result = first == second
+ for i in range(len(result)):
+ print '{0:4}: {1:2} == {2:1} = {3}'.format(i, first[i], second[i],
result[i])
+ return False
+ return True
+
+
+def exact_value(la, lb):
+ return np.log((np.exp(la + lb) + 1) / (np.exp(la) + np.exp(lb)))
+
+
+def approx_value(la, lb):
+ return np.sign(la) * np.sign(lb) * np.minimum(np.abs(la), np.abs(lb))
+
+
+def test_1024_rate_1_code():
+ # effectively a Monte-Carlo simulation for channel polarization.
+ ntests = 10000
+ n = 256
+ k = n
+ transition_prob = 0.11
+ num_transitions = int(k * transition_prob)
+ frozenbits = np.zeros(n - k)
+ frozenbitposition = np.array((), dtype=int)
+ encoder = PolarEncoder(n, k, frozenbitposition, frozenbits)
+ decoder = PolarDecoder(n, k, frozenbitposition, frozenbits)
+
+ channel_counter = np.zeros(k)
+ possible_indices = np.arange(n, dtype=int)
+ for i in range(ntests):
+ bits = np.random.randint(2, size=k)
+ tx = encoder.encode(bits)
+ np.random.shuffle(possible_indices)
+ tx[possible_indices[0:num_transitions]] =
(tx[possible_indices[0:num_transitions]] + 1) % 2
+ rx = tx
+ recv = decoder.decode(rx)
+ channel_counter += (bits == recv)
+
+ print channel_counter
+ print(np.min(channel_counter), np.max(channel_counter))
+
+ np.save('channel_counter_' + str(ntests) + '.npy', channel_counter)
+
+
+def find_good_indices(res, nindices):
+ channel_counter = np.copy(res)
+ good_indices = np.zeros(channel_counter.size)
+
+ for i in range(nindices):
+ idx = np.argmax(channel_counter)
+ good_indices[idx] = 1
+ channel_counter[idx] = 0
+ return good_indices
+
+
+def channel_analysis():
+ ntests = 10000
+ filename = 'channel_counter_' + str(ntests) + '.npy'
+ channel_counter = np.load(filename)
+ print(np.min(channel_counter), np.max(channel_counter))
+ channel_counter[0] = np.min(channel_counter)
+ good_indices = find_good_indices(channel_counter, channel_counter.size //
2)
+ frozen_bit_positions = np.where(good_indices > 0)
+ print(frozen_bit_positions)
+ np.save('frozen_bit_positions_n256_k128_p0.11.npy', frozen_bit_positions)
+ good_indices *= 2000
+ good_indices += 4000
+
+
+ plt.plot(channel_counter)
+ plt.plot(good_indices)
+ plt.show()
+
+def main():
+ # n = 16
+ # k = 8
+ # frozenbits = np.zeros(n - k)
+ # frozenbitposition8 = np.array((0, 1, 2, 4), dtype=int)
+ # frozenbitposition = np.array((0, 1, 2, 3, 4, 5, 8, 9), dtype=int)
+ # print frozenbitposition
+
+ test_enc_dec_chain()
+
+ # test_1024_rate_1_code()
+
+ # channel_analysis()
+
+if __name__ == '__main__':
+ main()
\ No newline at end of file
- [Commit-gnuradio] [gnuradio] 10/25: polar: rebased and fixed branch for merge again., (continued)
- [Commit-gnuradio] [gnuradio] 10/25: polar: rebased and fixed branch for merge again., git, 2015/09/23
- [Commit-gnuradio] [gnuradio] 19/25: polar: rearranged polar_code_example flowgraph, git, 2015/09/23
- [Commit-gnuradio] [gnuradio] 21/25: Merge branch 'maint', git, 2015/09/23
- [Commit-gnuradio] [gnuradio] 05/25: polar: VOLK kernel integration, git, 2015/09/23
- [Commit-gnuradio] [gnuradio] 24/25: Merge remote-tracking branch 'gnuradio-wg-grc/master_grcwg', git, 2015/09/23
- [Commit-gnuradio] [gnuradio] 07/25: polar: updated polar code example GRC, git, 2015/09/23
- [Commit-gnuradio] [gnuradio] 02/25: polar: encoder implemented in C++, git, 2015/09/23
- [Commit-gnuradio] [gnuradio] 25/25: Merge remote-tracking branch 'tom/socis/master', git, 2015/09/23
- [Commit-gnuradio] [gnuradio] 06/25: polar: refarctoring and clean-up, git, 2015/09/23
- [Commit-gnuradio] [gnuradio] 03/25: polar: SC and SC List decoder implemented in C++, git, 2015/09/23
- [Commit-gnuradio] [gnuradio] 01/25: polar: Python test code for polar codes,
git <=
- [Commit-gnuradio] [gnuradio] 04/25: polar: Python implementation for channel construction, git, 2015/09/23
- [Commit-gnuradio] [gnuradio] 13/25: polar: updated documentation and examples, git, 2015/09/23
- [Commit-gnuradio] [gnuradio] 15/25: fec: polar: moving polar codes into ::code:: namespace like the others., git, 2015/09/23
- [Commit-gnuradio] [gnuradio] 16/25: fec: polar: whitespace fixes/tabs to go along with namespace change from last commit., git, 2015/09/23
- [Commit-gnuradio] [gnuradio] 14/25: polar: updated documentation and examples, git, 2015/09/23