added Reed Solomon error correction

lifted from my FEC project:
https://github.com/xdsopl/FEC
This commit is contained in:
Ahmet Inan 2018-09-20 18:04:23 +02:00
commit 01db8f516e
10 changed files with 846 additions and 1 deletions

View file

@ -33,5 +33,14 @@ When dealing with unaligned and arbitrary-bit-sized elements in a data stream, t
### [galois_field.hh](galois_field.hh)
We have to thank [Évariste Galois](https://en.wikipedia.org/wiki/%C3%89variste_Galois)' contribution of the [Finite field](https://en.wikipedia.org/wiki/Finite_field) to mathematics, which laid the cornerstone for a variety of applications that we take for granted today.
We have to thank [Évariste Galois](https://en.wikipedia.org/wiki/%C3%89variste_Galois) for his contribution of the [Finite field](https://en.wikipedia.org/wiki/Finite_field) to mathematics, which laid the cornerstone for a variety of applications that we take for granted today.
One of them is [ReedSolomon error correction](https://en.wikipedia.org/wiki/Reed%E2%80%93Solomon_error_correction):
### [reed_solomon_error_correction.hh](reed_solmon_error_correction.hh)
Implemented are the following Encoders and Decoders:
* [reed_solomon_encoder.hh](reed_solomon_encoder.hh)
* [reed_solomon_decoder.hh](reed_solomon_decoder.hh)
* [bose_chaudhuri_hocquenghem_encoder.hh](bose_chaudhuri_hocquenghem_encoder.hh)
* [bose_chaudhuri_hocquenghem_decoder.hh](bose_chaudhuri_hocquenghem_decoder.hh)

View file

@ -0,0 +1,77 @@
/*
Bose Chaudhuri Hocquenghem Decoder
Copyright 2018 Ahmet Inan <inan@aicodix.de>
*/
#ifndef BOSE_CHAUDHURI_HOCQUENGHEM_DECODER_HH
#define BOSE_CHAUDHURI_HOCQUENGHEM_DECODER_HH
#include "reed_solomon_error_correction.hh"
namespace CODE {
template <int NR, int FCR, int K, typename GF>
class BoseChaudhuriHocquenghemDecoder
{
public:
typedef typename GF::value_type value_type;
typedef typename GF::ValueType ValueType;
typedef typename GF::IndexType IndexType;
static const int N = GF::N, NP = N - K;
BoseChaudhuriHocquenghemDecoder() {}
int compute_syndromes(ValueType *code, ValueType *syndromes)
{
// $syndromes_i = code(pe^{FCR+i})$
for (int i = 0; i < NR; ++i)
syndromes[i] = code[0];
for (int j = 1; j < N; ++j) {
IndexType root(FCR), pe(1);
for (int i = 0; i < NR; ++i) {
syndromes[i] = fma(root, syndromes[i], code[j]);
root *= pe;
}
}
int nonzero = 0;
for (int i = 0; i < NR; ++i)
nonzero += !!syndromes[i];
return nonzero;
}
int operator()(ValueType *code, IndexType *erasures = 0, int erasures_count = 0)
{
assert(0 <= erasures_count && erasures_count <= NR);
if (0) {
for (int i = 0; i < erasures_count; ++i)
code[(int)erasures[i]] = ValueType(0);
}
ValueType syndromes[NR];
if (!compute_syndromes(code, syndromes))
return 0;
IndexType locations[NR];
ValueType magnitudes[NR];
int count = ReedSolomonErrorCorrection<NR, FCR, GF>::algorithm(syndromes, locations, magnitudes, erasures, erasures_count);
if (count <= 0)
return count;
for (int i = 0; i < count; ++i)
if (1 < (int)magnitudes[i])
return -1;
for (int i = 0; i < count; ++i)
code[(int)locations[i]] += magnitudes[i];
int corrections_count = 0;
for (int i = 0; i < count; ++i)
corrections_count += !!magnitudes[i];
return corrections_count;
}
int operator()(value_type *code, value_type *erasures = 0, int erasures_count = 0)
{
return (*this)(reinterpret_cast<ValueType *>(code), reinterpret_cast<IndexType *>(erasures), erasures_count);
}
int compute_syndromes(value_type *code, value_type *syndromes)
{
return compute_syndromes(reinterpret_cast<ValueType *>(code), reinterpret_cast<ValueType *>(syndromes));
}
};
}
#endif

View file

@ -0,0 +1,86 @@
/*
Bose Chaudhuri Hocquenghem Encoder
Copyright 2018 Ahmet Inan <inan@aicodix.de>
*/
#ifndef BOSE_CHAUDHURI_HOCQUENGHEM_ENCODER_HH
#define BOSE_CHAUDHURI_HOCQUENGHEM_ENCODER_HH
#include <initializer_list>
namespace CODE {
template <int NR, int FCR, int K, typename GF>
class BoseChaudhuriHocquenghemEncoder
{
public:
typedef typename GF::value_type value_type;
typedef typename GF::ValueType ValueType;
typedef typename GF::IndexType IndexType;
static const int N = GF::N, NP = N - K;
ValueType generator[NP+1];
BoseChaudhuriHocquenghemEncoder(std::initializer_list<int> minimal_polynomials)
{
// $generator(x) = \prod_i(minpoly_i(x))$
int generator_degree = 1;
generator[0] = ValueType(1);
for (int i = 1; i <= NP; ++i)
generator[i] = ValueType(0);
for (auto m: minimal_polynomials) {
assert(0 < m && m < 1<<(GF::M+1));
int m_degree = GF::M;
while (!(m>>m_degree))
--m_degree;
assert(generator_degree + m_degree <= NP + 1);
for (int i = generator_degree; i >= 0; --i) {
if (!generator[i])
continue;
generator[i] = ValueType(m&1);
for (int j = 1; j <= m_degree; ++j)
generator[i+j] += ValueType((m>>j)&1);
}
generator_degree += m_degree;
}
assert(generator_degree == NP + 1);
if (0) {
IndexType root(FCR), pe(1);
for (int i = 0; i < NR; ++i) {
ValueType tmp(generator[NP]);
for (int j = 1; j <= NP; ++j)
tmp = fma(root, tmp, generator[NP-j]);
assert(!tmp);
root *= pe;
}
std::cout << "generator =";
for (int i = 0; i <= NP; ++i)
std::cout << " " << (int)generator[i];
std::cout << std::endl;
}
}
void operator()(ValueType *code)
{
// $code = data * x^{NP} + (data * x^{NP}) \mod{generator}$
for (int i = 0; i < NP; ++i)
code[K+i] = ValueType(0);
for (int i = 0; i < K; ++i) {
if (code[i] != code[K]) {
for (int j = 1; j < NP; ++j)
code[K+j-1] = generator[NP-j] + code[K+j];
code[N-1] = generator[0];
} else {
for (int j = 1; j < NP; ++j)
code[K+j-1] = code[K+j];
code[N-1] = ValueType(0);
}
}
}
void operator()(value_type *code)
{
(*this)(reinterpret_cast<ValueType *>(code));
}
};
}
#endif

74
reed_solomon_decoder.hh Normal file
View file

@ -0,0 +1,74 @@
/*
Reed Solomon Decoder
Copyright 2018 Ahmet Inan <inan@aicodix.de>
*/
#ifndef REED_SOLOMON_DECODER_HH
#define REED_SOLOMON_DECODER_HH
#include "reed_solomon_error_correction.hh"
namespace CODE {
template <int NR, int FCR, typename GF>
class ReedSolomonDecoder
{
public:
typedef typename GF::value_type value_type;
typedef typename GF::ValueType ValueType;
typedef typename GF::IndexType IndexType;
static const int N = GF::N, K = N - NR;
ReedSolomonDecoder() {}
int compute_syndromes(ValueType *code, ValueType *syndromes)
{
// $syndromes_i = code(pe^{FCR+i})$
for (int i = 0; i < NR; ++i)
syndromes[i] = code[0];
for (int j = 1; j < N; ++j) {
IndexType root(FCR), pe(1);
for (int i = 0; i < NR; ++i) {
syndromes[i] = fma(root, syndromes[i], code[j]);
root *= pe;
}
}
int nonzero = 0;
for (int i = 0; i < NR; ++i)
nonzero += !!syndromes[i];
return nonzero;
}
int operator()(ValueType *code, IndexType *erasures = 0, int erasures_count = 0)
{
assert(0 <= erasures_count && erasures_count <= NR);
if (0) {
for (int i = 0; i < erasures_count; ++i)
code[(int)erasures[i]] = ValueType(0);
}
ValueType syndromes[NR];
if (!compute_syndromes(code, syndromes))
return 0;
IndexType locations[NR];
ValueType magnitudes[NR];
int count = ReedSolomonErrorCorrection<NR, FCR, GF>::algorithm(syndromes, locations, magnitudes, erasures, erasures_count);
if (count <= 0)
return count;
for (int i = 0; i < count; ++i)
code[(int)locations[i]] += magnitudes[i];
int corrections_count = 0;
for (int i = 0; i < count; ++i)
corrections_count += !!magnitudes[i];
return corrections_count;
}
int operator()(value_type *code, value_type *erasures = 0, int erasures_count = 0)
{
return (*this)(reinterpret_cast<ValueType *>(code), reinterpret_cast<IndexType *>(erasures), erasures_count);
}
int compute_syndromes(value_type *code, value_type *syndromes)
{
return compute_syndromes(reinterpret_cast<ValueType *>(code), reinterpret_cast<ValueType *>(syndromes));
}
};
}
#endif

78
reed_solomon_encoder.hh Normal file
View file

@ -0,0 +1,78 @@
/*
Reed Solomon Encoder
Copyright 2018 Ahmet Inan <inan@aicodix.de>
*/
#ifndef REED_SOLOMON_ENCODER_HH
#define REED_SOLOMON_ENCODER_HH
namespace CODE {
template <int NR, int FCR, typename GF>
class ReedSolomonEncoder
{
public:
typedef typename GF::value_type value_type;
typedef typename GF::ValueType ValueType;
typedef typename GF::IndexType IndexType;
static const int N = GF::N, K = N - NR;
IndexType generator[NR+1];
ReedSolomonEncoder()
{
// $generator = \prod_{i=0}^{NR}(x-pe^{FCR+i})$
ValueType tmp[NR+1];
IndexType root(FCR), pe(1);
for (int i = 0; i < NR; ++i) {
tmp[i] = ValueType(1);
for (int j = i; j > 0; --j)
tmp[j] = fma(root, tmp[j], tmp[j-1]);
tmp[0] *= root;
root *= pe;
}
tmp[NR] = ValueType(1);
if (0) {
std::cout << "generator = ";
for (int i = NR; i > 0; --i) {
if (!tmp[i])
continue;
if (tmp[i] != ValueType(1))
std::cout << (int)tmp[i] << "*";
std::cout << "x";
if (i != 1)
std::cout << "^" << i;
std::cout << " + ";
}
std::cout << (int)tmp[0] << std::endl;
}
for (int i = 0; i <= NR; ++i)
generator[i] = index(tmp[i]);
}
void operator()(ValueType *code)
{
// $code = data * x^{NR} + (data * x^{NR}) \mod{generator}$
for (int i = 0; i < NR; ++i)
code[K+i] = ValueType(0);
for (int i = 0; i < K; ++i) {
ValueType feedback = code[i] + code[K];
if (feedback) {
IndexType fb = index(feedback);
for (int j = 1; j < NR; ++j)
code[K+j-1] = fma(fb, generator[NR-j], code[K+j]);
code[N-1] = value(generator[0] * fb);
} else {
for (int j = 1; j < NR; ++j)
code[K+j-1] = code[K+j];
code[N-1] = ValueType(0);
}
}
}
void operator()(value_type *code)
{
(*this)(reinterpret_cast<ValueType *>(code));
}
};
}
#endif

View file

@ -0,0 +1,252 @@
/*
Reed Solomon Error Correction
Copyright 2018 Ahmet Inan <inan@aicodix.de>
*/
#ifndef REED_SOLOMON_ERROR_CORRECTION_HH
#define REED_SOLOMON_ERROR_CORRECTION_HH
namespace CODE {
namespace RS {
template <int NR, typename GF>
struct Chien
{
typedef typename GF::value_type value_type;
typedef typename GF::ValueType ValueType;
typedef typename GF::IndexType IndexType;
static const int N = GF::N, K = N - NR;
static int search(ValueType *locator, int locator_degree, IndexType *locations)
{
ValueType tmp[locator_degree+1];
for (int i = 0; i <= locator_degree; ++i)
tmp[i] = locator[i];
int count = 0;
for (int i = 0; i < N; ++i) {
ValueType sum(tmp[0]);
for (int j = 1; j <= locator_degree; ++j)
sum += tmp[j] *= IndexType(j);
if (!sum)
locations[count++] = IndexType(i);
}
return count;
}
};
template <int NR, typename GF>
struct FindLocations
{
typedef typename GF::value_type value_type;
typedef typename GF::ValueType ValueType;
typedef typename GF::IndexType IndexType;
static int search(ValueType *locator, int locator_degree, IndexType *locations)
{
if (locator_degree == 1) {
locations[0] = (index(locator[0]) / index(locator[1])) / IndexType(1);
return 1;
}
#if 0
if (locator_degree == 2) {
if (!locator[1] || !locator[0])
return 0;
ValueType a(locator[2]), b(locator[1]), c(locator[0]);
ValueType ba(b/a), R(Artin_Schreier_imap(a*c/(b*b)));
if (!R)
return 0;
locations[0] = index(ba * R) / IndexType(1);
locations[1] = index(ba * R + ba) / IndexType(1);
return 2;
}
#endif
return Chien<NR, GF>::search(locator, locator_degree, locations);
}
};
template <int NR, int FCR, typename GF>
struct Forney
{
typedef typename GF::value_type value_type;
typedef typename GF::ValueType ValueType;
typedef typename GF::IndexType IndexType;
static const int N = GF::N, K = N - NR;
static int compute_evaluator(ValueType *syndromes, ValueType *locator, int locator_degree, ValueType *evaluator)
{
// $evaluator = (syndromes * locator) \bmod{x^{NR}}$
int tmp = std::min(locator_degree, NR-1);
int degree = -1;
for (int i = 0; i <= tmp; ++i) {
evaluator[i] = syndromes[i] * locator[0];
for (int j = 1; j <= i; ++j)
evaluator[i] += syndromes[i-j] * locator[j];
if (evaluator[i])
degree = i;
}
return degree;
}
static void compute_magnitudes(ValueType *locator, IndexType *locations, int count, ValueType *evaluator, int evaluator_degree, ValueType *magnitudes)
{
// $magnitude = root^{FCR-1} * \frac{evaluator(root)}{locator'(root)}$
for (int i = 0; i < count; ++i) {
IndexType root(locations[i] * IndexType(1)), tmp(root);
ValueType eval(evaluator[0]);
for (int j = 1; j <= evaluator_degree; ++j) {
eval += evaluator[j] * tmp;
tmp *= root;
}
if (!eval) {
magnitudes[i] = ValueType(0);
continue;
}
ValueType deriv(locator[1]);
IndexType root2(root * root), tmp2(root2);
for (int j = 3; j <= count; j += 2) {
deriv += locator[j] * tmp2;
tmp2 *= root2;
}
IndexType magnitude(index(eval) / index(deriv));
if (FCR == 0)
magnitude /= root;
if (FCR > 1)
for (int j = 1; j < FCR; ++j)
magnitude *= root;
magnitudes[i] = value(magnitude);
}
}
static int algorithm(ValueType *syndromes, ValueType *locator, IndexType *locations, int count, ValueType *evaluator, ValueType *magnitudes)
{
int evaluator_degree = compute_evaluator(syndromes, locator, count, evaluator);
compute_magnitudes(locator, locations, count, evaluator, evaluator_degree, magnitudes);
return evaluator_degree;
}
};
template <int NR, typename GF>
struct BerlekampMassey
{
typedef typename GF::value_type value_type;
typedef typename GF::ValueType ValueType;
typedef typename GF::IndexType IndexType;
static const int N = GF::N, K = N - NR;
static int algorithm(ValueType *s, ValueType *C, int count = 0)
{
ValueType B[NR+1];
for (int i = 0; i <= NR; ++i)
B[i] = C[i];
int L = count;
for (int n = count, m = 1; n < NR; ++n) {
ValueType d(s[n]);
for (int i = 1; i <= L; ++i)
d += C[i] * s[n-i];
if (!d) {
++m;
} else {
ValueType T[NR+1];
for (int i = 0; i < m; ++i)
T[i] = C[i];
for (int i = m; i <= NR; ++i)
T[i] = fma(d, B[i-m], C[i]);
if (2 * L <= n + count) {
L = n + count + 1 - L;
for (int i = 0; i <= NR; ++i)
B[i] = C[i] / d;
m = 1;
} else {
++m;
}
for (int i = 0; i <= NR; ++i)
C[i] = T[i];
}
}
return L;
}
};
}
template <int NR, int FCR, typename GF>
struct ReedSolomonErrorCorrection
{
typedef typename GF::value_type value_type;
typedef typename GF::ValueType ValueType;
typedef typename GF::IndexType IndexType;
static const int N = GF::N, K = N - NR;
static int algorithm(ValueType *syndromes, IndexType *locations, ValueType *magnitudes, IndexType *erasures = 0, int erasures_count = 0)
{
assert(0 <= erasures_count && erasures_count <= NR);
ValueType locator[NR+1];
locator[0] = ValueType(1);
for (int i = 1; i <= NR; ++i)
locator[i] = ValueType(0);
// $locator = \prod_{i=0}^{count}(1-x\,pe^{N-1-erasures_i})$
if (erasures_count)
locator[1] = value(IndexType(N-1) / erasures[0]);
for (int i = 1; i < erasures_count; ++i) {
IndexType tmp(IndexType(N-1) / erasures[i]);
for (int j = i; j >= 0; --j)
locator[j+1] += tmp * locator[j];
}
int locator_degree = RS::BerlekampMassey<NR, GF>::algorithm(syndromes, locator, erasures_count);
assert(locator_degree);
assert(locator_degree <= NR);
assert(locator[0] == ValueType(1));
while (!locator[locator_degree])
if (--locator_degree < 0)
return -1;
int count = RS::FindLocations<NR, GF>::search(locator, locator_degree, locations);
if (count < locator_degree)
return -1;
ValueType evaluator[NR];
int evaluator_degree = RS::Forney<NR, FCR, GF>::algorithm(syndromes, locator, locations, count, evaluator, magnitudes);
if (0) {
static bool once;
if (!once) {
once = true;
std::cout << "syndromes =";
for (int i = 0; i < NR; ++i)
std::cout << " " << (int)syndromes[i];
std::cout << std::endl;
std::cout << "locator = ";
for (int i = NR; i > 0; --i) {
if (!locator[i])
continue;
if (locator[i] != ValueType(1))
std::cout << (int)locator[i] << "*";
std::cout << "x";
if (i != 1)
std::cout << "^" << i;
std::cout << " + ";
}
std::cout << (int)locator[0] << std::endl;
std::cout << "locations =";
for (int i = 0; i < count; ++i)
std::cout << " " << (int)locations[i];
std::cout << std::endl;
std::cout << "evaluator = ";
for (int i = evaluator_degree; i > 0; --i) {
if (!evaluator[i])
continue;
if (evaluator[i] != ValueType(1))
std::cout << (int)evaluator[i] << "*";
std::cout << "x";
if (i != 1)
std::cout << "^" << i;
if (i != 1 || evaluator[0])
std::cout << " + ";
}
if (evaluator[0])
std::cout << (int)evaluator[0];
std::cout << std::endl;
std::cout << "magnitudes =";
for (int i = 0; i < count; ++i)
std::cout << " " << (int)magnitudes[i];
std::cout << std::endl;
}
}
return count;
}
};
}
#endif

65
tests/bch_decoder_test.cc Normal file
View file

@ -0,0 +1,65 @@
/*
Test for the Bose Chaudhuri Hocquenghem Decoder
Copyright 2018 Ahmet Inan <inan@aicodix.de>
*/
#include <cassert>
#include <random>
#include <iostream>
#include "galois_field.hh"
#include "bose_chaudhuri_hocquenghem_decoder.hh"
int main()
{
std::random_device rd;
std::default_random_engine generator(rd());
if (1) {
// NASA INTRO BCH(15, 5) T=3
typedef CODE::GaloisField<4, 0b10011, uint8_t> GF;
typedef CODE::BoseChaudhuriHocquenghemDecoder<6, 1, 5, GF> BCH;
GF instance;
BCH decode;
uint8_t target[15] = { 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0 };
uint8_t code[15];
for (int i = 0; i < 15; ++i)
code[i] = target[i];
std::uniform_int_distribution<uint8_t> distribution(0, 15);
auto noise = std::bind(distribution, generator);
for (int i = 0; i < 3; ++i)
code[noise()] ^= 1;
decode(code);
for (int i = 0; i < 15; ++i)
assert(code[i] == target[i]);
}
if (1) {
// DVB-S2 FULL BCH(65535, 65343) T=12
typedef CODE::GaloisField<16, 0b10000000000101101, uint16_t> GF;
typedef CODE::BoseChaudhuriHocquenghemDecoder<24, 1, 65343, GF> BCH;
GF *instance = new GF();
BCH *decode = new BCH();
uint16_t *target = new uint16_t[65535];
for (int i = 0, s = 0; i < 65343; ++i, s=(s*(s*s*51767+71287)+35149)&0xffffff)
target[i] = (s^=s>>7)&1;
uint16_t parity[192] = { 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0 };
for (int i = 0; i < 192; ++i)
target[65343+i] = parity[i];
uint16_t *code = new uint16_t[65535];
for (int i = 0; i < 65535; ++i)
code[i] = target[i];
std::uniform_int_distribution<uint16_t> distribution(0, 65535);
auto noise = std::bind(distribution, generator);
for (int i = 0; i < 12; ++i)
code[noise()] ^= 1;
(*decode)(code);
for (int i = 0; i < 65535; ++i)
assert(code[i] == target[i]);
delete[] target;
delete[] code;
delete decode;
delete instance;
}
std::cerr << "Bose Chaudhuri Hocquenghem Decoder test passed!" << std::endl;
return 0;
}

50
tests/bch_encoder_test.cc Normal file
View file

@ -0,0 +1,50 @@
/*
Test for the Bose Chaudhuri Hocquenghem Encoder
Copyright 2018 Ahmet Inan <inan@aicodix.de>
*/
#include <cassert>
#include <iostream>
#include "galois_field.hh"
#include "bose_chaudhuri_hocquenghem_encoder.hh"
int main()
{
if (1) {
// NASA INTRO BCH(15, 5) T=3
typedef CODE::GaloisField<4, 0b10011, uint8_t> GF;
typedef CODE::BoseChaudhuriHocquenghemEncoder<6, 1, 5, GF> BCH;
GF instance;
BCH encode({0b10011, 0b11111, 0b00111});
uint8_t code[15] = { 1, 1, 0, 0, 1 };
uint8_t target[15] = { 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0 };
encode(code);
for (int i = 0; i < 15; ++i)
assert(code[i] == target[i]);
}
if (1) {
// DVB-S2 FULL BCH(65535, 65343) T=12
typedef CODE::GaloisField<16, 0b10000000000101101, uint16_t> GF;
typedef CODE::BoseChaudhuriHocquenghemEncoder<24, 1, 65343, GF> BCH;
GF *instance = new GF();
BCH *encode = new BCH({0b10000000000101101, 0b10000000101110011, 0b10000111110111101, 0b10101101001010101, 0b10001111100101111, 0b11111011110110101, 0b11010111101100101, 0b10111001101100111, 0b10000111010100001, 0b10111010110100111, 0b10011101000101101, 0b10001101011100011});
uint16_t *code = new uint16_t[65535];
uint16_t *target = new uint16_t[65535];
for (int i = 0, s = 0; i < 65343; ++i, s=(s*(s*s*51767+71287)+35149)&0xffffff)
target[i] = code[i] = (s^=s>>7)&1;
uint16_t parity[192] = { 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0 };
for (int i = 0; i < 192; ++i)
target[65343+i] = parity[i];
(*encode)(code);
for (int i = 0; i < 65535; ++i)
assert(code[i] == target[i]);
delete[] target;
delete[] code;
delete encode;
delete instance;
}
std::cerr << "Bose Chaudhuri Hocquenghem Encoder test passed!" << std::endl;
return 0;
}

88
tests/rs_decoder_test.cc Normal file
View file

@ -0,0 +1,88 @@
/*
Test for the Reed Solomon Decoder
Copyright 2018 Ahmet Inan <inan@aicodix.de>
*/
#include <cassert>
#include <random>
#include <iostream>
#include "galois_field.hh"
#include "reed_solomon_decoder.hh"
int main()
{
std::random_device rd;
std::default_random_engine generator(rd());
if (1) {
// BBC WHP031 RS(15, 11) T=2
typedef CODE::GaloisField<4, 0b10011, uint8_t> GF;
typedef CODE::ReedSolomonDecoder<4, 0, GF> RS;
GF instance;
RS decode;
uint8_t target[15] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 3, 3, 12, 12 };
uint8_t code[15];
for (int i = 0; i < 15; ++i)
code[i] = target[i];
std::uniform_int_distribution<uint8_t> distribution(0, 15);
auto noise = std::bind(distribution, generator);
for (int i = 0; i < 2; ++i)
code[noise()] = noise();
decode(code);
for (int i = 0; i < 15; ++i)
assert(code[i] == target[i]);
}
if (1) {
// DVB-T RS(255, 239) T=8
typedef CODE::GaloisField<8, 0b100011101, uint8_t> GF;
typedef CODE::ReedSolomonDecoder<16, 0, GF> RS;
GF instance;
RS decode;
uint8_t target[255];
for (int i = 0; i < 239; ++i)
target[i] = i + 1;
uint8_t parity[16] = { 1, 126, 147, 48, 155, 224, 3, 157, 29, 226, 40, 114, 61, 30, 244, 75 };
for (int i = 0; i < 16; ++i)
target[239+i] = parity[i];
uint8_t code[255];
for (int i = 0; i < 255; ++i)
code[i] = target[i];
std::uniform_int_distribution<uint8_t> distribution(0, 255);
auto noise = std::bind(distribution, generator);
for (int i = 0; i < 8; ++i)
code[noise()] = noise();
decode(code);
for (int i = 0; i < 255; ++i)
assert(code[i] == target[i]);
}
if (1) {
// FUN RS(65535, 65471) T=32
typedef CODE::GaloisField<16, 0b10001000000001011, uint16_t> GF;
typedef CODE::ReedSolomonDecoder<64, 1, GF> RS;
GF *instance = new GF();
RS *decode = new RS();
uint16_t *target = new uint16_t[65535];
for (int i = 0; i < 65471; ++i)
target[i] = i + 1;
uint16_t parity[64] = { 25271, 26303, 22052, 31318, 31233, 6076, 40148, 29468, 47507, 32655, 12404, 13265, 23901, 38403, 50967, 50433, 40818, 226, 62296, 23636, 56393, 12952, 11476, 44416, 518, 50014, 10037, 57582, 33421, 42654, 54025, 7157, 4826, 52148, 17167, 23294, 6427, 40953, 11168, 35305, 18209, 1868, 39971, 54928, 27566, 1424, 4846, 25347, 34710, 42190, 56452, 21859, 49805, 28028, 41657, 25756, 22014, 24479, 28758, 17438, 12976, 61743, 46735, 1557 };
for (int i = 0; i < 64; ++i)
target[65471+i] = parity[i];
uint16_t *code = new uint16_t[65535];
for (int i = 0; i < 65535; ++i)
code[i] = target[i];
std::uniform_int_distribution<uint16_t> distribution(0, 65535);
auto noise = std::bind(distribution, generator);
for (int i = 0; i < 32; ++i)
code[noise()] = noise();
(*decode)(code);
for (int i = 0; i < 65535; ++i)
assert(code[i] == target[i]);
delete[] target;
delete[] code;
delete decode;
delete instance;
}
std::cerr << "Reed Solomon Decoder test passed!" << std::endl;
return 0;
}

66
tests/rs_encoder_test.cc Normal file
View file

@ -0,0 +1,66 @@
/*
Test for the Reed Solomon Encoder
Copyright 2018 Ahmet Inan <inan@aicodix.de>
*/
#include <cassert>
#include <iostream>
#include "galois_field.hh"
#include "reed_solomon_encoder.hh"
int main()
{
if (1) {
// BBC WHP031 RS(15, 11) T=2
typedef CODE::GaloisField<4, 0b10011, uint8_t> GF;
typedef CODE::ReedSolomonEncoder<4, 0, GF> RS;
GF instance;
RS encode;
uint8_t code[15] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
uint8_t target[15] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 3, 3, 12, 12 };
encode(code);
for (int i = 0; i < 15; ++i)
assert(code[i] == target[i]);
}
if (1) {
// DVB-T RS(255, 239) T=8
typedef CODE::GaloisField<8, 0b100011101, uint8_t> GF;
typedef CODE::ReedSolomonEncoder<16, 0, GF> RS;
GF instance;
RS encode;
uint8_t code[255], target[255];
for (int i = 0; i < 239; ++i)
target[i] = code[i] = i + 1;
uint8_t parity[16] = { 1, 126, 147, 48, 155, 224, 3, 157, 29, 226, 40, 114, 61, 30, 244, 75 };
for (int i = 0; i < 16; ++i)
target[239+i] = parity[i];
encode(code);
for (int i = 0; i < 255; ++i)
assert(code[i] == target[i]);
}
if (1) {
// FUN RS(65535, 65471) T=32
typedef CODE::GaloisField<16, 0b10001000000001011, uint16_t> GF;
typedef CODE::ReedSolomonEncoder<64, 1, GF> RS;
GF *instance = new GF();
RS *encode = new RS();
uint16_t *code = new uint16_t[65535];
uint16_t *target = new uint16_t[65535];
for (int i = 0; i < 65471; ++i)
target[i] = code[i] = i + 1;
uint16_t parity[64] = { 25271, 26303, 22052, 31318, 31233, 6076, 40148, 29468, 47507, 32655, 12404, 13265, 23901, 38403, 50967, 50433, 40818, 226, 62296, 23636, 56393, 12952, 11476, 44416, 518, 50014, 10037, 57582, 33421, 42654, 54025, 7157, 4826, 52148, 17167, 23294, 6427, 40953, 11168, 35305, 18209, 1868, 39971, 54928, 27566, 1424, 4846, 25347, 34710, 42190, 56452, 21859, 49805, 28028, 41657, 25756, 22014, 24479, 28758, 17438, 12976, 61743, 46735, 1557 };
for (int i = 0; i < 64; ++i)
target[65471+i] = parity[i];
(*encode)(code);
for (int i = 0; i < 65535; ++i)
assert(code[i] == target[i]);
delete[] target;
delete[] code;
delete encode;
delete instance;
}
std::cerr << "Reed Solomon Encoder test passed!" << std::endl;
return 0;
}