dec.c

用C语言编写的LDPC译码程序

/* PRPRP.C - Procedures for performing probability propagation. */

/* Copyright (c) 2000 by Radford M. Neal 
 *
 * Permission is granted for anyone to copy, use, or modify this program 
 * for purposes of research or education, provided this copyright notice 
 * is retained, and note is made of any changes that have been made. 
 *
 * This program is distributed without any warranty, express or implied.
 * As this program was written for research purposes only, it has not been
 * tested to the degree that would be advisable in any important application.
 * All use of this program is entirely at the user's own risk.
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "alloc.h"
#include "mod2sparse.h"
#include "mod2dense.h"
#include "mod2convert.h"
#include "rcode.h"
#include "dec.h"
#include "enc.h"


/* DECODE BY EXHAUSTIVE ENUMERATION.  Decodes by trying all possible source
   messages (and hence all possible codewords, unless the parity check matrix
   was redundant).  If the last argument is 1, it sets dblk to the most likely
   entire block; if this argument is 0, each bit of dblk is set to the most
   likely value for that bit.  The marginal probabilities of each bit being 1
   are returned in bitpr, if it is not null.  

   The function value returned is the total number of codewords tried (which 
   will be the same for all blocks).  The return valued is "unsigned" because
   it might conceivably be as big as 2^31.

   The parity check matrix and other data are taken from the global variables
   declared in rcode.h.

   The number of message bits should not be greater than 31 for this procedure.
 */

unsigned enum_decode
( double *lratio,	/* Likelihood ratios for bits */
  char *dblk, 		/* Place to stored decoded message */
  double *bitpr,	/* Place to store marginal bit probabilities, or null */
  int max_block		/* Maximize probability of whole block being correct? */
)
{
  mod2dense *u, *v;
  double lk, maxlk, tpr;
  double *bpr, *lk0, *lk1;
  char sblk[31];
  char *cblk;
  unsigned d;
  int i, j;

  if (N-M>31) abort();

  /* Allocate needed space. */

  bpr = bitpr;
  if (bpr==0 && max_block==0)
  { bpr = chk_alloc (N, sizeof *bpr);
  }

  cblk = chk_alloc (N, sizeof *cblk);

  if (type=='d')
  { u = mod2dense_allocate(N-M,1);
    v = mod2dense_allocate(M,1);
  }

  if (type=='m')
  { u = mod2dense_allocate(M,1);
    v = mod2dense_allocate(M,1);
  }

  lk0 = chk_alloc (N, sizeof *lk0);
  lk1 = chk_alloc (N, sizeof *lk1);

  /* Pre-compute likelihoods for bits. */

  for (j = 0; j<N; j++)
  { lk0[j] = 1/(1+lratio[j]);
    lk1[j] = 1 - lk0[j];
  }

  /* Initialize marginal bit probabilities. */

  if (bpr)
  { for (j = 0; j<N; j++) bpr[j] = 0.0;
  }

  /* Exhaustively try all possible decoded messages. */

  tpr = 0.0;

  for (d = 0; d<=(1<<(N-M))-1; d++)
  {
    /* Unpack message into source block. */

    for (i = N-M-1; i>=0; i--)
    { sblk[i] = (d>>i)&1;
    }

    /* Find full codeword for this message. */

    switch (type)
    { case 's':
      { sparse_encode (sblk, cblk);
        break;
      }
      case 'd':
      { dense_encode (sblk, cblk, u, v);
        break;
      }
      case 'm':
      { mixed_encode (sblk, cblk, u, v);
        break;
      }
    }

    /* Compute likelihood for this decoding. */

    lk = 1;
    for (j = 0; j<N; j++)
    { lk *= cblk[j]==0 ? lk0[j] : lk1[j];
    }

    /* Update maximum likelihood decoding. */

    if (max_block)
    { if (d==0 || lk>maxlk)
      { for (j = 0; j<N; j++)
        { dblk[j] = cblk[j];
        }
        maxlk = lk;
      }
    }

    /* Update bit probabilities. */
 
    if (bpr)
    { for (j = 0; j<N; j++) 
      { if (cblk[j]==1) 
        { bpr[j] += lk;
        }
      }
      tpr += lk;
    }
  }

  /* Normalize bit probabilities. */

  if (bpr)
  { for (j = 0; j<N; j++) bpr[j] /= tpr;
  }

  /* Decoding to maximize bit-by-bit success, if that's what's wanted. */

  if (!max_block)
  { for (j = 0; j<N; j++) 
    { dblk[j] = bpr[j]>0.5;
    }
  }

  /* Free space. */

  if (bpr!=0 && bpr!=bitpr) free(bpr);
  free(cblk);
  free(lk0);
  free(lk1);

  return 1<<(N-M);
}


