rmbnd.c

快速傅立叶变换程序代码,学信号的同学,可要注意了

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

   - belief net decoder for RM-based codes  - 

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

   ( RMbnd2 does codes with multiple trellises )
   ------
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 ) ;

}

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 ) ; 
  return ( p->count_viol ) ; /* zero if believed a 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 = p->lt->xo ;
  unsigned char *xo = p->xo ;
  alist_matrix *a = p->a ;
  int lmax , *mlist ; 
  linear_trellis *h = p->lt ; /* see RMbnd2.c for multiple trellis version */
  int verbose = h->verbose ;
  
  /* 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++ ) {   /*  load the right trellis */
    /* 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? */
void horizontal_pass ( bnd_param *p ) 
{
  int l , u ; 
  int m , n ;
  double **tq = p->lt->q ;
  double **tr = p->lt->r ;
  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 = p->lt ; 
  int verbose = h->verbose ;

  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++ ) {
    /* 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: Sat Jun  7 17:44:09 1997
<!-- hhmts end -->
*/