ldpc.h

Go to the documentation of this file.

// 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",
                    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 = new SOFTWORD[n_k / SWSIZE]; // Forbidden to statically allocate with non constant size
        encode(table, cw, k, n, expected, false);
        // Reverse the integrator mixing from the received check bits
        SOFTWORD *received = new SOFTWORD[n_k / SWSIZE]; // Forbidden to statically allocate with non constant size
        diff_bits(cw + k / SWSIZE, received, n_k / SWSIZE);

        // Compute initial scores
        score_t *score = new score_t[k]; // Forbidden to statically allocate with non constant size
        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
        }

        delete[] score;
        delete[] received;
        delete[] expected;

        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