rmbnd2.c

来自「快速傅立叶变换程序代码,学信号的同学,可要注意了」· C语言 代码 · 共 373 行

/* RMbnd2.c                                      (c) DJCM 95 08 31

   - belief net decoder for RM-based codes  - 
   RMbnd2 does codes with various trellises

   This code is (c) David J.C. MacKay 1995, 1997. It is free software 
   as defined by the free software foundation. 

   ------
see also   bnd.c:
   ------
   */

#include "./ansi/r.h"
#include "./ansi/mynr.h"
#include "./ansi/cmatrix.h"
#include "./RM.h"
#include "./RMbnd.h"

void bnd_allocate ( bnd_param *p , bnd_control *c ) {
  int M = p->M ; 
  int N = p->N ;

/* get z, xo to point to a vector allocated elsewhere; x is not used except in the printstate diagnostic */
  p->t = cvector ( 1 , M ) ; 

  p->qp1 = dvector ( 1 , N ) ; 

  p->r0 = dmatrix ( 1 , N , 1 , p->a->biggest_num_n ) ;
  p->r1 = dmatrix ( 1 , N , 1 , p->a->biggest_num_n ) ;

  p->q0 = dmatrix ( 1 , N , 0 , p->a->biggest_num_n + 1 ) ;
  p->q1 = dmatrix ( 1 , N , 0 , p->a->biggest_num_n + 1 ) ;

  p->qc0 = dvector ( 0 , p->a->biggest_num_n + 1 ) ;
  p->qc1 = dvector ( 0 , p->a->biggest_num_n + 1 ) ;
  p->qt0 = dvector ( 0 , p->a->biggest_num_n + 1 ) ;
  p->qt1 = dvector ( 0 , p->a->biggest_num_n + 1 ) ;
  p->qb0 = dvector ( 0 , p->a->biggest_num_n + 1 ) ;
  p->qb1 = dvector ( 0 , p->a->biggest_num_n + 1 ) ;

  p->trellism = ivector ( 1 , M ) ; 
  p->lts = ( linear_trellis ** ) malloc( (unsigned)(c->Ntrellis) * (sizeof(linear_trellis *)) )   ;
  p->lts -- ; 
  p->totms = ivector ( 1 , c->Ntrellis ) ; 
  set_ivector_const ( p->totms ,  1 , c->Ntrellis , 0 ) ; 
}

/* The following rather naff routine allocates rows to trellises on the 
   basis of the number of 1s per row of the alist. So no option of 
   having two trellises with same N!
*/
int set_up_trellism (  bnd_param *p, bnd_control *c ) {
  int m,M=p->M ;
  int *trellism = p->trellism ;
  int *totms = p->totms ;
  alist_matrix *a = p->a ;
  int *num_mlist = a->num_mlist ;    
  int status = 0 ; 

  for ( m = 1 ; m <= M ; m ++ ) {
    if ( c->trellisN == num_mlist[m] ) {
      trellism[m] = 1 ; totms[1] ++ ;
    }
    else if ( c->trellis2N == num_mlist[m] ) {
      trellism[m] = 2 ; totms[2] ++ ;
    }
    else if ( c->trellis3N == num_mlist[m] ) {
      trellism[m] = 3 ; totms[3] ++ ;
    }
    else {
      fprintf ( stderr , "error matching trellis m %d %d - fatal\n" , m , num_mlist[m] ) ; status -- ; 
    }
  }
  return status ;
}

int deduce_realM ( bnd_param *p , bnd_control *c ) {
  int bM = 0 , i , M=0 ;
  int *totms = p->totms ;
  linear_trellis *h ; 
  for ( i = 1 ; i <= c->Ntrellis ; i ++ ) {
    h =  p->lts[i] ;   /*  load the right trellis */
    bM += totms[i] * h->M ; 
    M += totms[i] ;
  }
  printf ( "deduced that %d lattice constraints correspond to %d binary parity checks\n" , M , bM ) ;
  return bM ;
}


void bnd_defaults ( bnd_param *p , bnd_control *bndc ) {
#include "RMbnd_var_def.c"
}

