diff options
| author | Johannes Demel <ufcsy@student.kit.edu> | 2015-05-28 12:39:36 +0200 |
|---|---|---|
| committer | Johannes Demel <ufcsy@student.kit.edu> | 2015-09-21 10:45:12 +0200 |
| commit | 8ac5a86f98d5a545cf65dda77c960cf273ee7aca (patch) | |
| tree | ce89f793a584e39842a13dedcd700a79f5f9aeed /gr-fec/python | |
| parent | a03893d51553eb74ac348ee88d4dd879b0a253c8 (diff) | |
polar: Python test code for polar codes
Diffstat (limited to 'gr-fec/python')
| -rw-r--r-- | gr-fec/python/fec/polar/README.md | 4 | ||||
| -rw-r--r-- | gr-fec/python/fec/polar/bit_reversed_polar_encoding_scheme.svg | 1160 | ||||
| -rw-r--r-- | gr-fec/python/fec/polar/channel_construction_bec.py | 110 | ||||
| -rwxr-xr-x | gr-fec/python/fec/polar/channel_construction_bsc.py | 206 | ||||
| -rw-r--r-- | gr-fec/python/fec/polar/common.py | 69 | ||||
| -rw-r--r-- | gr-fec/python/fec/polar/decoder.py | 278 | ||||
| -rw-r--r-- | gr-fec/python/fec/polar/encoder.py | 135 | ||||
| -rw-r--r-- | gr-fec/python/fec/polar/helper_functions.py | 84 | ||||
| -rw-r--r-- | gr-fec/python/fec/polar/polar_code.py | 493 | ||||
| -rwxr-xr-x | gr-fec/python/fec/polar/testbed.py | 143 |
10 files changed, 2682 insertions, 0 deletions
diff --git a/gr-fec/python/fec/polar/README.md b/gr-fec/python/fec/polar/README.md new file mode 100644 index 0000000000..2bd00dc3de --- /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 0000000000..84849c994e --- /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 0000000000..cd61f463f6 --- /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 0000000000..8cee7ceac5 --- /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 0000000000..b4b152de61 --- /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 0000000000..d74f1f9e1a --- /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 0000000000..6f87a22191 --- /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 0000000000..1d82457f17 --- /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 0000000000..f263c6d056 --- /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 0000000000..f34f8a6b96 --- /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 |
