From a8d963cdfd5c6f182cf200fde15ddb9ea5bf1af8 Mon Sep 17 00:00:00 2001 From: Karlis Goba Date: Wed, 14 Nov 2018 17:06:32 +0200 Subject: [PATCH] Added belief propagation decoder as per WSJT-X --- decode_ft8.cpp | 36 +++--- ft8/arrays.h | 20 +++- ft8/ldpc.cpp | 300 ++++++++++++++++++++++++++++--------------------- ft8/ldpc.h | 4 +- 4 files changed, 214 insertions(+), 146 deletions(-) diff --git a/decode_ft8.cpp b/decode_ft8.cpp index 38edb08..41ee7ac 100644 --- a/decode_ft8.cpp +++ b/decode_ft8.cpp @@ -192,15 +192,15 @@ uint8_t max2(uint8_t a, uint8_t b) { } -uint8_t max4(uint8_t a, uint8_t b, uint8_t c, uint8_t d) { - return max2(max2(a, b), max2(c, d)); +uint8_t max4(uint8_t a, uint8_t b, uint8_t cand, uint8_t d) { + return max2(max2(a, b), max2(cand, d)); } -// Compute log likelihood log(p(0) / p(1)) of 174 message bits +// Compute log likelihood log(p(1) / p(0)) of 174 message bits // for later use in soft-decision LDPC decoding -void extract_likelihood(const uint8_t * power, int num_bins, const Candidate & c, float * log174) { - int offset = (c.time_offset * 4 + c.time_sub * 2 + c.freq_sub) * num_bins + c.freq_offset; +void extract_likelihood(const uint8_t * power, int num_bins, const Candidate & cand, float * log174) { + int offset = (cand.time_offset * 4 + cand.time_sub * 2 + cand.freq_sub) * num_bins + cand.freq_offset; int k = 0; // Go over FSK tones and skip Costas sync symbols @@ -212,9 +212,9 @@ void extract_likelihood(const uint8_t * power, int num_bins, const Candidate & c // Extract bit significance (and convert them to float) // 8 FSK tones = 3 bits - log174[k + 0] = -(int)max4(ps[4], ps[5], ps[6], ps[7]) + (int)max4(ps[0], ps[1], ps[2], ps[3]); - log174[k + 1] = -(int)max4(ps[2], ps[3], ps[6], ps[7]) + (int)max4(ps[0], ps[1], ps[4], ps[5]); - log174[k + 2] = -(int)max4(ps[1], ps[3], ps[5], ps[7]) + (int)max4(ps[0], ps[2], ps[4], ps[6]); + log174[k + 0] = (int)max4(ps[4], ps[5], ps[6], ps[7]) - (int)max4(ps[0], ps[1], ps[2], ps[3]); + log174[k + 1] = (int)max4(ps[2], ps[3], ps[6], ps[7]) - (int)max4(ps[0], ps[1], ps[4], ps[5]); + log174[k + 2] = (int)max4(ps[1], ps[3], ps[5], ps[7]) - (int)max4(ps[0], ps[2], ps[4], ps[6]); // printf("%d %d %d %d %d %d %d %d : %.0f %.0f %.0f\n", // ps[0], ps[1], ps[2], ps[3], ps[4], ps[5], ps[6], ps[7], // log174[k + 0], log174[k + 1], log174[k + 2]); @@ -270,29 +270,31 @@ int main(int argc, char ** argv) { extract_power(signal, num_blocks, num_bins, power); - const int num_candidates = 250; + int num_candidates = 250; Candidate heap[num_candidates]; find_sync(power, num_blocks, num_bins, num_candidates, heap); - for (int i = 0; i < num_candidates; ++i) { - float freq_hz = (heap[i].freq_offset + heap[i].freq_sub / 2.0f) * fsk_dev; - float time_sec = (heap[i].time_offset + heap[i].time_sub / 2.0f) / fsk_dev; + for (int idx = 0; idx < num_candidates; ++idx) { + Candidate &cand = heap[idx]; + float freq_hz = (cand.freq_offset + cand.freq_sub / 2.0f) * fsk_dev; + float time_sec = (cand.time_offset + cand.time_sub / 2.0f) / fsk_dev; // printf("%03d: score = %d freq = %.1f time = %.2f\n", i, // heap[i].score, freq_hz, time_sec); float log174[3 * ND]; - extract_likelihood(power, num_bins, heap[i], log174); + extract_likelihood(power, num_bins, cand, log174); - const int num_iters = 10; + const int num_iters = 20; int plain[3 * ND]; int ok; - ldpc_decode(log174, num_iters, plain, &ok); + bp_decode(log174, num_iters, plain, &ok); + //ldpc_decode(log174, num_iters, plain, &ok); //printf("ldpc_decode() = %d\n", ok); if (ok == 87) { - printf("%03d: score = %d freq = %.1f time = %.2f\n", i, - heap[i].score, freq_hz, time_sec); + printf("%03d: score = %d freq = %.1f time = %.2f\n", idx, + cand.score, freq_hz, time_sec); } } diff --git a/ft8/arrays.h b/ft8/arrays.h index 9e0d38f..fe7ebb6 100644 --- a/ft8/arrays.h +++ b/ft8/arrays.h @@ -1,7 +1,9 @@ +#include + // LDPC(174,87) parameters from WSJT-X. // this is an indirection table that moves a // codeword's 87 systematic (message) bits to the end. -int colorder[] = { +uint8_t colorder[] = { 0, 1, 2, 3, 30, 4, 5, 6, 7, 8, 9, 10, 11, 32, 12, 40, 13, 14, 15, 16, 17, 18, 37, 45, 29, 19, 20, 21, 41, 22, 42, 31, 33, 34, 44, 35, 47, 51, 50, 43, 36, 52, 63, 46, 25, 55, 27, 24, 23, 53, 39, 49, 59, 38, @@ -21,7 +23,7 @@ int colorder[] = { // each number is an index into the codeword (1-origin). // the codeword bits mentioned in each row must xor to zero. // From WSJT-X's bpdecode174.f90. -int Nm[][7] = { +uint8_t Nm[][7] = { {1, 30, 60, 89, 118, 147, 0}, {2, 31, 61, 90, 119, 147, 0}, {3, 32, 62, 91, 120, 148, 0}, @@ -116,7 +118,7 @@ int Nm[][7] = { // the numbers indicate which three parity // checks (rows in Nm) refer to the codeword bit. // 1-origin. -int Mn[][3] = { +uint8_t Mn[][3] = { {1, 25, 69}, {2, 5, 73}, {3, 32, 68}, @@ -292,3 +294,15 @@ int Mn[][3] = { {54, 61, 71}, {46, 58, 69}, }; + +uint8_t nrw[] = { + 6,6,6,6,6,6,6,6,6,6, + 6,6,6,6,6,6,6,6,6,6, + 6,6,6,6,6,6,6,6,6,6, + 6,6,6,6,6,6,6,6,6,7, + 6,6,6,6,6,7,6,6,6,6, + 6,6,6,6,6,7,6,6,6,6, + 7,6,5,6,6,6,6,6,6,5, + 6,6,6,6,6,6,6,6,6,6, + 5,6,6,6,5,6,6 +}; diff --git a/ft8/ldpc.cpp b/ft8/ldpc.cpp index 0775207..698c09e 100644 --- a/ft8/ldpc.cpp +++ b/ft8/ldpc.cpp @@ -16,6 +16,13 @@ int ldpc_check(int codeword[]); +const int N = 174; +const int M = 87; + +// https://varietyofsound.wordpress.com/2011/02/14/efficient-tanh-computation-using-lamberts-continued-fraction/ +// http://functions.wolfram.com/ElementaryFunctions/ArcTanh/10/0001/ +// https://mathr.co.uk/blog/2017-09-06_approximating_hyperbolic_tangent.html + // thank you Douglas Bagnall // https://math.stackexchange.com/a/446411 @@ -27,69 +34,127 @@ float fast_tanh(float x) { return 1.0f; } float x2 = x * x; - float a = x * (135135.0f + x2 * (17325.0f + x2 * (378.0f + x2))); - float b = 135135.0f + x2 * (62370.0f + x2 * (3150.0f + x2 * 28.0f)); + //float a = x * (135135.0f + x2 * (17325.0f + x2 * (378.0f + x2))); + //float b = 135135.0f + x2 * (62370.0f + x2 * (3150.0f + x2 * 28.0f)); + //float a = x * (10395.0f + x2 * (1260.0f + x2 * 21.0f)); + //float b = 10395.0f + x2 * (4725.0f + x2 * (210.0f + x2)); + float a = x * (945.0f + x2 * (105.0f + x2)); + float b = 945.0f + x2 * (420.0f + x2 * 15.0f); return a / b; } +float fast_atanh(float x) { + float x2 = x * x; + //float a = x * (-15015.0f + x2 * (19250.0f + x2 * (-5943.0f + x2 * 256.0f))); + //float b = (-15015.0f + x2 * (24255.0f + x2 * (-11025.0f + x2 * 1225.0f))); + //float a = x * (-1155.0f + x2 * (1190.0f + x2 * -231.0f)); + //float b = (-1155.0f + x2 * (1575.0f + x2 * (-525.0f + x2 * 25.0f))); + float a = x * (945.0f + x2 * (-735.0f + x2 * 64.0f)); + float b = (945.0f + x2 * (-1050.0f + x2 * 225.0f)); + return a / b; +} + + +float pltanh(float x) { + float isign = +1; + if (x < 0) { + isign = -1; + x = -x; + } + if (x < 0.8f) { + return isign * 0.83 * x; + } + if (x < 1.6f) { + return isign * (0.322f * x + 0.4064f); + } + if (x < 3.0f) { + return isign * (0.0524f * x + 0.8378f); + } + if (x < 7.0f) { + return isign * (0.0012f * x + 0.9914f); + } + return isign*0.9998f; +} + + +float platanh(float x) { + float isign = +1; + if (x < 0) { + isign = -1; + x = -x; + } + if (x < 0.664f) { + return isign * x / 0.83f; + } + if (x < 0.9217f) { + return isign * (x - 0.4064f) / 0.322f; + } + if (x < 0.9951f) { + return isign * (x - 0.8378f) / 0.0524f; + } + if (x < 0.9998f) { + return isign * (x - 0.9914f) / 0.0012f; + } + return isign * 7.0f; +} + + // codeword is 174 log-likelihoods. // plain is a return value, 174 ints, to be 0 or 1. -// iters is how hard to try. +// max_iters is how hard to try. // ok == 87 means success. -void ldpc_decode(float codeword[], int iters, int plain[], int *ok) { - float m[87][174]; // ~60 kB - float e[87][174]; // ~60 kB +void ldpc_decode(float codeword[], int max_iters, int plain[], int *ok) { + float m[M][N]; // ~60 kB + float e[M][N]; // ~60 kB int best_score = -1; - int best_cw[174]; - for (int i = 0; i < 174; i++) { - for (int j = 0; j < 87; j++) { + for (int j = 0; j < M; j++) { + for (int i = 0; i < N; i++) { m[j][i] = codeword[i]; e[j][i] = 0.0f; } } - for (int iter = 0; iter < iters; iter++) { - for (int j = 0; j < 87; j++) { - for (int ii1 = 0; ii1 < 7; ii1++) { + for (int iter = 0; iter < max_iters; iter++) { + for (int j = 0; j < M; j++) { + for (int ii1 = 0; ii1 < nrw[j]; ii1++) { int i1 = Nm[j][ii1] - 1; - if (i1 < 0) { - continue; - } float a = 1.0f; - for (int ii2 = 0; ii2 < 7; ii2++) { + for (int ii2 = 0; ii2 < nrw[j]; ii2++) { int i2 = Nm[j][ii2] - 1; - if (i2 >= 0 && i2 != i1) { - a *= fast_tanh(m[j][i2] / 2.0f); + if (i2 != i1) { + a *= fast_tanh(-m[j][i2] / 2.0f); } } - e[j][i1] = logf((1 + a) / (1 - a)); + e[j][i1] = logf((1 - a) / (1 + a)); } } - int cw[174]; - for (int i = 0; i < 174; i++) { + int cw[N]; + for (int i = 0; i < N; i++) { float l = codeword[i]; for (int j = 0; j < 3; j++) l += e[Mn[i][j] - 1][i]; - cw[i] = (l <= 0.0f); + cw[i] = (l > 0) ? 1 : 0; } + int score = ldpc_check(cw); if (score > best_score) { - for (int i = 0; i < 174; i++) { - best_cw[i] = cw[i]; + for (int i = 0; i < N; i++) { + plain[i] = cw[colorder[i]]; } best_score = score; + + if (score == M) { + // Found a perfect answer + break; + } } - if (score == 87) { - // Found a perfect answer - break; - } - for (int i = 0; i < 174; i++) { + for (int i = 0; i < N; i++) { for (int ji1 = 0; ji1 < 3; ji1++) { int j1 = Mn[i][ji1] - 1; float l = codeword[i]; @@ -104,17 +169,6 @@ void ldpc_decode(float codeword[], int iters, int plain[], int *ok) { } } - // decode complete (perhaps partially) -#if 0 - for(int i = 0; i < 87; i++) { - plain[i] = best_cw[colorder[174-87+i]]; - } -#else - for (int i = 0; i < 174; i++) { - plain[i] = best_cw[colorder[i]]; - } -#endif - *ok = best_score; } @@ -128,108 +182,104 @@ int ldpc_check(int codeword[]) { int score = 0; // Nm[87][7] - for (int j = 0; j < 87; j++) { + for (int j = 0; j < M; ++j) { int x = 0; - for (int ii1 = 0; ii1 < 7; ii1++) { - int i1 = Nm[j][ii1] - 1; - if (i1 >= 0) { - x ^= codeword[i1]; - } + for (int ii1 = 0; ii1 < nrw[j]; ++ii1) { + x ^= codeword[Nm[j][ii1] - 1]; } if (x == 0) { - score++; + ++score; } } return score; } -/* -def bp_decode(codeword, max_iterations = 10): - ## 174 codeword bits - ## 87 parity checks +void bp_decode(float codeword[], int max_iters, int plain[], int *ok) { + float tov[N][3]; + float toc[M][7]; - mnx = numpy.array(Mn, dtype=numpy.int32) - nmx = numpy.array(Nm, dtype=numpy.int32) + int best_score = -1; - ncw = 3 - tov = numpy.zeros( (3, N) ) - toc = numpy.zeros( (7, M) ) - tanhtoc = numpy.zeros( (7, M) ) - zn = numpy.zeros(N) - nclast = 0 - ncnt = 0 + int nclast = 0; + int ncnt = 0; - # initialize messages to checks - for j in range(M): - for i in range(nrw[j]): - toc[i, j] = codeword[nmx[j, i] - 1] - - for iteration in range(max_iterations): - # Update bit log likelihood ratios (tov=0 in iteration 0). - #for i in range(N): - # zn[i] = codeword[i] + numpy.sum(tov[:,i]) - zn = codeword + numpy.sum(tov, axis = 0) - #print(numpy.sum(tov, axis=0)) + // initialize messages to checks + for (int i = 0; i < M; ++i) { + for (int j = 0; j < nrw[i]; ++j) { + toc[i][j] = codeword[Nm[i][j] - 1]; + } + } - # Check to see if we have a codeword (check before we do any iteration). - cw = numpy.zeros(N, dtype=numpy.int32) - cw[zn > 0] = 1 - ncheck = 0 - for i in range(M): - synd = numpy.sum(cw[ nmx[i, :nrw[i]]-1 ]) - if synd % 2 > 0: - ncheck += 1 + for (int i = 0; i < N; ++i) { + for (int j = 0; j < 3; ++j) { + tov[i][j] = 0; + } + } - if ncheck == 0: - # we have a codeword - reorder the columns and return it - codeword = cw[colorder] - #nerr = 0 - #for i in range(N): - # if (2*cw[i]-1)*codeword[i] < 0: - # nerr += 1 - - #print("DECODED!", nerr) - return codeword[M:N] + for (int iter = 0; iter < max_iters; ++iter) { + float zn[N]; + int cw[N]; - if iter > 0: - # this code block implements an early stopping criterion - nd = ncheck - nclast - if nd < 0: # of unsatisfied parity checks decreased - ncnt = 0 # reset counter - else: - ncnt += 1 - - if ncnt >= 5 and iter >= 10 and ncheck >= 15: - nharderror = -1 - #return numpy.array([]) + // Update bit log likelihood ratios (tov=0 in iter 0) + for (int i = 0; i < N; ++i) { + zn[i] = codeword[i] + tov[i][0] + tov[i][1] + tov[i][2]; + cw[i] = (zn[i] > 0) ? 1 : 0; + } - nclast = ncheck + // Check to see if we have a codeword (check before we do any iter) + int score = ldpc_check(cw); - # Send messages from bits to check nodes - for j in range(M): - for i in range(nrw[j]): - ibj = nmx[j, i] - 1 - toc[i, j] = zn[ibj] - for kk in range(ncw): # subtract off what the bit had received from the check - if mnx[ibj, kk] - 1 == j: - toc[i, j] -= tov[kk, ibj] + if (score > best_score) { + // we have a better guess - reorder the columns and store it + for (int i = 0; i < N; i++) { + plain[i] = cw[colorder[i]]; + } + best_score = score; - # send messages from check nodes to variable nodes - #for i in range(M): - # tanhtoc[:,i] = numpy.tanh(-toc[:,i] / 2) - tanhtoc = numpy.tanh(-toc / 2) + if (score == M) { + break; + } + } - for j in range(N): - for i in range(ncw): - ichk = mnx[j, i] - 1 # Mn(:,j) are the checks that include bit j - Tmn = 1.0 - for k in range(nrw[ichk]): - if nmx[ichk, k] - 1 == j: continue - Tmn *= tanhtoc[k, ichk] - y = numpy.arctanh(-Tmn) - #y = platanh(-Tmn) - tov[i, j] = 2*y + // Send messages from bits to check nodes + for (int i = 0; i < M; ++i) { + for (int j = 0; j < nrw[i]; ++j) { + int ibj = Nm[i][j] - 1; + toc[i][j] = zn[ibj]; + for (int kk = 0; kk < 3; ++kk) { + // subtract off what the bit had received from the check + if (Mn[ibj][kk] - 1 == i) { + toc[i][j] -= tov[ibj][kk]; + } + } + } + } + + // send messages from check nodes to variable nodes + for (int i = 0; i < M; ++i) { + for (int j = 0; j < nrw[i]; ++j) { + //toc[i][j] = pltanh(-toc[i][j] / 2); + toc[i][j] = fast_tanh(-toc[i][j] / 2); + //toc[i][j] = tanhf(-toc[i][j] / 2); + } + } - return numpy.array([]) -*/ \ No newline at end of file + for (int i = 0; i < N; ++i) { + for (int j = 0; j < 3; ++j) { + int ichk = Mn[i][j] - 1; // Mn(:,j) are the checks that include bit j + float Tmn = 1.0f; + for (int k = 0; k < nrw[ichk]; ++k) { + if (Nm[ichk][k] - 1 != i) { + Tmn *= toc[ichk][k]; + } + } + //tov[i][j] = 2 * platanh(-Tmn); + tov[i][j] = 2 * fast_atanh(-Tmn); + //tov[i][j] = 2 * atanhf(-Tmn); + } + } + } + + *ok = best_score; +} diff --git a/ft8/ldpc.h b/ft8/ldpc.h index ebb831f..5928b3e 100644 --- a/ft8/ldpc.h +++ b/ft8/ldpc.h @@ -4,4 +4,6 @@ // plain is a return value, 174 ints, to be 0 or 1. // iters is how hard to try. // ok == 87 means success. -void ldpc_decode(float codeword[], int iters, int plain[], int *ok); +void ldpc_decode(float codeword[], int max_iters, int plain[], int *ok); + +void bp_decode(float codeword[], int max_iters, int plain[], int *ok);