From 01db8f516e9b71abe8c2af5b5fb2472b36cc235b Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Thu, 20 Sep 2018 18:04:23 +0200 Subject: [PATCH] added Reed Solomon error correction lifted from my FEC project: https://github.com/xdsopl/FEC --- README.md | 11 +- bose_chaudhuri_hocquenghem_decoder.hh | 77 ++++++++ bose_chaudhuri_hocquenghem_encoder.hh | 86 +++++++++ reed_solomon_decoder.hh | 74 ++++++++ reed_solomon_encoder.hh | 78 ++++++++ reed_solomon_error_correction.hh | 252 ++++++++++++++++++++++++++ tests/bch_decoder_test.cc | 65 +++++++ tests/bch_encoder_test.cc | 50 +++++ tests/rs_decoder_test.cc | 88 +++++++++ tests/rs_encoder_test.cc | 66 +++++++ 10 files changed, 846 insertions(+), 1 deletion(-) create mode 100644 bose_chaudhuri_hocquenghem_decoder.hh create mode 100644 bose_chaudhuri_hocquenghem_encoder.hh create mode 100644 reed_solomon_decoder.hh create mode 100644 reed_solomon_encoder.hh create mode 100644 reed_solomon_error_correction.hh create mode 100644 tests/bch_decoder_test.cc create mode 100644 tests/bch_encoder_test.cc create mode 100644 tests/rs_decoder_test.cc create mode 100644 tests/rs_encoder_test.cc diff --git a/README.md b/README.md index 64e0de4..3a8103e 100644 --- a/README.md +++ b/README.md @@ -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 [Reed–Solomon 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) diff --git a/bose_chaudhuri_hocquenghem_decoder.hh b/bose_chaudhuri_hocquenghem_decoder.hh new file mode 100644 index 0000000..9c8e9f7 --- /dev/null +++ b/bose_chaudhuri_hocquenghem_decoder.hh @@ -0,0 +1,77 @@ +/* +Bose Chaudhuri Hocquenghem Decoder + +Copyright 2018 Ahmet Inan +*/ + +#ifndef BOSE_CHAUDHURI_HOCQUENGHEM_DECODER_HH +#define BOSE_CHAUDHURI_HOCQUENGHEM_DECODER_HH + +#include "reed_solomon_error_correction.hh" + +namespace CODE { + +template +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::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(code), reinterpret_cast(erasures), erasures_count); + } + int compute_syndromes(value_type *code, value_type *syndromes) + { + return compute_syndromes(reinterpret_cast(code), reinterpret_cast(syndromes)); + } +}; + +} + +#endif diff --git a/bose_chaudhuri_hocquenghem_encoder.hh b/bose_chaudhuri_hocquenghem_encoder.hh new file mode 100644 index 0000000..7c71ee7 --- /dev/null +++ b/bose_chaudhuri_hocquenghem_encoder.hh @@ -0,0 +1,86 @@ +/* +Bose Chaudhuri Hocquenghem Encoder + +Copyright 2018 Ahmet Inan +*/ + +#ifndef BOSE_CHAUDHURI_HOCQUENGHEM_ENCODER_HH +#define BOSE_CHAUDHURI_HOCQUENGHEM_ENCODER_HH + +#include + +namespace CODE { + +template +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 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(code)); + } +}; + +} + +#endif diff --git a/reed_solomon_decoder.hh b/reed_solomon_decoder.hh new file mode 100644 index 0000000..ca8392b --- /dev/null +++ b/reed_solomon_decoder.hh @@ -0,0 +1,74 @@ +/* +Reed Solomon Decoder + +Copyright 2018 Ahmet Inan +*/ + +#ifndef REED_SOLOMON_DECODER_HH +#define REED_SOLOMON_DECODER_HH + +#include "reed_solomon_error_correction.hh" + +namespace CODE { + +template +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::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(code), reinterpret_cast(erasures), erasures_count); + } + int compute_syndromes(value_type *code, value_type *syndromes) + { + return compute_syndromes(reinterpret_cast(code), reinterpret_cast(syndromes)); + } +}; + +} + +#endif diff --git a/reed_solomon_encoder.hh b/reed_solomon_encoder.hh new file mode 100644 index 0000000..508dfdb --- /dev/null +++ b/reed_solomon_encoder.hh @@ -0,0 +1,78 @@ +/* +Reed Solomon Encoder + +Copyright 2018 Ahmet Inan +*/ + +#ifndef REED_SOLOMON_ENCODER_HH +#define REED_SOLOMON_ENCODER_HH + +namespace CODE { + +template +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(code)); + } +}; + +} + +#endif diff --git a/reed_solomon_error_correction.hh b/reed_solomon_error_correction.hh new file mode 100644 index 0000000..37b9e22 --- /dev/null +++ b/reed_solomon_error_correction.hh @@ -0,0 +1,252 @@ +/* +Reed Solomon Error Correction + +Copyright 2018 Ahmet Inan +*/ + +#ifndef REED_SOLOMON_ERROR_CORRECTION_HH +#define REED_SOLOMON_ERROR_CORRECTION_HH + +namespace CODE { +namespace RS { + +template +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 +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::search(locator, locator_degree, locations); + } +}; + +template +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 +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 +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::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::search(locator, locator_degree, locations); + if (count < locator_degree) + return -1; + ValueType evaluator[NR]; + int evaluator_degree = RS::Forney::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 diff --git a/tests/bch_decoder_test.cc b/tests/bch_decoder_test.cc new file mode 100644 index 0000000..a5bc9d8 --- /dev/null +++ b/tests/bch_decoder_test.cc @@ -0,0 +1,65 @@ +/* +Test for the Bose Chaudhuri Hocquenghem Decoder + +Copyright 2018 Ahmet Inan +*/ + +#include +#include +#include +#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 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 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; +} + diff --git a/tests/bch_encoder_test.cc b/tests/bch_encoder_test.cc new file mode 100644 index 0000000..5a70980 --- /dev/null +++ b/tests/bch_encoder_test.cc @@ -0,0 +1,50 @@ +/* +Test for the Bose Chaudhuri Hocquenghem Encoder + +Copyright 2018 Ahmet Inan +*/ + +#include +#include +#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; +} + diff --git a/tests/rs_decoder_test.cc b/tests/rs_decoder_test.cc new file mode 100644 index 0000000..be1ad16 --- /dev/null +++ b/tests/rs_decoder_test.cc @@ -0,0 +1,88 @@ +/* +Test for the Reed Solomon Decoder + +Copyright 2018 Ahmet Inan +*/ + +#include +#include +#include +#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 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 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 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; +} + diff --git a/tests/rs_encoder_test.cc b/tests/rs_encoder_test.cc new file mode 100644 index 0000000..f3194dd --- /dev/null +++ b/tests/rs_encoder_test.cc @@ -0,0 +1,66 @@ +/* +Test for the Reed Solomon Encoder + +Copyright 2018 Ahmet Inan +*/ + +#include +#include +#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; +} +