void bnd_free ( bnd_param *p , bnd_control *c ) {
  int M = p->M ; 
  int N = p->N ;

  free_cvector ( p->t , 1, M ) ; 
  
  free_dvector (   p->qp1 , 1 , N ) ; 

  free_dmatrix ( p->r0 , 1 , N , 1 , p->a->biggest_num_n ) ;
  free_dmatrix ( p->r1 , 1 , N , 1 , p->a->biggest_num_n ) ;

  free_dmatrix ( p->q0 , 1 , N , 0 , p->a->biggest_num_n +1 ) ;
  free_dmatrix ( p->q1 , 1 , N , 0 , p->a->biggest_num_n +1 ) ;

  free_dvector ( p->qc0 , 0 , p->a->biggest_num_n + 1 ) ;
  free_dvector ( p->qc1 , 0 , p->a->biggest_num_n + 1 ) ;
  free_dvector ( p->qt0 , 0 , p->a->biggest_num_n + 1 ) ;
  free_dvector ( p->qt1 , 0 , p->a->biggest_num_n + 1 ) ;
  free_dvector ( p->qb0 , 0 , p->a->biggest_num_n + 1 ) ;
  free_dvector ( p->qb1 , 0 , p->a->biggest_num_n + 1 ) ;

}

void   bnd_load_dqc ( bnd_param *p , bnd_control *c ) {
  int u , n , N=p->N ; 
  int verbose = 0 ; 

  for ( n = 1 ; n <= N ; n++ )  {
    for ( u = 1 ; u <= p->a->biggest_num_n ; u++ )  {
      p->q0[n][u] = 1.0 - p->bias[n] ;
      p->q1[n][u] = p->bias[n] ; /* is this right ? */
    }
    if ( verbose >= 3 ) {
      for ( u = 1 ; u <= p->a->biggest_num_n ; u++ )  {
	printf ( "prior q0/1[%d][%d] = %9.5g,%9.5g\n" , n , u , p->q0[n][u] ,p->q1[n][u] ) ; 
      }
    }
  }

}

int bndecode ( bnd_param *p , bnd_control *c ) {

  bnd_load_dqc ( p , c ) ; /* sets all dqc entries to the prior to start off */

  for ( c->loop = 1 , p->not_decoded = 1 ; 
       ( c->loop <= c->loops ) && p->not_decoded ; 
       c->loop ++ ) {

    horizontal_pass ( p ) ;  /* turns dqc into dpr and dpf 
			     then reads out dpc -> pc0 and pc1 */
    vertical_pass ( p , c ) ;    /* turns pc0 and pc1 into qt0, qb0, etc, 
			     and works out qc0, qc1 -> dqc.
			     Also works out the latest marginal belief */
    bnd_score_state ( p ) ;  /* takes belief vector, turns it into an
			     x, and runs it through Ax =? z .
			     If equal then decoding is halted.
			     */
    if ( c->writelog ) 
      bnd_fprint_state ( p , c ) ; 

  }
  c->loop -- ; 
  printf( "\t RM: \t viols %3d \t recon %3d \t its %2d\n", 
	   p->count_viol ,  p->count_high  , c->loop ) ; fflush(stdout);
  return ( p->count_viol ) ; /* zero if success */
}

void bnd_fprint_state ( bnd_param *p , bnd_control *c ) {
  FILE *fp ;
  int n ;
  double *q1 = p->qp1 ;  /* posteriors */

  fp = fopen ( c->logfile , "a" ) ;
  if ( fp ) {
    for ( n = 1 ; n <= c->writelog ; n++ ) 
      fprintf ( fp ,  "%-2.0f " , 99.0 * q1[n] ) ; 
    fprintf ( fp ,  "\n" ) ; 
    for ( n = 1 ; n <= c->writelog ; n++ ) 
      fprintf ( fp ,  "%-2d " , p->x[n] ) ; 
    fprintf ( fp ,  "\n" ) ; 
    fclose ( fp ) ; 
  }
  else 
    fprintf ( stderr , "can't open %s \n" , c->logfile ) ;
}