/* DECODE BY USING PROBABILITY PROPAGATION.  Tries to find the most probable 
   values for the bits of the codeword, given a parity check matrix (H), and
   likelihood ratios (lratio) for each bit.  If max_iter is positive, up to 
   that many iterations of probability propagation are done, stopping before 
   then if the tentative decoding is a valid codeword.  If max_iter is 
   negative, abs(max_iter) iterations are done, regardless of whether a 
   codeword was found earlier.

   Returns the number of iterations done (as an "unsigned" for consistency
   with enum_decode).  Regardless of whether or not a valid codeword was 
   reached, the bit vector from thresholding the bit-by-bit probabilities is 
   stored in dblk, and the resulting parity checks are stored in pchk (all 
   will be zero if the codeword is valid).  The final probability for each 
   bit being a 1 are also stored in bprb, unless this is null.
*/

unsigned prprp_decode
( mod2sparse *H,	/* Parity check matrix */
  double *lratio,	/* Likelihood ratios for bits */
  char *dblk,		/* Place to store decoding */
  char *pchk,		/* Place to store parity checks */
  double *bprb,		/* Place to store bit probabilities, 0 if not wanted  */
  int max_iter		/* Maximum number of iterations to perform, negative  */
)			/*  if should always perform abs(max_iter) iterations */
{ 
  int n;

  /* Initialize probability and likelihood ratios, and find initial guess. */

  initprp(H,lratio,dblk,bprb);

  /* Do up to abs(max_iter) iterations of probability propagation, stopping
     early if a codeword is found, unless max_iter is negative. */

  for (n = 0; 
       n!=max_iter && n!=-max_iter && (max_iter<0 || !check(H,dblk,pchk)); 
       n++)
  { 
    iterprp(H,lratio,dblk,bprb,(char*)0,(char*)0);
  }

  return n;
}


/* INITIALIZE PROBABILITY PROPAGATION.  Stores initial ratios, probabilities,
   and guess at decoding. */

void initprp
( mod2sparse *H,	/* Parity check matrix */
  double *lratio,	/* Likelihood ratios for bits */
  char *dblk,		/* Place to store decoding */
  double *bprb		/* Place to store bit probabilities, 0 if not wanted */
)
{ 
  mod2entry *e;
  int N;
  int j;

  N = mod2sparse_cols(H);

  for (j = 0; j<N; j++)
  { for (e = mod2sparse_first_in_col(H,j);
         !mod2sparse_at_end(e);
         e = mod2sparse_next_in_col(e))
    { e->pr = lratio[j];
      e->lr = 1;
    }
    if (bprb) bprb[j] = 1 - 1/(1+lratio[j]);
    dblk[j] = lratio[j]>1;
  }
}


/* DO ONE ITERATION OF PROBABILITY PROPAGATION. */

void iterprp
( mod2sparse *H,	/* Parity check matrix */
  double *lratio,	/* Likelihood ratios for bits */
  char *dblk,		/* Place to store decoding */
  double *bprb,		/* Place to store bit probabilities, 0 if not wanted */
  char *look_row,	/* Which rows to look at, look at all if null */
  char *look_col	/* Which columns to look at, look at all if null */
)
{
  double pr, dl, t;
  mod2entry *e;
  int N, M;
  int i, j;

  M = mod2sparse_rows(H);
  N = mod2sparse_cols(H);

  /* Recompute likelihood ratios. */

  for (i = 0; i<M; i++)
  { if (look_row!=0 && !look_row[i]) continue;
    dl = 1;
    for (e = mod2sparse_first_in_row(H,i);
         !mod2sparse_at_end(e);
         e = mod2sparse_next_in_row(e))
    { e->lr = dl;
      dl *= 2/(1+e->pr) - 1;
    }
    dl = 1;
    for (e = mod2sparse_last_in_row(H,i);
         !mod2sparse_at_end(e);
         e = mod2sparse_prev_in_row(e))
    { t = e->lr * dl;
      e->lr = (1-t)/(1+t);
      dl *= 2/(1+e->pr) - 1;
    }
  }

  /* Recompute probability ratios.  Also find the next guess based on the
     individually most likely values. */

  for (j = 0; j<N; j++)
  { if (look_col!=0 && !look_col[j]) continue;
    pr = lratio[j];
    for (e = mod2sparse_first_in_col(H,j);
         !mod2sparse_at_end(e);
         e = mod2sparse_next_in_col(e))
    { e->pr = pr;
      pr *= e->lr;
    }
    if (bprb) bprb[j] = 1 - 1/(1+pr);
    dblk[j] = pr>1;
    pr = 1;
    for (e = mod2sparse_last_in_col(H,j);
         !mod2sparse_at_end(e);
         e = mod2sparse_prev_in_col(e))
    { e->pr *= pr;
      pr *= e->lr;
    }
  }
}


/* SEE IF A GUESS CHECKS OUT.  Returns 1 if dblk contains a valid codeword.
   Whether or not it does, the results of all the parity checks are stored
   in pchk. */

int check
( mod2sparse *H,	/* Parity check matrix */
  char *dblk,		/* Guess for codeword */
  char *pchk		/* Place to store parity checks */
)
{
  int M;
  int i;

  M = mod2sparse_rows(H);

  mod2sparse_mulvec (H, dblk, pchk);

  for (i = 0; i<M; i++) 
  { if (pchk[i]!=0) return 0;
  }

  return 1;
}