gad.c

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

/* GAd.c, was bnd.c                                      (c) DJCM 95 08 31

   - belief net decoder  -

   Fri 16/2/01 - added gallager A decoding options to this and
   saved as GAd.c

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

   ------
   bnd.c:
   ------

   This is a set of routines to do Gallager / MS like  iterative 
 solution of 

 A x = z 

 defined by 

 A, represented in an alist_matrix. 
 a prior on x, given in a vector bias
 an initial condition

 No use is made of the true answer.
 
*/
/*
   How a bit in the A matrix is labelled.

   A one bit is labelled in two ways. 

   1) Its (m,n) pair, i.e. its row (m) and column (n).

   2) Its l and u values, which counts how many 1 bits it is along the row (l)
       and how many 1 bits it is down the column (u).

   A bit can be identified by its (m,l) value, its (n,u) value, or by its (m,n)
   value.

   e.g. the last bit in the bottom row of the (6*12) matrix below 
        has m=6 and n=10, and u=2, and l=4. 

   A = [ 0 1 0 0 0 1 0 1 0 0 0 1 ]
       [ 1 0 0 0 1 0 1 0 0 0 1 0 ]
       [ 0 0 1 1 0 0 0 1 0 0 0 1 ]
       [ 0 1 0 0 1 0 0 0 0 1 1 0 ]
       [ 1 0 0 1 0 1 0 0 1 0 0 0 ]
       [ 0 0 1 0 0 0 1 0 1 1 0 0 ]

   Notice that all rows in this matrix have equal weight lmax=4 and all 
   columns have equal weight umax=2 .

   A = [   1       1   1       1 ]
       [ 1       1   1       1   ]
       [     1 1       1       1 ]
       [   1     1         1 1   ]
       [ 1     1   1     1       ]
       [     1       1   1 1     ]

*/
/* CONTENTS

void GAd_allocate ( GAd_param *p , GAd_control *c ) {
void GAd_defaults ( GAd_param *p , GAd_control *GAdc ) {
void GAd_free ( GAd_param *p , GAd_control *c ) {
void GAd_load_dqc ( GAd_param *p , GAd_control *c ) {
void GAd_fprint_state ( GAd_param *p , GAd_control *c ) {
int GAdecode ( GAd_param *p , GAd_control *c ) {
void GAd_score_state ( GAd_param *p ) {
void GAhorizontal_pass ( GAd_param *p ) 
void GAvertical_pass ( GAd_param *p , GAd_control *c  ) 
void GAd_tweak_prior ( GAd_param *p , GAd_control *c  ) 

*/


#include "./ansi/r.h"
#include "./ansi/mynr.h"
#include "./ansi/cmatrix.h"
#include "./GAd.h"

/* malloc memory space */
void GAd_allocate ( GAd_param *p , GAd_control *c ) {
  int M = p->M ; 
  int N = p->N ;
  int m ; 

/* 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->bias = dvector ( 1 , N ) ; */ /* cut out in mnc7 onwards; this is
                                         allocated outside -- NB this makes
					 this code incompatible with mnc6.c */
  /*  p->q1 = cvector ( 1 , N ) ;  */ /* done by hook in gallagerA */
  p->pc1 = imatrix ( 1 , N , 1 , p->a->biggest_num_n ) ;
  p->dqc = cmatrix ( 1 , N , 0 , p->a->biggest_num_n + 1 ) ;

  p->qt1 = ivector ( 0 , p->a->biggest_num_n + 1 ) ;
  p->qb1 = ivector ( 0 , p->a->biggest_num_n + 1 ) ;

  p->dpf = cmatrix ( 1 , M , 0 , p->a->biggest_num_m + 1 ) ;
  p->dpr = cmatrix ( 1 , M , 0 , p->a->biggest_num_m + 1 ) ;

  for ( m = 1 ; m <= M ; m ++ ) {
    p->dpf[m][0] = 1 ; 
    p->dpr[m][p->a->num_mlist[m]+1 ] = 1 ; 
  }
  /* vertical passes must be initialized too  - done in the vertical pass */
}

void GAd_defaults ( GAd_param *p , GAd_control *GAdc ) {
#include "GAd_var_def.c"
}

