kopia lustrzana https://github.com/pabr/leansdr
FEC for DVB-S2 support
pabr
rodzic
3ece9cf9aa
commit
c49abe2707
|
|
@ -0,0 +1,216 @@
|
|||
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
|
||||
// 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
|
||||
#ifndef LEANSDR_BCH_H
|
||||
#define LEANSDR_BCH_H
|
||||
|
||||
#include "leansdr/discrmath.h"
|
||||
|
||||
namespace leansdr {
|
||||
|
||||
// Interface to hide the template parameters
|
||||
|
||||
struct bch_interface {
|
||||
virtual void encode(const uint8_t *msg, size_t msgbytes, uint8_t *out) = 0;
|
||||
virtual int decode(uint8_t *cw, size_t cwbytes) = 0;
|
||||
}; // bch_interface
|
||||
|
||||
// BCH error correction.
|
||||
// T: Unsigned type for packing binary polynomials.
|
||||
// N: Number of parity bits.
|
||||
// NP: Width of the polynomials supplied.
|
||||
// DP: Actual degree of the minimum polynomials (all must be same).
|
||||
// TGF: Unsigned type for syndromes (must be wider than DP).
|
||||
// GFTRUNCGEN: Generator polynomial for GF(2^DP), with X^DP omitted.
|
||||
|
||||
template<typename T, int N, int NP, int DP, typename TGF, int GFTRUNCGEN>
|
||||
struct bch_engine : bch_interface {
|
||||
bch_engine(const bitvect<T,NP> *polys, int _npolys)
|
||||
: npolys(_npolys)
|
||||
{
|
||||
// Build the generator polynomial (product of polys[]).
|
||||
g = 1;
|
||||
for ( int i=0; i<npolys; ++i )
|
||||
g = g * polys[i];
|
||||
// Convert the polynomials to truncated representation
|
||||
// (with X^DP omitted) for use with divmod().
|
||||
truncpolys = new bitvect<T,DP>[npolys];
|
||||
for ( int i=0; i<npolys; ++i )
|
||||
truncpolys[i].copy(polys[i]);
|
||||
// Check which polynomial contains each root.
|
||||
// Note: The DVB-S2 polynomials are numbered so that
|
||||
// syndpoly[2*i]==i, but we don't use that property.
|
||||
syndpolys = new int[2*npolys];
|
||||
for ( int i=0; i<2*npolys; ++i ) {
|
||||
int j;
|
||||
for ( j=0; j<npolys; ++j )
|
||||
if ( ! eval_poly(truncpolys[j], true, 1+i) ) break;
|
||||
if ( j == npolys ) fail("Bad polynomials/root");
|
||||
syndpolys[i] = j;
|
||||
}
|
||||
}
|
||||
|
||||
// Generate BCH parity bits.
|
||||
|
||||
void encode(const uint8_t *msg, size_t msgbytes, uint8_t *out) {
|
||||
bitvect<T,N> parity = shiftdivmod(msg, msgbytes, g);
|
||||
// Output as bytes, coefficient of highest degree first
|
||||
for ( int i=N/8; i--; ++out )
|
||||
*out = parity.v[i/sizeof(T)] >> ((i&(sizeof(T)-1))*8);
|
||||
}
|
||||
|
||||
// Decode BCH.
|
||||
// Return number of bits corrected, or -1 on failure.
|
||||
|
||||
int decode(uint8_t *cw, size_t cwbytes) {
|
||||
//again:
|
||||
bool corrupted = false;
|
||||
// Divide by individual polynomials.
|
||||
// TBD Maybe do in parallel, scanning cw only once.
|
||||
bitvect<T,DP> rem[npolys];
|
||||
for ( int j=0; j<npolys; ++j ) {
|
||||
rem[j] = divmod(cw, cwbytes, truncpolys[j]);
|
||||
}
|
||||
// Compute syndromes.
|
||||
TGF S[2*npolys];
|
||||
for ( int i=0; i<2*npolys; ++i ) {
|
||||
// Compute R(alpha^(1+i)), exploiting the fact that
|
||||
// R(x)=Q(x)g_j(X)+rem_j(X) and g_j(alpha^(1+i))=0
|
||||
// for some j that we already determined.
|
||||
// TBD Compute even exponents using conjugates.
|
||||
S[i] = eval_poly(rem[syndpolys[i]], false, 1+i);
|
||||
if ( S[i] ) corrupted = true;
|
||||
}
|
||||
if ( ! corrupted ) return 0;
|
||||
#if 0
|
||||
fprintf(stderr, "synd:");
|
||||
for ( int i=0; i<2*npolys; ++i ) fprintf(stderr, " %04x", S[i]);
|
||||
fprintf(stderr, "\n");
|
||||
#endif
|
||||
|
||||
// S_j = R(alpha_j) = 0+E(alpha_j) = sum(l=1..L)((alpha^j)^i_l)
|
||||
// where i_1 .. i_L are the degrees of the non-zero coefficients of E.
|
||||
// S_j = sum(l=1..L)((alpha^i_l)^j) = sum(l=1..L)(X_l^j)
|
||||
|
||||
// Berlekamp - Massey
|
||||
// http://en.wikipedia.org/wiki/Berlekamp%E2%80%93Massey_algorithm
|
||||
// TBD More efficient to work with logs of syndromes ?
|
||||
|
||||
int NN = 2*npolys;
|
||||
TGF C[NN]={1,}, B[NN]={1,};
|
||||
int L=0, m=1;
|
||||
TGF b = 1;
|
||||
for ( int n=0; n<NN; ++n ) {
|
||||
TGF d = S[n];
|
||||
for ( int i=1; i<=L; ++i ) d = GF.add(d, GF.mul(C[i],S[n-i]));
|
||||
if ( d == 0 )
|
||||
++m;
|
||||
else {
|
||||
TGF d_div_b = GF.mul(d, GF.inv(b));
|
||||
if ( 2*L <= n ) {
|
||||
TGF tmp[NN];
|
||||
memcpy(tmp, C, sizeof(tmp));
|
||||
for ( int i=0; i<NN-m; ++i )
|
||||
C[m+i] = GF.sub(C[m+i], GF.mul(d_div_b, B[i]));
|
||||
L = n + 1 - L;
|
||||
memcpy(B, tmp, sizeof(B));
|
||||
b = d;
|
||||
m = 1;
|
||||
} else {
|
||||
for ( int i=0; i<NN-m; ++i )
|
||||
C[m+i] = GF.sub(C[m+i], GF.mul(d_div_b, B[i]));
|
||||
++m;
|
||||
}
|
||||
}
|
||||
}
|
||||
// L is the number of errors.
|
||||
// C of degree L is the error locator polynomial (Lambda).
|
||||
// C(X) = sum(l=1..L)(1-X_l*X).
|
||||
#if 0
|
||||
fprintf(stderr, "C[%d]=", L);
|
||||
for ( int i=0; i<NN; ++i ) fprintf(stderr, " %04x", C[i]);
|
||||
fprintf(stderr, "\n");
|
||||
#endif
|
||||
|
||||
// Forney
|
||||
// http://en.wikipedia.org/wiki/Forney_algorithm
|
||||
// Simplified because coefficients are in GF(2).
|
||||
|
||||
// Find zeroes of C by exhaustive search.
|
||||
// TODO Chien method
|
||||
int roots_found = 0;
|
||||
for ( int i=0; i<(1<<DP)-1; ++i ) {
|
||||
// Candidate root ALPHA^i
|
||||
TGF v = eval_poly(C, L, i);
|
||||
if ( ! v ) {
|
||||
// ALPHA^i is a root of C, i.e. the inverse of an X_l.
|
||||
int loc = (i ? (1<<DP)-1-i : 0); // exponent of inverse
|
||||
// Reverse because cw[0..cwbytes-1] is stored MSB first
|
||||
int rloc = cwbytes*8-1 - loc;
|
||||
if ( rloc < 0 ) {
|
||||
// This may happen if the code is used truncated.
|
||||
return -1;
|
||||
}
|
||||
cw[rloc/8] ^= 128 >> (rloc&7);
|
||||
++roots_found;
|
||||
if ( roots_found == L ) break;
|
||||
}
|
||||
}
|
||||
if ( roots_found != L ) return -1;
|
||||
return L;
|
||||
}
|
||||
|
||||
private:
|
||||
|
||||
// Eval a GF(2)[X] polynomial at a power of ALPHA.
|
||||
|
||||
TGF eval_poly(const bitvect<T,DP> &poly, bool is_trunc, int rootexp) {
|
||||
TGF acc = 0;
|
||||
int re = 0;
|
||||
for ( int i=0; i<DP; ++i ) {
|
||||
if ( poly[i] ) acc = GF.add(acc, GF.exp(re));
|
||||
re += rootexp;
|
||||
if ( re >= (1<<DP)-1 ) re -= (1<<DP)-1; // mod 2^DP-1 incrementally
|
||||
}
|
||||
if ( is_trunc ) acc = GF.add(acc, GF.exp(re));
|
||||
return acc;
|
||||
}
|
||||
|
||||
// Eval a GF(2^16)[X] polynomial at a power of ALPHA.
|
||||
|
||||
TGF eval_poly(const TGF *poly, int deg, int rootexp) {
|
||||
TGF acc = 0;
|
||||
int re = 0;
|
||||
for ( int i=0; i<=deg; ++i ) {
|
||||
acc = GF.add(acc, GF.mul(poly[i],GF.exp(re)));
|
||||
re += rootexp;
|
||||
if ( re >= (1<<DP)-1 ) re -= (1<<DP)-1; // mod 2^DP-1 incrementally
|
||||
}
|
||||
return acc;
|
||||
}
|
||||
|
||||
bitvect<T,DP> *truncpolys;
|
||||
int npolys;
|
||||
int *syndpolys;
|
||||
bitvect<T,N> g;
|
||||
// Finite group for syndrome calculations
|
||||
gf2n<TGF, DP, 2, GFTRUNCGEN> GF;
|
||||
}; // bch_engine
|
||||
|
||||
} // namespace
|
||||
|
||||
#endif // LEANSDR_BCH_H
|
||||
|
|
@ -0,0 +1,432 @@
|
|||
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
|
||||
// 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
|
||||
#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<typename Taddr>
|
||||
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<typename SOFTBIT, typename SOFTWORD, int SWSIZE, typename Taddr>
|
||||
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<Taddr> *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<Taddr>::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<Taddr> *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<n_k/SWSIZE; ++i ) softword_zero(&parity[i]);
|
||||
|
||||
// Iterate over rows
|
||||
for ( const typename ldpc_table<Taddr>::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 ( ! msgbit ) continue; // TBD Generic soft version
|
||||
const Taddr *pa = prow->cols;
|
||||
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<Taddr> *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<Taddr>::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<SOFTBIT,SOFTWORD>(m,v) * PPCM / vnodes[v].nedges;
|
||||
//fprintf(stderr, "c[%d] bad => v[%d] += %d (%d*%d)\n",
|
||||
///c, v, s, err, softwords_weight<SOFTBIT,SOFTWORD>(m,*pe));
|
||||
score[v] += s;
|
||||
total += s;
|
||||
}
|
||||
}
|
||||
}
|
||||
return total;
|
||||
}
|
||||
|
||||
int decode_bitflip(const ldpc_table<Taddr> *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<SOFTBIT,SOFTWORD>(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<SOFTBIT,SOFTWORD>(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<SOFTBIT,SOFTWORD>(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<SOFTBIT,SOFTWORD>(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
|
||||
|
||||
#endif // LEANSDR_LDPC_H
|
||||
Ładowanie…
Reference in New Issue