Added belief propagation decoder as per WSJT-X

This commit is contained in:
Karlis Goba 2018-11-14 17:06:32 +02:00
parent db8f28dd8b
commit a8d963cdfd
4 changed files with 214 additions and 146 deletions

View file

@ -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);
}
}

View file

@ -1,7 +1,9 @@
#include <stdint.h>
// 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
};

View file

@ -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([])
*/
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;
}

View file

@ -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);