added polar bit freezers, encoders and decoders

This commit is contained in:
Ahmet Inan 2021-07-10 23:11:05 +02:00
commit cc1f412f5d
8 changed files with 985 additions and 0 deletions

View file

@ -84,6 +84,27 @@ Encoder for [augmented Hadamard codes](https://en.wikipedia.org/wiki/Hadamard_co
[SIMD](https://en.wikipedia.org/wiki/SIMD) intra-frame accelerated [Low-density parity-check](https://en.wikipedia.org/wiki/Low-density_parity-check_code) layered decoder.
### [polar_freezer.hh](polar_freezer.hh)
Bit freezers for the construction of [polar codes](https://en.wikipedia.org/wiki/Polar_code_(coding_theory)).
* PolarFreezer: Constructs code for a given erasure probability without the need for storing nor sorting of erasure probabilities for the virtual channels.
* PolarCodeConst0: Constructs code by choosing the K best virtual channels computed from a given erasure probability.
### [polar_encoder.hh](polar_encoder.hh)
Encoders for [non-systematic and systematic](https://en.wikipedia.org/wiki/Systematic_code) [polar codes](https://en.wikipedia.org/wiki/Polar_code_(coding_theory)).
### [polar_decoder.hh](polar_decoder.hh)
Successive cancellation decoding of [polar codes](https://en.wikipedia.org/wiki/Polar_code_(coding_theory)).
### [polar_list_decoder.hh](polar_list_decoder.hh)
Successive cancellation [list decoding](https://en.wikipedia.org/wiki/List_decoding) of [polar codes](https://en.wikipedia.org/wiki/Polar_code_(coding_theory)).
List size depends on used SIMD type. Decoding performance of the fixed-point implementation is better with shorter codes, as the list size is larger but a deterioration can be observed with larger codes.
### [osd.hh](osd.hh)
Ordered statistics decoding allows the practical [soft-decision decoding](https://en.wikipedia.org/wiki/Soft-decision_decoder) of short to medium sized [linear codes](https://en.wikipedia.org/wiki/Linear_code).

99
polar_decoder.hh Normal file
View file

@ -0,0 +1,99 @@
/*
Successive cancellation decoding of polar codes
Copyright 2020 Ahmet Inan <inan@aicodix.de>
*/
#pragma once
#include <algorithm>
#include "polar_helper.hh"
namespace CODE {
template <typename TYPE, int M>
struct PolarTree
{
typedef PolarHelper<TYPE> PH;
static const int N = 1 << M;
static void decode(TYPE **message, TYPE *hard, TYPE *soft, const uint8_t *frozen)
{
for (int i = 0; i < N/2; ++i)
soft[i+N/2] = PH::prod(soft[i+N], soft[i+N/2+N]);
PolarTree<TYPE, M-1>::decode(message, hard, soft, frozen);
for (int i = 0; i < N/2; ++i)
soft[i+N/2] = PH::madd(hard[i], soft[i+N], soft[i+N/2+N]);
PolarTree<TYPE, M-1>::decode(message, hard+N/2, soft, frozen+N/2);
for (int i = 0; i < N/2; ++i)
hard[i] = PH::qmul(hard[i], hard[i+N/2]);
}
};
template <typename TYPE>
struct PolarTree<TYPE, 0>
{
typedef PolarHelper<TYPE> PH;
static void decode(TYPE **message, TYPE *hard, TYPE *soft, const uint8_t *frozen)
{
if (*frozen) {
*hard = PH::one();
} else {
*hard = PH::signum(soft[1]);
*(*message)++ = *hard;
}
}
};
template <typename TYPE, int MAX_M>
class PolarDecoder
{
typedef PolarHelper<TYPE> PH;
static const int MAX_N = 1 << MAX_M;
TYPE soft[2*MAX_N];
TYPE hard[MAX_N];
public:
void operator()(TYPE *message, const TYPE *codeword, const uint8_t *frozen, int level)
{
assert(level <= MAX_M);
int length = 1 << level;
for (int i = 0; i < length; ++i)
soft[length+i] = codeword[i];
switch (level) {
case 0: PolarTree<TYPE, 0>::decode(&message, hard, soft, frozen); break;
case 1: PolarTree<TYPE, 1>::decode(&message, hard, soft, frozen); break;
case 2: PolarTree<TYPE, 2>::decode(&message, hard, soft, frozen); break;
case 3: PolarTree<TYPE, 3>::decode(&message, hard, soft, frozen); break;
case 4: PolarTree<TYPE, 4>::decode(&message, hard, soft, frozen); break;
case 5: PolarTree<TYPE, 5>::decode(&message, hard, soft, frozen); break;
case 6: PolarTree<TYPE, 6>::decode(&message, hard, soft, frozen); break;
case 7: PolarTree<TYPE, 7>::decode(&message, hard, soft, frozen); break;
case 8: PolarTree<TYPE, 8>::decode(&message, hard, soft, frozen); break;
case 9: PolarTree<TYPE, 9>::decode(&message, hard, soft, frozen); break;
case 10: PolarTree<TYPE, 10>::decode(&message, hard, soft, frozen); break;
case 11: PolarTree<TYPE, 11>::decode(&message, hard, soft, frozen); break;
case 12: PolarTree<TYPE, 12>::decode(&message, hard, soft, frozen); break;
case 13: PolarTree<TYPE, 13>::decode(&message, hard, soft, frozen); break;
case 14: PolarTree<TYPE, 14>::decode(&message, hard, soft, frozen); break;
case 15: PolarTree<TYPE, 15>::decode(&message, hard, soft, frozen); break;
case 16: PolarTree<TYPE, 16>::decode(&message, hard, soft, frozen); break;
case 17: PolarTree<TYPE, 17>::decode(&message, hard, soft, frozen); break;
case 18: PolarTree<TYPE, 18>::decode(&message, hard, soft, frozen); break;
case 19: PolarTree<TYPE, 19>::decode(&message, hard, soft, frozen); break;
case 20: PolarTree<TYPE, 20>::decode(&message, hard, soft, frozen); break;
case 21: PolarTree<TYPE, 21>::decode(&message, hard, soft, frozen); break;
case 22: PolarTree<TYPE, 22>::decode(&message, hard, soft, frozen); break;
case 23: PolarTree<TYPE, 23>::decode(&message, hard, soft, frozen); break;
case 24: PolarTree<TYPE, 24>::decode(&message, hard, soft, frozen); break;
case 25: PolarTree<TYPE, 25>::decode(&message, hard, soft, frozen); break;
case 26: PolarTree<TYPE, 26>::decode(&message, hard, soft, frozen); break;
case 27: PolarTree<TYPE, 27>::decode(&message, hard, soft, frozen); break;
case 28: PolarTree<TYPE, 28>::decode(&message, hard, soft, frozen); break;
case 29: PolarTree<TYPE, 29>::decode(&message, hard, soft, frozen); break;
default: assert(false);
}
}
};
}

66
polar_encoder.hh Normal file
View file

@ -0,0 +1,66 @@
/*
Polar encoder for non-systematic and systematic codes
Copyright 2020 Ahmet Inan <inan@aicodix.de>
*/
#pragma once
#include "polar_helper.hh"
namespace CODE {
template <typename TYPE>
class PolarEncoder
{
typedef PolarHelper<TYPE> PH;
public:
void operator()(TYPE *codeword, const TYPE *message, const uint8_t *frozen, int level)
{
int length = 1 << level;
for (int i = 0; i < length; i += 2) {
TYPE msg0 = frozen[i] ? PH::one() : *message++;
TYPE msg1 = frozen[i+1] ? PH::one() : *message++;
codeword[i] = PH::qmul(msg0, msg1);
codeword[i+1] = msg1;
}
for (int h = 2; h < length; h *= 2)
for (int i = 0; i < length; i += 2 * h)
for (int j = i; j < i + h; ++j)
codeword[j] = PH::qmul(codeword[j], codeword[j+h]);
}
};
template <typename TYPE>
class PolarSysEnc
{
typedef PolarHelper<TYPE> PH;
public:
void operator()(TYPE *codeword, const TYPE *message, const uint8_t *frozen, int level)
{
int length = 1 << level;
for (int i = 0; i < length; i += 2) {
TYPE msg0 = frozen[i] ? PH::one() : *message++;
TYPE msg1 = frozen[i+1] ? PH::one() : *message++;
codeword[i] = PH::qmul(msg0, msg1);
codeword[i+1] = msg1;
}
for (int h = 2; h < length; h *= 2)
for (int i = 0; i < length; i += 2 * h)
for (int j = i; j < i + h; ++j)
codeword[j] = PH::qmul(codeword[j], codeword[j+h]);
for (int i = 0; i < length; i += 2) {
TYPE msg0 = frozen[i] ? PH::one() : codeword[i];
TYPE msg1 = frozen[i+1] ? PH::one() : codeword[i+1];
codeword[i] = PH::qmul(msg0, msg1);
codeword[i+1] = msg1;
}
for (int h = 2; h < length; h *= 2)
for (int i = 0; i < length; i += 2 * h)
for (int j = i; j < i + h; ++j)
codeword[j] = PH::qmul(codeword[j], codeword[j+h]);
}
};
}

67
polar_freezer.hh Normal file
View file

@ -0,0 +1,67 @@
/*
Bit freezers for polar codes
Copyright 2020 Ahmet Inan <inan@aicodix.de>
*/
#pragma once
#include <algorithm>
namespace CODE {
class PolarFreezer
{
static void freeze(uint8_t *bits, long double pe, long double th, int i, int h)
{
if (h) {
freeze(bits, pe * (2-pe), th, i, h/2);
freeze(bits, pe * pe, th, i+h, h/2);
} else {
bits[i] = pe > th;
}
}
public:
int operator()(uint8_t *frozen_bits, int level, long double erasure_probability = 0.5L, long double freezing_threshold = 0.5L)
{
int length = 1 << level;
freeze(frozen_bits, erasure_probability, freezing_threshold, 0, length / 2);
int K = 0;
for (int i = 0; i < length; ++i)
K += !frozen_bits[i];
return K;
}
};
template <int MAX_M>
class PolarCodeConst0
{
void compute(long double pe, int i, int h)
{
if (h) {
compute(pe * (2-pe), i, h/2);
compute(pe * pe, i+h, h/2);
} else {
prob[i] = pe;
}
}
long double prob[1<<MAX_M];
int index[1<<MAX_M];
public:
void operator()(uint8_t *frozen_bits, int level, int K, long double erasure_probability = std::exp(-1.L))
{
assert(level <= MAX_M);
int length = 1 << level;
compute(erasure_probability, 0, length / 2);
for (int i = 0; i < length; ++i)
index[i] = i;
std::nth_element(index, index+K, index+length, [this](int a, int b){ return prob[a] < prob[b]; });
for (int i = 0; i < K; ++i)
frozen_bits[index[i]] = 0;
for (int i = K; i < length; ++i)
frozen_bits[index[i]] = 1;
}
};
}

206
polar_helper.hh Normal file
View file

@ -0,0 +1,206 @@
/*
SIMD wrapper used by polar encoder and decoder
Copyright 2020 Ahmet Inan <inan@aicodix.de>
*/
#pragma once
#include "simd.hh"
namespace CODE {
template <typename TYPE>
struct PolarHelper
{
typedef TYPE PATH;
static TYPE one()
{
return 1;
}
static TYPE zero()
{
return 0;
}
static TYPE signum(TYPE v)
{
return (v > 0) - (v < 0);
}
template <typename IN>
static TYPE quant(IN in)
{
return in;
}
static TYPE qabs(TYPE a)
{
return std::abs(a);
}
static TYPE qmin(TYPE a, TYPE b)
{
return std::min(a, b);
}
static TYPE qadd(TYPE a, TYPE b)
{
return a + b;
}
static TYPE qmul(TYPE a, TYPE b)
{
return a * b;
}
static TYPE prod(TYPE a, TYPE b)
{
return signum(a) * signum(b) * qmin(qabs(a), qabs(b));
}
static TYPE madd(TYPE a, TYPE b, TYPE c)
{
return a * b + c;
}
};
template <typename VALUE, int WIDTH>
struct PolarHelper<SIMD<VALUE, WIDTH>>
{
typedef SIMD<VALUE, WIDTH> TYPE;
typedef VALUE PATH;
typedef SIMD<typename TYPE::uint_type, WIDTH> MAP;
static TYPE one()
{
return vdup<TYPE>(1);
}
static TYPE zero()
{
return vzero<TYPE>();
}
static TYPE signum(TYPE a)
{
return vsignum(a);
}
static TYPE qabs(TYPE a)
{
return vabs(a);
}
static TYPE qmin(TYPE a, TYPE b)
{
return vmin(a, b);
}
static TYPE qadd(TYPE a, TYPE b)
{
return vadd(a, b);
}
static TYPE qmul(TYPE a, TYPE b)
{
return vmul(a, b);
}
static TYPE prod(TYPE a, TYPE b)
{
return vmul(vmul(vsignum(a), vsignum(b)), vmin(vabs(a), vabs(b)));
}
static TYPE madd(TYPE a, TYPE b, TYPE c)
{
return vadd(vmul(a, b), c);
}
};
template <int WIDTH>
struct PolarHelper<SIMD<int8_t, WIDTH>>
{
typedef SIMD<int8_t, WIDTH> TYPE;
typedef int PATH;
typedef SIMD<uint8_t, WIDTH> MAP;
static TYPE one()
{
return vdup<TYPE>(1);
}
static TYPE zero()
{
return vzero<TYPE>();
}
static TYPE signum(TYPE a)
{
return vsignum(a);
}
static TYPE qabs(TYPE a)
{
return vqabs(a);
}
static TYPE qadd(TYPE a, TYPE b)
{
return vqadd(a, b);
}
static TYPE qmul(TYPE a, TYPE b)
{
#ifdef __ARM_NEON__
return vmul(a, b);
#else
return vsign(a, b);
#endif
}
static TYPE prod(TYPE a, TYPE b)
{
#ifdef __ARM_NEON__
return vmul(vmul(vsignum(a), vsignum(b)), vmin(vqabs(a), vqabs(b)));
#else
return vsign(vmin(vqabs(a), vqabs(b)), vsign(vsignum(a), b));
#endif
}
static TYPE madd(TYPE a, TYPE b, TYPE c)
{
#ifdef __ARM_NEON__
return vqadd(vmul(a, vmax(b, vdup<TYPE>(-127))), c);
#else
return vqadd(vsign(vmax(b, vdup<TYPE>(-127)), a), c);
#endif
}
};
template <>
struct PolarHelper<int8_t>
{
typedef int PATH;
static int8_t one()
{
return 1;
}
static int8_t zero()
{
return 0;
}
static int8_t signum(int8_t v)
{
return (v > 0) - (v < 0);
}
template <typename IN>
static int8_t quant(IN in)
{
return std::min<IN>(std::max<IN>(std::nearbyint(in), -128), 127);
}
static int8_t qabs(int8_t a)
{
return std::abs(std::max<int8_t>(a, -127));
}
static int8_t qmin(int8_t a, int8_t b)
{
return std::min(a, b);
}
static int8_t qadd(int8_t a, int8_t b)
{
return std::min<int16_t>(std::max<int16_t>(int16_t(a) + int16_t(b), -128), 127);
}
static int8_t qmul(int8_t a, int8_t b)
{
// return std::min<int16_t>(std::max<int16_t>(int16_t(a) * int16_t(b), -128), 127);
// only used for hard decision values anyway
return a * b;
}
static int8_t prod(int8_t a, int8_t b)
{
return signum(a) * signum(b) * qmin(qabs(a), qabs(b));
}
static int8_t madd(int8_t a, int8_t b, int8_t c)
{
return std::min<int16_t>(std::max<int16_t>(int16_t(a) * int16_t(b) + int16_t(c), -128), 127);
}
};
}

146
polar_list_decoder.hh Normal file
View file

@ -0,0 +1,146 @@
/*
Successive cancellation list decoding of polar codes
Copyright 2020 Ahmet Inan <inan@aicodix.de>
*/
#pragma once
#include <algorithm>
#include "polar_helper.hh"
namespace CODE {
template <typename TYPE, int M>
struct PolarListTree
{
typedef PolarHelper<TYPE> PH;
typedef typename PH::PATH PATH;
typedef typename PH::MAP MAP;
static const int N = 1 << M;
static MAP decode(PATH *metric, TYPE *message, MAP *maps, int *count, TYPE *hard, TYPE *soft, const uint8_t *frozen)
{
for (int i = 0; i < N/2; ++i)
soft[i+N/2] = PH::prod(soft[i+N], soft[i+N/2+N]);
MAP lmap = PolarListTree<TYPE, M-1>::decode(metric, message, maps, count, hard, soft, frozen);
for (int i = 0; i < N/2; ++i)
soft[i+N/2] = PH::madd(hard[i], vshuf(soft[i+N], lmap), vshuf(soft[i+N/2+N], lmap));
MAP rmap = PolarListTree<TYPE, M-1>::decode(metric, message, maps, count, hard+N/2, soft, frozen+N/2);
for (int i = 0; i < N/2; ++i)
hard[i] = PH::qmul(vshuf(hard[i], rmap), hard[i+N/2]);
return vshuf(lmap, rmap);
}
};
template <typename TYPE>
struct PolarListTree<TYPE, 0>
{
typedef PolarHelper<TYPE> PH;
typedef typename PH::PATH PATH;
typedef typename PH::MAP MAP;
static MAP decode(PATH *metric, TYPE *message, MAP *maps, int *count, TYPE *hard, TYPE *soft, const uint8_t *frozen)
{
MAP map;
TYPE hrd, sft = soft[1];
if (*frozen) {
for (int k = 0; k < TYPE::SIZE; ++k)
if (sft.v[k] < 0)
metric[k] -= sft.v[k];
hrd = PH::one();
for (int k = 0; k < TYPE::SIZE; ++k)
map.v[k] = k;
} else {
PATH fork[2*TYPE::SIZE];
for (int k = 0; k < TYPE::SIZE; ++k)
fork[k] = fork[k+TYPE::SIZE] = metric[k];
for (int k = 0; k < TYPE::SIZE; ++k)
if (sft.v[k] < 0)
fork[k] -= sft.v[k];
else
fork[k+TYPE::SIZE] += sft.v[k];
int perm[2*TYPE::SIZE];
for (int k = 0; k < 2*TYPE::SIZE; ++k)
perm[k] = k;
std::nth_element(perm, perm+TYPE::SIZE, perm+2*TYPE::SIZE, [fork](int a, int b){ return fork[a] < fork[b]; });
for (int k = 0; k < TYPE::SIZE; ++k)
metric[k] = fork[perm[k]];
for (int k = 0; k < TYPE::SIZE; ++k)
map.v[k] = perm[k] % TYPE::SIZE;
for (int k = 0; k < TYPE::SIZE; ++k)
hrd.v[k] = perm[k] < TYPE::SIZE ? 1 : -1;
message[*count] = hrd;
maps[*count] = map;
++*count;
}
*hard = hrd;
return map;
}
};
template <typename TYPE, int MAX_M>
class PolarListDecoder
{
typedef PolarHelper<TYPE> PH;
typedef typename TYPE::value_type VALUE;
typedef typename PH::PATH PATH;
typedef typename PH::MAP MAP;
static const int MAX_N = 1 << MAX_M;
TYPE soft[2*MAX_N];
TYPE hard[MAX_N];
MAP maps[MAX_N];
public:
void operator()(PATH *metric, TYPE *message, const VALUE *codeword, const uint8_t *frozen, int level)
{
assert(level <= MAX_M);
int count = 0;
metric[0] = 0;
for (int k = 1; k < TYPE::SIZE; ++k)
metric[k] = 1000;
int length = 1 << level;
for (int i = 0; i < length; ++i)
soft[length+i] = vdup<TYPE>(codeword[i]);
switch (level) {
case 0: PolarListTree<TYPE, 0>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 1: PolarListTree<TYPE, 1>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 2: PolarListTree<TYPE, 2>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 3: PolarListTree<TYPE, 3>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 4: PolarListTree<TYPE, 4>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 5: PolarListTree<TYPE, 5>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 6: PolarListTree<TYPE, 6>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 7: PolarListTree<TYPE, 7>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 8: PolarListTree<TYPE, 8>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 9: PolarListTree<TYPE, 9>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 10: PolarListTree<TYPE, 10>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 11: PolarListTree<TYPE, 11>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 12: PolarListTree<TYPE, 12>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 13: PolarListTree<TYPE, 13>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 14: PolarListTree<TYPE, 14>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 15: PolarListTree<TYPE, 15>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 16: PolarListTree<TYPE, 16>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 17: PolarListTree<TYPE, 17>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 18: PolarListTree<TYPE, 18>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 19: PolarListTree<TYPE, 19>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 20: PolarListTree<TYPE, 20>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 21: PolarListTree<TYPE, 21>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 22: PolarListTree<TYPE, 22>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 23: PolarListTree<TYPE, 23>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 24: PolarListTree<TYPE, 24>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 25: PolarListTree<TYPE, 25>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 26: PolarListTree<TYPE, 26>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 27: PolarListTree<TYPE, 27>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 28: PolarListTree<TYPE, 28>::decode(metric, message, maps, &count, hard, soft, frozen); break;
case 29: PolarListTree<TYPE, 29>::decode(metric, message, maps, &count, hard, soft, frozen); break;
default: assert(false);
}
MAP acc = maps[count-1];
for (int i = count-2; i >= 0; --i) {
message[i] = vshuf(message[i], acc);
acc = vshuf(maps[i], acc);
}
}
};
}

View file

@ -0,0 +1,202 @@
/*
Regression Test for the Polar Encoder and List Decoders
Copyright 2020 Ahmet Inan <inan@aicodix.de>
*/
#include <limits>
#include <random>
#include <chrono>
#include <cassert>
#include <iomanip>
#include <iostream>
#include <algorithm>
#include <functional>
#include "polar_helper.hh"
#include "polar_list_decoder.hh"
#include "polar_encoder.hh"
#include "polar_freezer.hh"
int main()
{
const int M = 16;
const int N = 1 << M;
const bool systematic = true;
#if 0
typedef int8_t code_type;
double SCALE = 2;
#else
typedef float code_type;
double SCALE = 1;
#endif
#ifdef __AVX2__
const int SIZEOF_SIMD = 32;
#else
const int SIZEOF_SIMD = 16;
#endif
const int SIMD_WIDTH = SIZEOF_SIMD / sizeof(code_type);
typedef SIMD<code_type, SIMD_WIDTH> simd_type;
std::random_device rd;
typedef std::default_random_engine generator;
typedef std::uniform_int_distribution<int> distribution;
auto data = std::bind(distribution(0, 1), generator(rd()));
auto frozen = new uint8_t[N];
auto codeword = new code_type[N];
auto temp = new simd_type[N];
long double erasure_probability = 1. / 3.;
int K = (1 - erasure_probability) * N;
double design_SNR = 10 * std::log10(-std::log(erasure_probability));
std::cerr << "design SNR: " << design_SNR << std::endl;
if (0) {
CODE::PolarFreezer freeze;
long double freezing_threshold = 0 ? 0.5 : std::numeric_limits<float>::epsilon();
K = freeze(frozen, M, erasure_probability, freezing_threshold);
} else {
auto freeze = new CODE::PolarCodeConst0<M>;
std::cerr << "sizeof(PolarCodeConst0<M>) = " << sizeof(CODE::PolarCodeConst0<M>) << std::endl;
double better_SNR = design_SNR + 1.59175;
std::cerr << "better SNR: " << better_SNR << std::endl;
long double probability = std::exp(-pow(10.0, better_SNR / 10));
(*freeze)(frozen, M, K, probability);
delete freeze;
}
std::cerr << "Polar(" << N << ", " << K << ")" << std::endl;
auto message = new code_type[K];
auto decoded = new simd_type[K];
CODE::PolarHelper<simd_type>::PATH metric[SIMD_WIDTH];
std::cerr << "sizeof(PolarListDecoder<simd_type, M>) = " << sizeof(CODE::PolarListDecoder<simd_type, M>) << std::endl;
auto decode = new CODE::PolarListDecoder<simd_type, M>;
auto orig = new code_type[N];
auto noisy = new code_type[N];
auto symb = new double[N];
double low_SNR = std::floor(design_SNR-3);
double high_SNR = std::ceil(design_SNR+5);
double min_SNR = high_SNR, max_mbs = 0;
int count = 0;
std::cerr << "SNR BER Mbit/s Eb/N0" << std::endl;
for (double SNR = low_SNR; count <= 3 && SNR <= high_SNR; SNR += 0.1, ++count) {
//double mean_signal = 0;
double sigma_signal = 1;
double mean_noise = 0;
double sigma_noise = std::sqrt(sigma_signal * sigma_signal / (2 * std::pow(10, SNR / 10)));
typedef std::normal_distribution<double> normal;
auto awgn = std::bind(normal(mean_noise, sigma_noise), generator(rd()));
int64_t awgn_errors = 0;
int64_t quantization_erasures = 0;
int64_t uncorrected_errors = 0;
int64_t ambiguity_erasures = 0;
double avg_mbs = 0;
int64_t loops = 0;
while (uncorrected_errors < 1000 && ++loops < 100) {
for (int i = 0; i < K; ++i)
message[i] = 1 - 2 * data();
if (systematic) {
CODE::PolarSysEnc<code_type> sysenc;
sysenc(codeword, message, frozen, M);
for (int i = 0, j = 0; i < N; ++i)
if (!frozen[i])
assert(codeword[i] == message[j++]);
} else {
CODE::PolarEncoder<code_type> encode;
encode(codeword, message, frozen, M);
}
for (int i = 0; i < N; ++i)
orig[i] = codeword[i];
for (int i = 0; i < N; ++i)
symb[i] = codeword[i];
for (int i = 0; i < N; ++i)
symb[i] += awgn();
// $LLR=log(\frac{p(x=+1|y)}{p(x=-1|y)})$
// $p(x|\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$
double DIST = 2; // BPSK
double fact = SCALE * DIST / (sigma_noise * sigma_noise);
for (int i = 0; i < N; ++i)
codeword[i] = CODE::PolarHelper<code_type>::quant(fact * symb[i]);
for (int i = 0; i < N; ++i)
noisy[i] = codeword[i];
auto start = std::chrono::system_clock::now();
(*decode)(metric, decoded, codeword, frozen, M);
auto end = std::chrono::system_clock::now();
auto usec = std::chrono::duration_cast<std::chrono::microseconds>(end - start);
double mbs = (double)K / usec.count();
avg_mbs += mbs;
if (systematic) {
CODE::PolarEncoder<simd_type> encode;
encode(temp, decoded, frozen, M);
for (int i = 0, j = 0; i < N; ++i)
if (!frozen[i])
decoded[j++] = temp[i];
}
int best = 0;
if (1) {
for (int k = 0; k < SIMD_WIDTH; ++k)
if (metric[k] < metric[best])
best = k;
} else {
int errs[SIMD_WIDTH] = { 0 };
for (int i = 0; i < K; ++i)
for (int k = 0; k < SIMD_WIDTH; ++k)
errs[k] += message[i] != decoded[i].v[k];
for (int k = 0; k < SIMD_WIDTH; ++k)
if (errs[k] < errs[best])
best = k;
}
for (int i = 0; i < N; ++i)
awgn_errors += noisy[i] * (orig[i] < 0);
for (int i = 0; i < N; ++i)
quantization_erasures += !noisy[i];
for (int i = 0; i < K; ++i)
uncorrected_errors += decoded[i].v[best] * message[i] <= 0;
for (int i = 0; i < K; ++i)
ambiguity_erasures += !decoded[i].v[best];
}
avg_mbs /= loops;
max_mbs = std::max(max_mbs, avg_mbs);
double bit_error_rate = (double)uncorrected_errors / (double)(K * loops);
if (!uncorrected_errors)
min_SNR = std::min(min_SNR, SNR);
else
count = 0;
int MOD_BITS = 1; // BPSK
double code_rate = (double)K / (double)N;
double spectral_efficiency = code_rate * MOD_BITS;
double EbN0 = 10 * std::log10(sigma_signal * sigma_signal / (spectral_efficiency * 2 * sigma_noise * sigma_noise));
if (0) {
std::cerr << SNR << " Es/N0 => AWGN with standard deviation of " << sigma_noise << " and mean " << mean_noise << std::endl;
std::cerr << EbN0 << " Eb/N0, using spectral efficiency of " << spectral_efficiency << " from " << code_rate << " code rate and " << MOD_BITS << " bits per symbol." << std::endl;
std::cerr << awgn_errors << " errors caused by AWGN." << std::endl;
std::cerr << quantization_erasures << " erasures caused by quantization." << std::endl;
std::cerr << uncorrected_errors << " errors uncorrected." << std::endl;
std::cerr << ambiguity_erasures << " ambiguity erasures." << std::endl;
std::cerr << bit_error_rate << " bit error rate." << std::endl;
std::cerr << avg_mbs << " megabit per second." << std::endl;
} else {
std::cout << SNR << " " << bit_error_rate << " " << avg_mbs << " " << EbN0 << std::endl;
}
}
std::cerr << "QEF at: " << min_SNR << " SNR, speed: " << max_mbs << " Mb/s." << std::endl;
double QEF_SNR = design_SNR + 0.2;
assert(min_SNR < QEF_SNR);
std::cerr << "Polar list regression test passed!" << std::endl;
return 0;
}

View file

@ -0,0 +1,178 @@
/*
Regression Test for the Polar Encoder and Decoders
Copyright 2020 Ahmet Inan <inan@aicodix.de>
*/
#include <limits>
#include <random>
#include <chrono>
#include <cassert>
#include <iomanip>
#include <iostream>
#include <algorithm>
#include <functional>
#include "polar_helper.hh"
#include "polar_decoder.hh"
#include "polar_encoder.hh"
#include "polar_freezer.hh"
int main()
{
const int M = 20;
const int N = 1 << M;
const bool systematic = true;
#if 1
typedef int8_t code_type;
double SCALE = 2;
#else
typedef float code_type;
double SCALE = 1;
#endif
std::random_device rd;
typedef std::default_random_engine generator;
typedef std::uniform_int_distribution<int> distribution;
auto data = std::bind(distribution(0, 1), generator(rd()));
auto frozen = new uint8_t[N];
auto codeword = new code_type[N];
auto temp = new code_type[N];
long double erasure_probability = 1. / 3.;
int K = (1 - erasure_probability) * N;
double design_SNR = 10 * std::log10(-std::log(erasure_probability));
std::cerr << "design SNR: " << design_SNR << std::endl;
if (0) {
CODE::PolarFreezer freeze;
long double freezing_threshold = 0 ? 0.5 : std::numeric_limits<float>::epsilon();
K = freeze(frozen, M, erasure_probability, freezing_threshold);
} else {
auto freeze = new CODE::PolarCodeConst0<M>;
std::cerr << "sizeof(PolarCodeConst0<M>) = " << sizeof(CODE::PolarCodeConst0<M>) << std::endl;
double better_SNR = design_SNR + 0.5;//1.59175;
std::cerr << "better SNR: " << better_SNR << std::endl;
long double probability = std::exp(-pow(10.0, better_SNR / 10));
(*freeze)(frozen, M, K, probability);
delete freeze;
}
std::cerr << "Polar(" << N << ", " << K << ")" << std::endl;
auto message = new code_type[K];
auto decoded = new code_type[K];
std::cerr << "sizeof(PolarDecoder<code_type, M>) = " << sizeof(CODE::PolarDecoder<code_type, M>) << std::endl;
auto decode = new CODE::PolarDecoder<code_type, M>;
auto orig = new code_type[N];
auto noisy = new code_type[N];
auto symb = new double[N];
double low_SNR = std::floor(design_SNR-3);
double high_SNR = std::ceil(design_SNR+5);
double min_SNR = high_SNR, max_mbs = 0;
int count = 0;
std::cerr << "SNR BER Mbit/s Eb/N0" << std::endl;
for (double SNR = low_SNR; count <= 3 && SNR <= high_SNR; SNR += 0.1, ++count) {
//double mean_signal = 0;
double sigma_signal = 1;
double mean_noise = 0;
double sigma_noise = std::sqrt(sigma_signal * sigma_signal / (2 * std::pow(10, SNR / 10)));
typedef std::normal_distribution<double> normal;
auto awgn = std::bind(normal(mean_noise, sigma_noise), generator(rd()));
int64_t awgn_errors = 0;
int64_t quantization_erasures = 0;
int64_t uncorrected_errors = 0;
int64_t ambiguity_erasures = 0;
double avg_mbs = 0;
int64_t loops = 0;
while (uncorrected_errors < 1000 && ++loops < 100) {
for (int i = 0; i < K; ++i)
message[i] = 1 - 2 * data();
if (systematic) {
CODE::PolarSysEnc<code_type> sysenc;
sysenc(codeword, message, frozen, M);
for (int i = 0, j = 0; i < N; ++i)
if (!frozen[i])
assert(codeword[i] == message[j++]);
} else {
CODE::PolarEncoder<code_type> encode;
encode(codeword, message, frozen, M);
}
for (int i = 0; i < N; ++i)
orig[i] = codeword[i];
for (int i = 0; i < N; ++i)
symb[i] = codeword[i];
for (int i = 0; i < N; ++i)
symb[i] += awgn();
// $LLR=log(\frac{p(x=+1|y)}{p(x=-1|y)})$
// $p(x|\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$
double DIST = 2; // BPSK
double fact = SCALE * DIST / (sigma_noise * sigma_noise);
for (int i = 0; i < N; ++i)
codeword[i] = CODE::PolarHelper<code_type>::quant(fact * symb[i]);
for (int i = 0; i < N; ++i)
noisy[i] = codeword[i];
auto start = std::chrono::system_clock::now();
(*decode)(decoded, codeword, frozen, M);
auto end = std::chrono::system_clock::now();
auto usec = std::chrono::duration_cast<std::chrono::microseconds>(end - start);
double mbs = (double)K / usec.count();
avg_mbs += mbs;
if (systematic) {
CODE::PolarEncoder<code_type> encode;
encode(temp, decoded, frozen, M);
for (int i = 0, j = 0; i < N; ++i)
if (!frozen[i])
decoded[j++] = temp[i];
}
for (int i = 0; i < N; ++i)
awgn_errors += noisy[i] * (orig[i] < 0);
for (int i = 0; i < N; ++i)
quantization_erasures += !noisy[i];
for (int i = 0; i < K; ++i)
uncorrected_errors += decoded[i] * message[i] <= 0;
for (int i = 0; i < K; ++i)
ambiguity_erasures += !decoded[i];
}
avg_mbs /= loops;
max_mbs = std::max(max_mbs, avg_mbs);
double bit_error_rate = (double)uncorrected_errors / (double)(K * loops);
if (!uncorrected_errors)
min_SNR = std::min(min_SNR, SNR);
else
count = 0;
int MOD_BITS = 1; // BPSK
double code_rate = (double)K / (double)N;
double spectral_efficiency = code_rate * MOD_BITS;
double EbN0 = 10 * std::log10(sigma_signal * sigma_signal / (spectral_efficiency * 2 * sigma_noise * sigma_noise));
if (0) {
std::cerr << SNR << " Es/N0 => AWGN with standard deviation of " << sigma_noise << " and mean " << mean_noise << std::endl;
std::cerr << EbN0 << " Eb/N0, using spectral efficiency of " << spectral_efficiency << " from " << code_rate << " code rate and " << MOD_BITS << " bits per symbol." << std::endl;
std::cerr << awgn_errors << " errors caused by AWGN." << std::endl;
std::cerr << quantization_erasures << " erasures caused by quantization." << std::endl;
std::cerr << uncorrected_errors << " errors uncorrected." << std::endl;
std::cerr << ambiguity_erasures << " ambiguity erasures." << std::endl;
std::cerr << bit_error_rate << " bit error rate." << std::endl;
std::cerr << avg_mbs << " megabit per second." << std::endl;
} else {
std::cout << SNR << " " << bit_error_rate << " " << avg_mbs << " " << EbN0 << std::endl;
}
}
std::cerr << "QEF at: " << min_SNR << " SNR, speed: " << max_mbs << " Mb/s." << std::endl;
double QEF_SNR = design_SNR + 0.2;
assert(min_SNR < QEF_SNR);
std::cerr << "Polar regression test passed!" << std::endl;
return 0;
}