void bnd_score_state ( bnd_param *p ) { /* no longer assumes zero codeword */
  double *q1 = p->qp1 ; /* but it does assume zero syndrome */
  int N = p->N ;
  int M = p->M ;
  int l , m , n ;
  unsigned char *txo ; 
  unsigned char *xo = p->xo ;
  alist_matrix *a = p->a ;
  int lmax , *mlist ; 
  linear_trellis *h ; 
  int *trellism = p->trellism ;
  
  /* first threshold the posterior to get a decoding */
  for ( p->count_high = 0 , n = 1 ; n <= N ; n++ )  {
    xo[n] = ( q1[n] >= 0.5 ) ? 1 : 0 ;
    p->count_high += xo[n] ;
  }

  /* then behave like a horizontal pass, putting the relevant bits into
     the trellis xo */

  p->count_viol = 0 ;
  
  for ( m = 1 ; m <= M ; m++ ) {
    h =  p->lts[trellism[m]] ;   /*  load the right trellis */
    txo = h->xo ;
    /* load xo */
    lmax = a->num_mlist[m] ;            /* number of bits in this row     */
    if ( lmax != h->N ) {
      fprintf ( stderr , "%d != %d ! \n" , lmax, h->N ) ; fflush ( stderr ) ; 
    }
    mlist = a->mlist[m] ;               /* list of bits */

    for ( l = 1 ; l <= lmax ; l ++ ) {  /* Run along row              */
      n = mlist[l] ;                    /* what col are we on ?         */
      txo[l] = xo[n] ;                 /* */
    }                                  /* It's all loaded up */
    p->count_viol += lt_syndrome( h ) ? 1 : 0 ; 
  }

  if ( p->count_viol == 0 ) p->not_decoded = 0 ; 
}

/* in the future, give every row a pointer 
   to a trellis of its own - so as to have a variety of constraints? 
the future is now */
void horizontal_pass ( bnd_param *p ) 
{
  int l , u ; 
  int m , n ;
  double **tq ;
  double **tr ;
  alist_matrix *a = p->a ;
  int *u_up_to = a->u_up_to ; 
  int lmax , *mlist ; 
  double **q0 = p->q0 ; 
  double **q1 = p->q1 ; 
  double **r0 = p->r0 ; 
  double **r1 = p->r1 ; 
  linear_trellis *h ;
  int verbose ;
  int *trellism = p->trellism ;

  for ( n = 1 ; n <= p->N ; n++ )   
    u_up_to[n] = 1 ;                      /* check all columns are at start    */

  for ( m = 1 ; m <= p->M ; m++ ) {
    h =  p->lts[trellism[m]] ;   /*  load the right trellis */
    tq = h->q ;
    tr = h->r ;
    verbose = h->verbose ;

    /* load q */
    lmax = a->num_mlist[m] ;            /* number of bits in this row     */
    /* check that this matches */
    if ( lmax != h->N ) {
      fprintf ( stderr , "%d != %d ! \n" , lmax, h->N ) ; 
      fflush ( stderr ) ; 
    }
    mlist = a->mlist[m] ;               /* list of bits */

    for ( l = 1 ; l <= lmax ; l ++ ) {  /* Run along row              */
      n = mlist[l] ;                    /* what col are we on ?         */
      /* Find what the value of u is in this column */
      u = u_up_to[n] ;
      tq[l][0] = q0[n][u] ;  /* */
      tq[l][1] = q1[n][u] ;  /* */
    }                                  /* It's all loaded up */

    if(verbose >= 2) {
      printf ( "sending to lt_fb %d:\n" , m ) ;
      for ( l = 1 ; l <= lmax ; l ++ ) 
	printf ( "q[l] = %9.5g,%9.5g\n" , tq[l][0] , tq[l][1] ) ; 
    }

    lt_fb( h ) ;                         /* DO FORWARD BACKWARD */

    /* read out the r's */
    for ( l = 1 ; l <= lmax ; l ++ ) {  /* Run along row              */
      n = mlist[l] ;                    /* what col are we on ?         */
      u = u_up_to[n] ; 
      if ( u > 3 ) { fprintf ( stderr , "? u > 3 %d ?\n", u ) ; }
      u_up_to[n] ++ ; 
      r0[n][u] = tr[l][0] ; 
      r1[n][u] = tr[l][1] ; 
    }
  }
}