void GAd_free ( GAd_param *p , GAd_control *c ) {
  int M = p->M ; 
  int N = p->N ;

  free_cvector ( p->t , 1, M ) ; 

  /*  free_cvector (   p->q1 , 1 , N ) ;  */

  free_cmatrix ( p->dpf , 1 , M , 0 , p->a->biggest_num_m + 1 ) ;
  free_cmatrix ( p->dpr , 1 , M , 0 , p->a->biggest_num_m + 1 ) ;

  free_imatrix ( p->pc1 , 1 , N , 1 , p->a->biggest_num_n ) ;

  free_cmatrix ( p->dqc , 1 , N , 0 , p->a->biggest_num_n + 1 ) ;

  free_ivector ( p->qt1 , 0 , p->a->biggest_num_n + 1 ) ;
  free_ivector ( p->qb1 , 0 , p->a->biggest_num_n + 1 ) ;

}

void   GAd_load_dqc ( GAd_param *p , GAd_control *c ) {
  int u , n , N=p->N ; 

  for ( n = 1 ; n <= N ; n++ )  {
    for ( u = 1 ; u <= p->a->biggest_num_n ; u++ )  {
      p->dqc[n][u] =  (p->bias[n]>0) ? 1 : 0 ; /* assume bias is +/- 1 */
    }
  }

}

int GAdecode ( GAd_param *p , GAd_control *c ) {

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

  c->include_bias = 2 ; /* standard value (assuming j=4) */

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

    GAhorizontal_pass ( p ) ;  /* turns binary dq vector
				  into dpr and dpf 
				  then reads out dpc -> pc1 */
    GAvertical_pass ( p , c ) ;    /* turns  pc1 into qt1, qb1, etc, 
				      and works out qc0, qc1 -> dqc.
				 Also works out the latest marginal belief */
    GAd_score_state (  p ) ;

    if ( c->GallagerB && ( c->loop ==  c->GallagerB ) ) {
      /* the time has come to change the bias contribution */
      c->include_bias = 0 ;
      fprintf(stderr,".B");
    }

    if ( c->writelog ) 
      GAd_fprint_state ( p , c ) ; 
    /* added 1997 02 */
    if ( c->wis && !(c->loop%c->wis) ) 
      GAd_write_state ( p , c ) ;  /* write the whole intermediate state */
  }
  c->loop -- ; 
  printf( "\tGAd:\tviols %3d  recon %3d  its %2d\n", p->count_viol , p->count_high  , c->loop ) ; 

  return ( p->count_viol ) ; /* zero if success */
}

void GAd_write_state ( GAd_param *p , GAd_control *c ) {
  FILE *fp ;
  int n ;
  unsigned char *q1 = p->q1 ; 
  int l = c->loop ; 
  char junk[1000] ; 
 
  sprintf ( junk , "%s.%02d" , c->isfile , l ) ;
  fp = fopen ( junk , "w" ) ;
  fprintf ( stdout , "writing intermediate state %d to %s\n" , l , junk  ) ; 
  if ( fp ) {
    for ( n = 1 ; n <= p->N ; n++ ) {
      fprintf ( fp ,  "%-1d " , q1[n] ) ; 
      if (!(n%c->wisline)) fprintf ( fp ,  "\n" ) ; 
    }
    fprintf ( fp ,  "\n" ) ; 
    fclose ( fp ) ; 
  }
  else 
    fprintf ( stderr , "can't open %s \n" , junk ) ;
}

