// This file is part of LeanSDR Copyright (C) 2016-2018 . // See the toplevel README for more information. // // This program is free software: you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // This program is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with this program. If not, see . #ifndef LEANSDR_LDPC_H #define LEANSDR_LDPC_H #define lfprintf(...) \ { \ } namespace leansdr { // LDPC sparse matrix specified like in the DVB-S2 standard // Taddr must be wide enough to index message bits and address bits. template struct ldpc_table { // TBD Save space static const int MAX_ROWS = 162; // 64800 * (9/10) / 360 static const int MAX_COLS = 13; int q; int nrows; struct row { int ncols; Taddr cols[MAX_COLS]; } rows[MAX_ROWS]; }; // LDPC ENGINE // SOFTBITs can be hard (e.g. bool) or soft (e.g. llr_t). // They are stored as SOFTWORDs containing SWSIZE SOFTBITs. // See interface in softword.h. template struct ldpc_engine { ldpc_engine() : vnodes(NULL), cnodes(NULL) { } // vnodes: Value/variable nodes (message bits) // cnodes: Check nodes (parity bits) int k; // Message size in bits int n; // Codeword size in bits struct node { Taddr *edges; int nedges; static const int CHUNK = 4; // Grow edges[] in steps of CHUNK. void append(Taddr a) { if (nedges % CHUNK == 0) { // Full ? edges = (Taddr *)realloc(edges, (nedges + CHUNK) * sizeof(Taddr)); if (!edges) fatal("realloc"); } edges[nedges++] = a; } }; node *vnodes; // [k] node *cnodes; // [n-k] // Initialize from a S2-style table. ldpc_engine(const ldpc_table *table, int _k, int _n) : k(_k), n(_n) { // Sanity checks if (360 % SWSIZE) fatal("Bad LDPC word size"); if (k % SWSIZE) fatal("Bad LDPC k"); if (n % SWSIZE) fatal("Bad LDPC n"); if (k != table->nrows * 360) fatal("Bad table"); int n_k = n - k; if (table->q * 360 != n_k) fatal("Bad q"); vnodes = new node[k]; memset(vnodes, 0, sizeof(node) * k); cnodes = new node[n_k]; memset(cnodes, 0, sizeof(node) * n_k); // Expand the graph. int m = 0; // Iterate over rows for (const typename ldpc_table::row *prow = table->rows; prow < table->rows + table->nrows; ++prow) { // Process 360 bits per row. int q = table->q; int qoffs = 0; for (int mw = 360; mw--; ++m, qoffs += q) { const Taddr *pa = prow->cols; for (int nc = prow->ncols; nc--; ++pa) { int a = (int)*pa + qoffs; if (a >= n_k) a -= n_k; // Modulo n-k. Note qoffs<360*q. if (a >= n_k) fail("Invalid LDPC table"); vnodes[m].append(a); cnodes[a].append(m); } } } } void print_node_stats() { int nedges = count_edges(vnodes, k); fprintf(stderr, "LDPC(%5d,%5d)(%.2f)" " %5.2f edges/vnode, %5.2f edges/cnode\n", k, n - k, (float)k / n, (float)nedges / k, (float)nedges / (n - k)); } int count_edges(node *nodes, int nnodes) { int c = 0; for (int i = 0; i < nnodes; ++i) c += nodes[i].nedges; return c; } // k: Message size in bits // n: Codeword size in bits // integrate: Optional S2-style post-processing #if 0 void encode_hard(const ldpc_table *table, const uint8_t *msg, int k, int n, uint8_t *parity, bool integrate=true) { // Sanity checks if ( 360 % SWSIZE ) fatal("Bad LDPC word size"); if ( k % SWSIZE ) fatal("Bad LDPC k"); if ( n % SWSIZE ) fatal("Bad LDPC n"); if ( k != table->nrows*360 ) fatal("Bad table"); int n_k = n - k; if ( table->q*360 != n_k ) fatal("Bad q"); for ( int i=0; i::row *prow = table->rows; // quirk prow < table->rows+table->nrows; ++prow ) { // Process 360 bits per row, in words of SWSIZE bits int q = table->q; int qoffs = 0; for ( int mw=360/SWSIZE; mw--; ++msg ) { SOFTWORD msgword = *msg; for ( int wbit=0; wbitcols; for ( int nc=prow->ncols; nc--; ++pa ) { // Don't wrap modulo range of Taddr int a = (int)*pa + qoffs; // Note: qoffs < 360*q=n-k if ( a >= n_k ) a -= n_k; // TBD not predictable softwords_flip(parity, a); } } } } if ( integrate ) integrate_bits(parity, parity, n_k/SWSIZE); } #endif void encode(const ldpc_table *table, const SOFTWORD *msg, int k, int n, SOFTWORD *parity, int integrate = true) { // Sanity checks if (360 % SWSIZE) fatal("Bad LDPC word size"); if (k % SWSIZE) fatal("Bad LDPC k"); if (n % SWSIZE) fatal("Bad LDPC n"); if (k != table->nrows * 360) fatal("Bad table"); int n_k = n - k; if (table->q * 360 != n_k) fatal("Bad q"); for (int i = 0; i < n_k / SWSIZE; ++i) softword_zero(&parity[i]); // Iterate over rows for (const typename ldpc_table::row *prow = table->rows; // quirk prow < table->rows + table->nrows; ++prow) { // Process 360 bits per row, in words of SWSIZE bits int q = table->q; int qoffs = 0; for (int mw = 360 / SWSIZE; mw--; ++msg) { SOFTWORD msgword = *msg; for (int wbit = 0; wbit < SWSIZE; ++wbit, qoffs += q) { SOFTBIT msgbit = softword_get(msgword, wbit); if (!softbit_harden(msgbit)) continue; const Taddr *pa = prow->cols; for (int nc = prow->ncols; nc--; ++pa) { int a = (int)*pa + qoffs; // Note: qoffs < 360*q=n-k if (a >= n_k) a -= n_k; // TBD not predictable softwords_flip(parity, a); } } } } if (integrate) integrate_bits(parity, parity, n_k / SWSIZE); } // Flip bits connected to parity errors, one at a time, // as long as things improve and max_bitflips is not exceeded. // cw: codeword (k value bits followed by n-k check bits) static const int PPCM = 39; typedef int64_t score_t; score_t compute_scores(SOFTWORD *m, SOFTWORD *p, SOFTWORD *q, int nc, score_t *score, int k) { int total = 0; memset(score, 0, k * sizeof(*score)); for (int c = 0; c < nc; ++c) { SOFTBIT err = softwords_xor(p, q, c); if (softbit_harden(err)) { Taddr *pe = cnodes[c].edges; for (int e = cnodes[c].nedges; e--; ++pe) { int v = *pe; int s = err * softwords_weight(m, v) * PPCM / vnodes[v].nedges; //fprintf(stderr, "c[%d] bad => v[%d] += %d (%d*%d)\n", ///c, v, s, err, softwords_weight(m,*pe)); score[v] += s; total += s; } } } return total; } int decode_bitflip(const ldpc_table *table, SOFTWORD *cw, int k, int n, int max_bitflips) { if (!vnodes) fail("LDPC graph not initialized"); int n_k = n - k; // Compute the expected check bits (without the final mixing) SOFTWORD expected[n_k / SWSIZE]; encode(table, cw, k, n, expected, false); // Reverse the integrator mixing from the received check bits SOFTWORD received[n_k / SWSIZE]; diff_bits(cw + k / SWSIZE, received, n_k / SWSIZE); // Compute initial scores score_t score[k]; score_t tots = compute_scores(cw, expected, received, n_k, score, k); lfprintf(stderr, "Initial score %d\n", (int)tots); int nflipped = 0; score_t score_threshold; { SOFTBIT one; softbit_set(&one, true); score_threshold = (int)one * 2; } bool progress = true; while (progress && nflipped < max_bitflips) { progress = false; // Try to flip parity bits. // Due to differential decoding, they appear as consecutive errors. SOFTBIT prev_err = softwords_xor(expected, received, 0); for (int b = 0; b < n - k - 1; ++b) { prev_err = softwords_xor(expected, received, b); //TBD SOFTBIT err = softwords_xor(expected, received, b + 1); if (softbit_harden(prev_err) && softbit_harden(err)) { lfprintf(stderr, "flip parity %d\n", b); softwords_flip(received, b); softwords_flip(received, b + 1); ++nflipped; // Counts as one flip before differential decoding. progress = true; int dtot = 0; // Depenalize adjacent message bits. { Taddr *pe = cnodes[b].edges; for (int e = cnodes[b].nedges; e--; ++pe) { int d = prev_err * softwords_weight(cw, *pe) * PPCM / vnodes[*pe].nedges; score[*pe] -= d; dtot -= d; } } { Taddr *pe = cnodes[b + 1].edges; for (int e = cnodes[b + 1].nedges; e--; ++pe) { int d = err * softwords_weight(cw, *pe) * PPCM / vnodes[*pe].nedges; score[*pe] -= d; dtot -= d; } } tots += dtot; #if 1 // Also update the codeword in-place. // TBD Useful for debugging only. softwords_flip(cw, k + b); #endif // TBD version soft. err = ! err; } prev_err = err; } // c nodes score_t maxs = -(1 << 30); for (int v = 0; v < k; ++v) if (score[v] > maxs) maxs = score[v]; if (!maxs) break; lfprintf(stderr, "maxs %d\n", (int)maxs); // Try to flip each message bits with maximal score for (int v = 0; v < k; ++v) { if (score[v] < score_threshold) continue; // if ( score[v] < maxs*9/10 ) continue; if (score[v] < maxs - 4) continue; lfprintf(stderr, " flip %d score=%d\n", (int)v, (int)score[v]); // Update expected parities and scores that depend on them. score_t dtot = 0; for (int commit = 0; commit <= 1; ++commit) { Taddr *pe = vnodes[v].edges; for (int e = vnodes[v].nedges; e--; ++pe) { Taddr c = *pe; SOFTBIT was_bad = softwords_xor(expected, received, c); if (softbit_harden(was_bad)) { Taddr *pe = cnodes[c].edges; for (int e = cnodes[c].nedges; e--; ++pe) { int d = was_bad * softwords_weight(cw, *pe) * PPCM / vnodes[*pe].nedges; if (commit) score[*pe] -= d; else dtot -= d; } } softwords_flip(expected, c); SOFTBIT is_bad = softwords_xor(expected, received, c); if (softbit_harden(is_bad)) { Taddr *pe = cnodes[c].edges; for (int e = cnodes[c].nedges; e--; ++pe) { int d = is_bad * softwords_weight(cw, *pe) * PPCM / vnodes[*pe].nedges; if (commit) score[*pe] += d; else dtot += d; } } if (!commit) softwords_flip(expected, c); } if (!commit) { if (dtot >= 0) { lfprintf(stderr, " abort %d\n", v); break; // Next v } } else { softwords_flip(cw, v); ++nflipped; tots += dtot; progress = true; v = k - 1; // Force exit to update maxs ? } } // commit } // v lfprintf(stderr, "progress %d\n", progress); #if 0 fprintf(stderr, "CHECKING TOTS INCREMENT (slow) %d\n", tots); score_t tots2 = compute_scores(cw, expected, received, n_k, score, k); if ( tots2 != tots ) fail("bad tots update"); #endif } return nflipped; } // EN 302 307-1 5.3.2.1 post-processing of parity bits. // In-place allowed. #if 1 static void integrate_bits(const SOFTWORD *in, SOFTWORD *out, int nwords) { SOFTBIT sum; softbit_clear(&sum); for (int i = 0; i < nwords; ++i) { SOFTWORD w = in[i]; for (int b = 0; b < SWSIZE; ++b) { sum = softbit_xor(sum, softword_get(w, b)); softword_write(w, b, sum); } out[i] = w; } } #else // Optimized for hard_sb static void integrate_bits(const uint8_t *in, uint8_t *out, int nwords) { // TBD Optimize uint8_t prev = 0; for (int i = 0; i < nwords; ++i) { uint8_t c = in[i]; for (int j = SWSIZE; j--;) { c ^= prev << j; prev = (c >> j) & 1; } out[i] = c; } } #endif // Undo EN 302 307-1 5.3.2.1, post-processing of parity bits. // In-place allowed. #if 1 static void diff_bits(const SOFTWORD *in, SOFTWORD *out, int nwords) { SOFTBIT prev; softbit_clear(&prev); for (int i = 0; i < nwords; ++i, ++in, ++out) { SOFTWORD w = *in; for (int b = 0; b < SWSIZE; ++b) { SOFTBIT n = softword_get(w, b); softword_write(w, b, softbit_xor(prev, n)); prev = n; } *out = w; } } #else // Optimized for hard_sb static void diff_bits(const uint8_t *in, uint8_t *out, int nwords) { uint8_t prev = 0; for (int i = 0; i < nwords; ++i) { uint8_t c = in[i]; out[i] = c ^ (prev | (c >> 1)); prev = (c & 1) << (SWSIZE - 1); } } #endif }; // ldpc_engine } // namespace leansdr #endif // LEANSDR_LDPC_H