void vertical_pass ( bnd_param *p , bnd_control *c  ) 
{
  int u , umax ; 
  int n ;
  double  *q0 , *q1 , *qt0 = p->qt0 , *qt1 = p->qt1 ;
  double              *qb0 = p->qb0 , *qb1 = p->qb1 ; 
  double *qc0 = p->qc0 ;   double *qc1 = p->qc1 ;
  double **r0 = p->r0 ;  double **r1 = p->r1 ;
  double sum ; 
  alist_matrix *a = (p->a) ;
  int clipwarn=0 ; /* whether we have spat a clip warning this time */
  int fwarn=0 ; /* whether we have spat a division warning this time */

  for ( n = 1 ; n <= p->N ; n++ ) {     /* Run through bits x[n] left->right */
    umax = a->num_nlist[n] ;              /* number of bits in this column   */
    q0 = p->q0[n] ;                   /* new contributions to be computed    */
    q1 = p->q1[n] ;                     

    qt0[0] = 1.0 - p->bias[n] ;
    qt1[0] = p->bias[n] ;
    qb0[umax+1] = 1.0 ;
    qb1[umax+1] = 1.0 ;
                                        /* * * * * * DOWNWARD PASS * * * * * */
    for ( u = 1 ; u <= umax ; u ++ ) {  /* Run down column                   */
      qt0[u] = qt0[u-1] * r0[n][u] ;   /*     downward step                 */
      qt1[u] = qt1[u-1] * r1[n][u] ;   /*     downward step                 */
    }
    sum = ( qt0[umax] + qt1[umax] ) ;
    if ( sum > c->tinydiv ) {
      p->qp1[n] = qt1[umax] / sum ;      /* pseudoposterior prob, bit n */
    } else { 
      if ( !fwarn ) {
	fprintf ( stderr , "*" ) ; fflush ( stderr ) ; 
	fwarn = 1 ;
      }
    }
                                        /* * * * * * UPWARD PASS * * * * * * */
    for ( u = umax ; u >= 1 ; u -- ) {  /* Run up column                     */
      qb0[u] = qb0[u+1] * r0[n][u] ;   /*    upward step                    */
      qb1[u] = qb1[u+1] * r1[n][u] ;   /*    upward step                    */

      qc0[u] = qt0[u-1] * qb0[u+1] ;    /*    everyone else's messages       */
      qc1[u] = qt1[u-1] * qb1[u+1] ;    /*     -- not normalized         */

      sum = qc0[u] + qc1[u] ;
      if ( sum > c->tinydiv ) {
	q0[u] = qc0[u] / sum ;      
	q1[u] = qc1[u] / sum ;      
	if ( c->doclip ) { /* pragmatic clipping procedure to stop overflow */
	  if ( q0[u] > c->clip ) { 
	    q0[u] = c->clip ;
	    q1[u] = 1.0 - c->clip ;
	    if ( !clipwarn ) { fprintf ( stderr , "c" ) ; fflush ( stderr ) ; clipwarn=1; }
	  }
	  else if ( q1[u] > c->clip ) { 
	    q1[u] = c->clip ;
	    q0[u] = 1.0 - c->clip ;
	    if ( !clipwarn ) { 
	      fprintf ( stderr , "c" ) ; fflush ( stderr ) ; clipwarn=1; 
	    }
	  }
	}
	
      } else { 
	q0[u] = 0.5 ;
	q1[u] = 0.5 ;
	if ( !fwarn ) {
	  fprintf ( stderr , "*" ) ; fflush ( stderr ) ; 
	  fwarn = 1 ;
	}
      }
    }
  }
}


/*
<!-- hhmts start -->
Last modified: Tue Nov  7 15:55:46 1995
<!-- hhmts end -->
*/