diff --git a/short_bch_code_decoder.hh b/short_bch_code_decoder.hh index 841bedc..e4df202 100644 --- a/short_bch_code_decoder.hh +++ b/short_bch_code_decoder.hh @@ -50,6 +50,84 @@ public: { return inp ^ err[(par[inp>>P] ^ inp) & (R-1)]; } + int operator()(const int8_t *code) + { + // maximum likelihood + if (0) { + int word = 0, best = 0; + for (int msg = 0; msg < W; ++msg) { + int sum = 0; + int enc = (msg << P) | par[msg]; + int lol = 8 * sizeof(enc) - 1; + for (int i = 0; i < N; ++i) + sum += ((enc << (lol-i)) >> lol) * code[i]; + if (sum > best) { + best = sum; + word = enc; + } + } + return word; + } + int cw = 0; + for (int i = 0; i < N; ++i) + cw |= (code[i] < 0) << i; + // hard decision + if (0) + return (*this)(cw); + int word = 0, best = 0; + // flip each bit and see .. + if (0) { + for (int j = 0; j < N; ++j) { + int tmp = 1 << j; + int dec = (*this)(cw ^ tmp); + int lol = 8 * sizeof(dec) - 1; + int sum = 0; + for (int i = 0; i < N; ++i) + sum += ((dec << (lol-i)) >> lol) * code[i]; + if (sum > best) { + best = sum; + word = dec; + } + } + } + // Chase algorithm + if (1) { + const int num = 4; + int worst[num] = { 0 }; + for (int i = 0; i < N; ++i) { + if (std::abs(code[i]) < std::abs(code[worst[0]])) { + worst[3] = worst[2]; + worst[2] = worst[1]; + worst[1] = worst[0]; + worst[0] = i; + } else if (std::abs(code[i]) < std::abs(code[worst[1]])) { + worst[3] = worst[2]; + worst[2] = worst[1]; + worst[1] = i; + } else if (std::abs(code[i]) < std::abs(code[worst[2]])) { + worst[3] = worst[2]; + worst[2] = i; + } else if (std::abs(code[i]) < std::abs(code[worst[3]])) { + worst[3] = i; + } + } + for (int j = 0; j < (1 << num); ++j) { + int tmp = 0; + for (int i = 0; i < num; ++i) + tmp |= ((j>>i)&1) << worst[i]; + int dec = (*this)(cw ^ tmp); + int lol = 8 * sizeof(dec) - 1; + int sum = 0; + for (int i = 0; i < N; ++i) + sum += ((dec << (lol-i)) >> lol) * code[i]; + if (sum > best) { + best = sum; + word = dec; + } + } + } + return word; + } }; } diff --git a/short_bch_code_encoder.hh b/short_bch_code_encoder.hh index b658a69..173f6e1 100644 --- a/short_bch_code_encoder.hh +++ b/short_bch_code_encoder.hh @@ -34,6 +34,12 @@ public: { return (msg << P) | par[msg]; } + void operator()(int8_t *code, int msg) + { + int cw = (*this)(msg); + for (int i = 0; i < N; ++i) + code[i] = 1 - 2 * ((cw >> i) & 1); + } }; } diff --git a/tests/soft_bch_regression_test.cc b/tests/soft_bch_regression_test.cc new file mode 100644 index 0000000..e6d2cdb --- /dev/null +++ b/tests/soft_bch_regression_test.cc @@ -0,0 +1,119 @@ +/* +Regression Test for the short BCH code encoder and soft decoder + +Copyright 2020 Ahmet Inan +*/ + +#include +#include +#include +#include +#include +#include +#include "short_bch_code_decoder.hh" +#include "short_bch_code_encoder.hh" + +template +int popcnt(TYPE x) +{ + int cnt = 0; + while (x) { + ++cnt; + x &= x-1; + } + return cnt; +} + +#if 1 + // Perfect binary Golay code using x^11+x^9+x^7+x^6+x^5+x+1 + const int LOOPS = 100000; + const float QEF_SNR = 3.0; + const int CODE_LEN = 23; + const int DATA_LEN = 12; + const int RADIUS_T = 3; + const int GEN_POLY = 0b101011100011; +#endif + +int main() +{ + CODE::ShortBCHCodeEncoder encode(GEN_POLY); + CODE::ShortBCHCodeDecoder decode(GEN_POLY, RADIUS_T); + + std::random_device rd; + std::default_random_engine generator(rd()); + typedef std::uniform_int_distribution uniform; + typedef std::normal_distribution normal; + + float min_SNR = 20; + + for (float SNR = -10; SNR <= 10; SNR += 0.1) { + //float mean_signal = 0; + float sigma_signal = 1; + float mean_noise = 0; + float sigma_noise = std::sqrt(sigma_signal * sigma_signal / (2 * std::pow(10, SNR / 10))); + + auto data = std::bind(uniform(0, (1 << DATA_LEN) - 1), generator); + auto awgn = std::bind(normal(mean_noise, sigma_noise), generator); + + int awgn_errors = 0; + int quantization_erasures = 0; + int uncorrected_errors = 0; + int decoder_errors = 0; + for (int loop = 0; loop < LOOPS; ++loop) { + int8_t code[CODE_LEN], orig[CODE_LEN], noisy[CODE_LEN]; + float symb[CODE_LEN]; + + int dat = data(); + encode(code, dat); + + for (int i = 0; i < CODE_LEN; ++i) + orig[i] = code[i]; + + for (int i = 0; i < CODE_LEN; ++i) + symb[i] = code[i]; + + for (int i = 0; i < CODE_LEN; ++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}}$ + float DIST = 2; // BPSK + float fact = DIST / (sigma_noise * sigma_noise); + for (int i = 0; i < CODE_LEN; ++i) + code[i] = std::min(std::max(std::nearbyint(fact * symb[i]), -128), 127); + + for (int i = 0; i < CODE_LEN; ++i) + noisy[i] = code[i]; + + int dec = decode(code) >> (CODE_LEN - DATA_LEN); + + for (int i = 0; i < CODE_LEN; ++i) + awgn_errors += noisy[i] * orig[i] < 0; + for (int i = 0; i < CODE_LEN; ++i) + quantization_erasures += !noisy[i]; + uncorrected_errors += popcnt(dat^dec); + for (int i = 0; i < DATA_LEN; ++i) + decoder_errors += ((dec^dat)&(1< 0; + } + float bit_error_rate = (float)uncorrected_errors / (float)(DATA_LEN * LOOPS); + if (bit_error_rate < 0.0001) + min_SNR = std::min(min_SNR, SNR); + + if (0) { + std::cerr << SNR << " Es/N0 => AWGN with standard deviation of " << sigma_noise << " and mean " << mean_noise << std::endl; + std::cerr << awgn_errors << " errors caused by AWGN." << std::endl; + std::cerr << quantization_erasures << " erasures caused by quantization." << std::endl; + std::cerr << decoder_errors << " errors caused by decoder." << std::endl; + std::cerr << uncorrected_errors << " errors uncorrected." << std::endl; + std::cerr << bit_error_rate << " bit error rate." << std::endl; + } else { + std::cout << SNR << " " << bit_error_rate << std::endl; + } + } + + std::cerr << "QEF at: " << min_SNR << " SNR" << std::endl; + assert(min_SNR < QEF_SNR); + std::cerr << "Soft BCH regression test passed!" << std::endl; + return 0; +} +