void GAd_fprint_state ( GAd_param *p , GAd_control *c ) {
  FILE *fp ;
  int n ;
  unsigned char *q1 = p->q1 ; 

  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 GAhorizontal_pass ( GAd_param *p ) {
  int l , u , umax ; 
  int m , n ;
  unsigned char *dpf , *dpr ;
  unsigned char *dqc ;
  unsigned char dpc ; 
  unsigned char **pdpf = p->dpf ;
  unsigned char **pdpr = p->dpr ;
  int *lupto  = p->a->l_up_to ;     /* keeps track of what l we are up to */
  int *nlist ; 
  alist_matrix *a = p->a ;
  unsigned char *z = p->z ; 
  int **pc1 = p->pc1 ;

                                  /* * * * * * FORWARD PASS  * * * * * */
  for ( m = 1 ; m <= p->M ; m++ )   
    lupto[m] = 0 ;                      /* check all columns are at start    */

  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     */
    nlist = a->nlist[n] ;               /* list of bits in this column       */
    dqc = p->dqc[n] ;                   /* contributions to be made          */

    for ( u = 1 ; u <= umax ; u ++ ) {  /* Run down column              */
      m = nlist[u] ;                    /* what row are we on ?         */
      dpf = pdpf[m] ;                   /* -- for efficiency --         */

      l = lupto[m] ;                    /* what value of l = left-right */
      lupto[m] ++ ;                     /*   are we at on this row ?    */

      dpf[l+1] = dqc[u] ^ dpf[l] ;      /* the forward step             */
    }
  }
                                        /* * * * * * BACKWARD PASS * * * * * */

  for ( n = p->N ; n >= 1 ; n-- ) {     /* Run back through bits n           */
    umax = a->num_nlist[n] ;            /* number of bits in this column     */
    nlist = a->nlist[n] ;               /* list of bits in this column       */
    dqc = p->dqc[n] ;                   /* contributions to be made          */
    

    for ( u = 1 ; u <= umax ; u ++ ) {  /* Run down column              */
      m = nlist[u] ;                    /* what row are we on ?         */
      dpr = pdpr[m] ;                   /* -- for efficiency --         */
      dpf = pdpf[m] ;                   

      l = lupto[m] ;                    /* what value of l = left-right */
      lupto[m] -- ;                     /*   are we at on this row ?    */

      dpr[l] = dqc[u] ^ dpr[l+1] ;      /* the backward step             */
      dpc = dpf[l-1] ^ dpr[l+1] ^ z[m] ; /* the desired difference        */
      pc1[n][u] = ( dpc ) ?	 1 : -1 ;
    }
  }
} 


   
void GAd_score_state ( GAd_param *p ) { 
  unsigned char *q1 = p->q1 ; 
  int N = p->N ;
  int M = p->M ;
  int m , n ;

  for ( p->count_high = 0 , n = 1 ; n <= N ; n++ )  {
    p->count_high += ( q1[n] > 0 ) ? 1 : 0 ; 
  } 
  
  alist_times_cvector_mod2 ( (p->a) , p->q1 , p->t ) ; 
  
  for ( p->count_viol = 0 , m = 1 ; m <= M ; m ++ ) {
    if ( p->t[m] != p->z[m] ) p->count_viol ++ ; 
  }
  if ( p->count_viol == 0 ) p->not_decoded = 0 ; 
  
}


void GAvertical_pass ( GAd_param *p , GAd_control *c  ) 
{
  int u , umax ; 
  int n ;
  unsigned char *dqc ;
  int  *qt1 = p->qt1 ;
  int         *qb1 = p->qb1 ; 
  int **pc1 = p->pc1 ;
  int sum  ; 
  alist_matrix *a = (p->a) ;

  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   */
    dqc = p->dqc[n] ;                     /* contributions to be computed    */

    if ( c->include_bias ) { /* note bias is +/-1 */
      qt1[0] = p->bias[n] * c->include_bias ;
    } else {
      qt1[0] = 0 ;
    }
    qb1[umax+1] = 0 ;
                                        /* * * * * * DOWNWARD PASS * * * * * */
    for ( u = 1 ; u <= umax ; u ++ ) {  /* Run down column                   */
      qt1[u] = qt1[u-1] + pc1[n][u] ;   /*     downward step                 */
    }
    sum =  qt1[umax]  ;
    p->q1[n] =  (sum  > 0 ) ? 1 : 0 ;

                                        /* * * * * * UPWARD PASS * * * * * * */
    for ( u = umax ; u >= 1 ; u -- ) {  /* Run up column                     */
      qb1[u] = qb1[u+1] + pc1[n][u] ;   /*    upward step                    */

      sum = qt1[u-1] + qb1[u+1] ;    /*    everyone else's messages       */
      /* this is the net vote */
      dqc[u] = (sum>0)? 1 : 0;

    }
  }
}