From 33052aae8df1a0bff5422067da6d041c3d038841 Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Sat, 3 Feb 2024 10:29:54 +0100 Subject: [PATCH] experimenting with list decoding --- osd.hh | 153 +++++++++++++++++++++++++++ tests/osld_regression_test.cc | 190 ++++++++++++++++++++++++++++++++++ 2 files changed, 343 insertions(+) create mode 100644 tests/osld_regression_test.cc diff --git a/osd.hh b/osd.hh index 9548cf2..91b9dcf 100644 --- a/osd.hh +++ b/osd.hh @@ -237,5 +237,158 @@ public: } }; +template +class OrderedStatisticsListDecoder +{ + static const int S = sizeof(size_t); + static const int W = (N+S-1) & ~(S-1); + int8_t G[W*K]; + int8_t codeword[W], candidate[W*2*L]; + int8_t softperm[W]; + int16_t perm[W]; + int score[2*L], scoreperm[2*L]; + void row_echelon() + { + for (int k = 0; k < K; ++k) { + // find pivot in this column + for (int j = k; j < K; ++j) { + if (G[W*j+k]) { + for (int i = k; j != k && i < N; ++i) + std::swap(G[W*j+i], G[W*k+i]); + break; + } + } + // keep searching for suitable column for pivot + // beware: this will use columns >= K if necessary. + for (int j = k + 1; !G[W*k+k] && j < N; ++j) { + for (int h = k; h < K; ++h) { + if (G[W*h+j]) { + // account column swap + std::swap(perm[k], perm[j]); + for (int i = 0; i < K; ++i) + std::swap(G[W*i+k], G[W*i+j]); + for (int i = k; h != k && i < N; ++i) + std::swap(G[W*h+i], G[W*k+i]); + break; + } + } + } + assert(G[W*k+k]); + // zero out column entries below pivot + for (int j = k + 1; j < K; ++j) + if (G[W*j+k]) + for (int i = k; i < N; ++i) + G[W*j+i] ^= G[W*k+i]; + } + } + void systematic() + { + for (int k = K-1; k; --k) + for (int j = 0; j < k; ++j) + if (G[W*j+k]) + for (int i = k; i < N; ++i) + G[W*j+i] ^= G[W*k+i]; + } + void encode() + { + for (int i = K; i < N; ++i) + codeword[i] = codeword[0] & G[i]; + for (int j = 1; j < K; ++j) + for (int i = K; i < N; ++i) + codeword[i] ^= codeword[j] & G[W*j+i]; + } + void flip(int j) + { + for (int i = 0; i < W; ++i) + codeword[i] ^= G[W*j+i]; + } + static int metric(const int8_t *hard, const int8_t *soft) + { + int sum = 0; + for (int i = 0; i < W; ++i) + sum += (1 - 2 * hard[i]) * soft[i]; + return sum; + } +public: + void operator()(int *rank, uint8_t *hard, const int8_t *soft, const int8_t *genmat) + { + for (int i = 0; i < N; ++i) + perm[i] = i; + for (int i = 0; i < N; ++i) + softperm[i] = std::abs(std::max(soft[i], -127)); + std::sort(perm, perm+N, [this](int a, int b){ return softperm[a] > softperm[b]; }); + for (int j = 0; j < K; ++j) + for (int i = 0; i < N; ++i) + G[W*j+i] = genmat[N*j+perm[i]]; + row_echelon(); + systematic(); + for (int i = 0; i < N; ++i) + softperm[i] = std::max(soft[perm[i]], -127); + for (int i = N; i < W; ++i) + softperm[i] = 0; + for (int i = 0; i < K; ++i) + codeword[i] = softperm[i] < 0; + encode(); + for (int i = 0; i < W; ++i) + candidate[i] = codeword[i]; + score[0] = metric(codeword, softperm); + int count = 1; + int worst = -1; + for (int i = 0; i < 2*L; ++i) + scoreperm[i] = i; + auto update = [this, &count, &worst]() { + int met = metric(codeword, softperm); + if (met > worst) { + score[scoreperm[count]] = met; + for (int i = 0; i < W; ++i) + candidate[scoreperm[count]*W+i] = codeword[i]; + if (++count >= 2*L) { + std::nth_element(scoreperm, scoreperm+(L-1), scoreperm+2*L, [this](int a, int b){ return score[a] > score[b]; }); + worst = score[scoreperm[L-1]]; + count = L; + } + } + }; + for (int a = 0; O >= 1 && a < K; ++a) { + flip(a); + update(); + for (int b = a + 1; O >= 2 && b < K; ++b) { + flip(b); + update(); + for (int c = b + 1; O >= 3 && c < K; ++c) { + flip(c); + update(); + for (int d = c + 1; O >= 4 && d < K; ++d) { + flip(d); + update(); + for (int e = d + 1; O >= 5 && e < K; ++e) { + flip(e); + update(); + for (int f = e + 1; O >= 6 && f < K; ++f) { + flip(f); + update(); + flip(f); + } + flip(e); + } + flip(d); + } + flip(c); + } + flip(b); + } + flip(a); + } + std::sort(scoreperm, scoreperm+count, [this](int a, int b){ return score[a] > score[b]; }); + for (int j = 0, r = 0; j < L; ++j) { + if (j > 0 && score[scoreperm[j-1]] != score[scoreperm[j]]) + ++r; + rank[j] = r; + for (int i = 0; i < N; ++i) + set_be_bit(hard+j*((N+7)/8), perm[i], candidate[scoreperm[j]*W+i]); + } + } +}; + } diff --git a/tests/osld_regression_test.cc b/tests/osld_regression_test.cc new file mode 100644 index 0000000..aed2cd5 --- /dev/null +++ b/tests/osld_regression_test.cc @@ -0,0 +1,190 @@ +/* +Regression Test for ordered statistics list decoding + +Copyright 2024 Ahmet Inan +*/ + +#include +#include +#include +#include +#include +#include "osd.hh" +#include "galois_field.hh" +#include "bose_chaudhuri_hocquenghem_encoder.hh" +#include "bose_chaudhuri_hocquenghem_decoder.hh" + +int main() +{ + const int L = 16; + const bool oracle = true; +#if 1 + // BCH(127, 64) T=10 + const int O = 2; + const int N = 127; + const int K = 64; + const int NR = 20; + const int loops = 1000; + const double low_SNR = -5; + const double high_SNR = 5; + const double QEF_SNR = 1.5; + typedef CODE::GaloisField<7, 0b10001001, uint8_t> GF; + std::initializer_list minpols { + 0b10001001, 0b10001111, 0b10011101, + 0b11110111, 0b10111111, 0b11010101, + 0b10000011, 0b11101111, 0b11001011, + }; +#endif +#if 0 + // NASA INTRO BCH(15, 5) T=3 + const int O = 1; + const int N = 15; + const int K = 5; + const int NR = 6; + const int loops = 100000; + const double low_SNR = -7; + const double high_SNR = 7; + const double QEF_SNR = 3.5; + typedef CODE::GaloisField<4, 0b10011, uint8_t> GF; + std::initializer_list minpols { + 0b10011, 0b11111, 0b00111 + }; +#endif + GF instance; + const int NW = (N+7)/8; + const int KW = (K+7)/8; + const int PW = (N-K+7)/8; + int8_t genmat[N*K]; + CODE::BoseChaudhuriHocquenghemGenerator::matrix(genmat, true, minpols); + CODE::BoseChaudhuriHocquenghemEncoder bchenc(minpols); + CODE::LinearEncoder linenc; + CODE::OrderedStatisticsListDecoder osddec; + CODE::BoseChaudhuriHocquenghemDecoder bchdec; + GF::value_type erasures[NR]; + uint8_t message[KW], parity[PW], decoded[NW*L], codeword[NW]; + int8_t orig[N], noisy[N]; + + std::random_device rd; + typedef std::default_random_engine generator; + typedef std::uniform_int_distribution distribution; + auto data = std::bind(distribution(0, 255), generator(rd())); + + int rank[L]; + double symb[N]; + double min_SNR = high_SNR; + for (double SNR = low_SNR; SNR <= high_SNR; SNR += 0.1) { + //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 normal; + auto awgn = std::bind(normal(mean_noise, sigma_noise), generator(rd())); + + int awgn_errors = 0; + int quantization_erasures = 0; + int uncorrected_errors = 0; + int ambiguity_erasures = 0; + int bchdec_errors = 0; + for (int l = 0; l < loops; ++l) { + for (int i = 0; i < KW; ++i) + message[i] = data(); + + linenc(codeword, message, genmat); + + for (int i = 0; i < K; ++i) + assert(CODE::get_be_bit(codeword, i) == CODE::get_be_bit(message, i)); + + bchenc(message, parity); + + for (int i = K; i < N; ++i) + assert(CODE::get_be_bit(codeword, i) == CODE::get_be_bit(parity, i-K)); + + for (int i = 0; i < N; ++i) + orig[i] = 1 - 2 * CODE::get_be_bit(codeword, i); + + for (int i = 0; i < N; ++i) + symb[i] = orig[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 = DIST / (sigma_noise * sigma_noise); + for (int i = 0; i < N; ++i) + noisy[i] = std::min(std::max(std::nearbyint(fact * symb[i]), -128), 127); + + osddec(rank, decoded, noisy, genmat); + bool unique = rank[0] != rank[1]; + + if (oracle) { + for (int j = 0; j < L; ++j) { + bool found = true; + for (int i = 0; i < N; ++i) + found &= CODE::get_be_bit(decoded+NW*j, i) == CODE::get_be_bit(codeword, i); + if (found) { + for (int i = 0; j && i < NW; ++i) + decoded[i] = decoded[NW*j+i]; + unique = true; + break; + } + } + } + + 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 < N; ++i) + uncorrected_errors += CODE::get_be_bit(decoded, i) != CODE::get_be_bit(codeword, i); + if (!unique) + ambiguity_erasures += N; + + for (int i = 0; i < K; ++i) + CODE::set_be_bit(message, i, noisy[i] < 0); + for (int i = K; i < N; ++i) + CODE::set_be_bit(parity, i-K, noisy[i] < 0); + int erasures_count = 0; + for (int i = 0; erasures_count < NR && i < N; ++i) + if (!noisy[i]) + erasures[erasures_count++] = i; + + bchdec(message, parity, erasures, erasures_count); + + for (int i = 0; i < K; ++i) + bchdec_errors += CODE::get_be_bit(message, i) != CODE::get_be_bit(codeword, i); + for (int i = K; i < N; ++i) + bchdec_errors += CODE::get_be_bit(parity, i-K) != CODE::get_be_bit(codeword, i); + } + + double bit_error_rate = (double)uncorrected_errors / (double)(N * loops); + if (!uncorrected_errors && !ambiguity_erasures) + min_SNR = std::min(min_SNR, SNR); + + 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)); + + double bchdec_ber = (double)bchdec_errors / (double)(N * loops); + + 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 << bchdec_ber << " BCH decoder bit error rate." << std::endl; + } else { + std::cout << SNR << " " << bit_error_rate << " " << bchdec_ber << " " << EbN0 << std::endl; + } + } + std::cerr << "QEF at: " << min_SNR << " SNR" << std::endl; + assert(min_SNR < QEF_SNR); + std::cerr << "Ordered statistics list decoding regression test passed!" << std::endl; + return 0; +}