ra.c

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

/*
   RA.c
              (c) DJCM 98 09 28

   Repeat-accumulate code simulator

 read in code definition
 loop {
   encode source string
   add noise
   decode
 }

 Code definition: (stored in "alist")

			 Use of alist allows arbitrary numbers of repetitions
			 of each bit.

  K                      source block length
  n_1 n_2 ... n_K        number of repetitions of each source bit
  N = sum n_k
  alist  defines         permutation of N encoded bits
                         note, an additional permutation of the N accumulated
			 bits may be a good idea. (for non-memoryless channels)

 transmitted bits are integral of encoded bits

 Future plans: 
  clump source bits into clumps. Have multiple parallel accumulated streams.
  Have little sub-matrices (like GF(q) ) defining response of accumulator to 
  clumps.
 
*/

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

#include "./RA.h" /* this defines  data_creation_param ; RA_control */

int RA_encode ( unsigned char * , RA_control * , unsigned char * ) ;
static int t_to_b (  unsigned char * ,	RA_control * ) ;
int RA_decode ( RA_control * ) ;
int RA_horizontal_pass ( RA_control * ) ;

static void dc_defaults ( RA_control * ) ; 
static int  process_command ( int , char **  , RA_control * ) ; 
static void print_usage ( char ** , FILE * ,  RA_control *  );
static int  make_sense ( RA_control * ) ; 
static int  make_space ( RA_control * ) ;
static int  score ( RA_control *  ) ; 

static void finalline (  FILE * , RA_control * , int ) ; /* int = 1 to get loads of info */

static double bern ( int  , int  , double * , double * ,  double * , double );
static void histo ( FILE * ,  RA_control * ) ;
static void snappyline ( RA_control * ) ;
static void RA_free ( RA_control * ) ; 
static int  check_alist_MN ( alist_matrix * , RA_control * ) ; 
static double h2 ( double ) ;

void   main ( int , char ** ) ;

/*
        MAIN
                     */
void main ( int argc, char *argv[] )
{
  FILE   *fp  ; 

  int k ; 
  RA_control         c ; 
  dc_defaults   ( &c ) ; 

  if ( process_command (argc, argv, &c ) < 0 ) exit (0) ;
  if ( read_allocate_alist ( &(c.a) , c.afile ) < 0 ) exit (0) ; 
  if ( check_alist_MN ( &(c.a) , &c ) < 0 ) exit (0) ; 
  if ( make_sense ( &c  ) < 0 ) exit (0) ;

  fprintf(stderr,"RA N=%d, K=%d, x=%6.3g xass=%6.3g fn=%6.3g fnass=%6.3g\n",
	  c.N , c.K , c.gcx , c.gcxass , c.fn , c.fnass ) ;
  fflush(stderr);

  if ( make_space ( &c ) < 0 ) exit (0) ;

  if ( c.writelog ) {
    fp = fopen ( c.logfile , "w" ) ;
    if ( !fp ) {
      fprintf ( stderr , " couldn't open logfile %s\n" , c.logfile ) ; 
      c.writelog = 0 ; 
    } else fclose (fp ) ; 
  }

  if ( c.writelog ) {
    fp = fopen ( c.logfile , "w" ) ;
    if ( !fp ) {
      fprintf ( stderr , " couldn't open logfile %s\n" , c.logfile ) ; 
      c.writelog = 0 ; 
    } else fclose (fp ) ; 
  }

  if ( c.error_log ) {
    fp = fopen ( c.error_logfile , "w" ) ;
    if ( !fp ) {
      fprintf ( stderr , " couldn't open logfile %s\n" , c.error_logfile ) ; 
      c.error_log = 0 ; 
    } else fclose (fp ) ; 
  }

  /*
      MAIN LOOP
                 */

  ran_seed ( c.vseed ) ; 
  c.message = 1 ; 
  for ( ;       ( c.message <= c.MESSAGE ) && 
	  ( ( c.failures==0 ) || ( c.failcount < c.failures ) ) ; 
	c.message ++ ) {

    snappyline( &c ) ;
    /* force parity bit at end */
    c.sourceweight = random_cvector ( c.s , c.fs , 1 , c.K ) ; 

    if ( c.verbose > 2) {
      printf ( "source vector:\n" ) ;
      for ( k = 1 ; k <= c.K ; k ++ ) {
	if ( c.s[k] ) printf ("1 ") ; else printf ( "0 " ); 
      }
      printf ( "\n" ) ; 
    }
    RA_encode ( c.s , &c , c.t ) ; 
    if ( c.verbose > 2) {
      printf ( "transmitted vector:\n" ) ;
      for ( k = 1 ; k <= c.N ; k ++ ) {
	if ( c.t[k] ) printf ("1 ") ; else printf ( "0 " ); 
      }
      printf ( "\n" ) ; 
    }
    c.flipped = t_to_b ( c.t , &c ) ; 
    if ( c.verbose > 2 ) {
      printf ( "received likelihoods:\n" ) ;
      for ( k = 1 ; k <= c.N ; k ++ ) {
	printf ("%1d " , (int) ( c.bias[k][1] * 10.0 ) ) ;
      }
      printf ( "\n" ) ; 
      for ( k = 1 ; k <= c.N ; k ++ ) {
	printf ("%1d " , (int) ( c.bias[k][1] * 2.0 ) ) ;
      }
      printf ( "\n" ) ; 
    }
    RA_decode ( &c ) ;

    if ( score ( &c ) < 0 ) exit ( 0 ) ; 
    if ( c.verbose > 0 ) finalline ( stdout , &c , 0 ) ;
    if ( c.printout ) { /* append */
      fp = fopen ( c.outfile , ((c.outappend)? "a":"w") ) ;
      if( !fp ) {
	fprintf( stderr, "No such file: %s\n", c.outfile ) ;
	finalline ( stderr , &c , 0 ) ;
      }      else 	{
	finalline ( fp ,  &c , 0 ) ;
	fclose ( fp ) ;
      }
    }
    if (  (!( ( c.message+1 <= c.MESSAGE ) && 
	  ( ( c.failures==0 )
	    || ( c.failcount < c.failures ) ) ) )
	  || (!(c.message % c.big_write_period ) ) ) {
      if ( c.printtot ) { /* write */
	fp = fopen ( c.totoutfile , "w" ) ;
	if( !fp ) {
	  fprintf( stderr, "No such file: %s\n", c.totoutfile ) ;
	  finalline ( stderr , &c , 1 ) ; /* totalline */
	}      else 	{
	  finalline ( fp , &c , 1 ) ;
	  fclose ( fp ) ;
	}
      }
      if ( c.printhisto && c.block_valid ) { /* update histogram file */
	fp = fopen ( c.histofile , "w" ) ;
	if( !fp ) {
	  fprintf( stderr, "No such file: %s\n", c.histofile ) ;
	}      else 	{
	  histo ( fp , &c ) ;
	  fclose ( fp ) ;
	}
      }
    }
  }
  snappyline( &c ) ; printf("\n") ; 
  RA_free (  &c ) ;
}

static void histo ( FILE *fp ,  RA_control *c ) {
  int l;
  double t , cum = 0.0 ;
  double  tot = (double) c->block_valid ; 

  fprintf ( fp , "# total valid blocks %d\n" , c->block_valid ) ; 
  for ( l = 1 ; l <= c->loops ; l ++ ) {
    t= (double) c->histo[l] ;
    cum += t ;
    fprintf ( fp , "%d\t%d\t%d\t%9.4g\t%9.4g\n" , l , (int)(t) , (int)(cum) , 
	      t/tot , cum/tot ) ; 
  }

}
static void snappyline ( RA_control *c ) {
  printf ( "%d:%du%dd%dl/%d\t" , c->block_errs , c->block_undet, c->block_det, c->block_detlw , c->message - 1 ) ; fflush ( stdout ) ; 

}

/*                                                  <<<<N>>>>
  Encoding method: 
  source bits d[1]..d[K] are mapped via an alist  ^ 1    1 1
  into a pre-transmission vector                  K  1 1  1 
  s[1]..s[N] .                                    V   1 1   1 

  t[n] = t[n-1] ^ s[n]

     s[n]      s[n+1]
      |         |
0 .. -+-> t[n] -+-> t[n+1] ..... t[N]
           |        |            |     
           V        V		 V     
          y[n]      y[n+1]	 y[N]

  */

int RA_encode ( unsigned char *d , RA_control *c , unsigned char *t ) {
  int n , k ; int status = 0 ; 
  alist_transpose_cvector_sparse_mod2 ( &c->a , d , t ) ; /* here 't' doubles as 's' */
  /* accumulate */

  if ( c->verbose > 2) {
    printf ( "extended source vector:\n" ) ;
    for ( k = 1 ; k <= c->N ; k ++ ) {
      if ( t[k] ) printf ("1 ") ; else printf ( "0 " ); 
    }
    printf ( "\n" ) ; 
  }

  for ( n = 2 ; n <= c->N ; n ++ ) {
    t[n] = t[n]^t[n-1] ; 
  }

  if ( c->verbose > 2) {
    printf ( "accumulated transmission:\n" ) ;
    for ( k = 1 ; k <= c->N ; k ++ ) {
      if ( t[k] ) printf ("1 ") ; else printf ( "0 " ); 
    }
    printf ( "\n" ) ; 
  }
  return status ; 
}

/* 
   The channel outputs a normalized likelihood vector 
   bias[n] = P( yn | tn = 1 )

   The state of the decoder is  q[1..N][0/1] and r[1..N][0/1]

   q[n][s] = pseudoprior( s[n] = s )     s=0/1      initially 0.5

   Use f/b algorithm to find:
                                       initial conditions: 
   f[n][t] = P( y1...yn , tn=t )       f[0][0] = 1 ; f[0][1] = 0 ; 
   b[n][t] = P( yn...yN | tn=t )       b[N+1][0] = 1 ; b[N+1][1] = 1 ;

   Using
   f[n][t] = bias[n][t] * sum_{t': t'+s=t} ( f[n-1][t'] pi[n][s] )
   b[n][t] = bias[n][t] * sum_{t': t'+s=t} ( b[n+1][t'] pi[n+1][s] )

   Find likelihood contribution at n:
      r[n][s] `=' P(y1..yN|s[n]=s) = sum_{s: t'+s=t} f[n-1][t] b[n][t'] 
      
   Then in the vertical step we visit each incarnation of the source bit
   for( r = 1 .. repetitions[k] ) (number on mlist)  n=nlist[r]
      d[r][s] = d[r-1][s] * r[n][s] ;
   reverse,
      u[r][s] = u[r+1][s] * r[n][s] ;
   ->
   pi[n][s] `=' P( s[n] = s | y1..yN ) = prior0(omit) * d[r-1][s] * u[r+1][s]

   If all l's have the same sign then the locally mp edges in the trellis
   form a valid codeword
   
   */

int RA_decode ( RA_control *c ) {
  int u , n , N=c->N ; 

  for ( n = 1 ; n <= N ; n++ )  {
    for ( u = 0 ; u <= 1 ; u++ )  {
      c->q[n][u] = 1.0 ; /* assume all data strings equiprobable */
    }
  }

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

    c->viols = RA_horizontal_pass ( c ) ;
    /*    if ( c->writelog )  RA_fprint_state ( c ) ;  */

  } 
  c->loop -- ; 
  printf( "\tviols %3d its %2d\n", 
	  c->viols ,  c->loop ) ; 

  return ( c->viols ) ; /* zero if success */
  
}

/* for efficiency, the backward pass variables are not stored; 
 they are immediately multiplied by the forward pass variables 
 to give likelihood signals.
 The likelihood signals are (at present) written into the same 
 array as the f's
 */

int RA_horizontal_pass ( RA_control *c ) {
  int  u , umax ; 
  int m , n , N=c->N ; /* ? */
  double **f = c->f ; 
  /* f has been initialized with f[0][0] = 1 ; f[0][1] = 0 ; */

  /*  double **b=c->b ;        b[N+1][0] = 1 ; b[N+1][1] = 1 ; */
  double **bias = c->bias ; /* likelihood */
  double **q = c->q ; 
  double **r = c->r ; 
  double p_t0 , p_t1 ; 
  int *mlist ; 
  int v0 , v1 , viols  ; /* votes for 0 and 1 ; viols sees if they disagree */
  /*  int TERMINATED = 0 ;  */
  double BOT = 0.99 ; 
  double TOP = 1.01 ; 
  alist_matrix *a = &(c->a) ;
  unsigned char *so = c->so ; 
  unsigned char *to = c->to ; 
  double f1 = 0.0 , f0 = 1.0 , b1 , b0 , p0 , p1 , q0 , q1 , s ; 
  int clipwarn=0 ; /* whether we have spat a clip warning this time */

  for ( n = 1 ; n <= N ; n++ )  { /* n=N is not really needed */
    q0 = q[n][0] ;    q1 = q[n][1] ; /* latest prior probabilities */
    p0 = q0 * f0 + q1 * f1 ; 
    p1 = q1 * f0 + q0 * f1 ; 
    if ( c->verbose > 4) {
      printf ( "n = %d\n" , n ) ;
      printf ( "q: %10.4g %10.4g\n" , q0 , q1 ) ; 
      printf ( "f: %10.4g %10.4g\n" , f0 , f1 ) ; 
      printf ( "p: %10.4g %10.4g\n" , p0 , p1 ) ; 
      printf ( "b: %10.4g %10.4g\n" , bias[n][0] , bias[n][1] ) ; 
      printf ( "\n" ) ; 
      pause_for_return() ; 
    }
    f0 = p0 * bias[n][0] ;    f1 = p1 * bias[n][1] ;  s=f0+f1 ;
    if ( s < BOT || ( s > TOP ) ) { s = 1.0/s ; f0 *= s ; f1 *= s ; }
    f[n][0] = f0 ;    f[n][1] = f1 ;
    /*
    f[n][1] = bias[n][1] * ( f[n-1][0] * q[n][1] + f[n-1][1] * q[n][0]  );
    f[n][0] = bias[n][0] * ( f[n-1][0] * q[n][0] + f[n-1][1] * q[n][1]  );

    b[n][1] = bias[n][1] * ( b[n+1][0] * q[n+1][1] + b[n+1][1] * q[n+1][0]  );
    b[n][0] = bias[n][0] * ( b[n+1][0] * q[n+1][0] + b[n+1][1] * q[n+1][1]  );
    */

  } 
  if ( c->verbose > 3) {
    printf ( "state after forward pass:\n%10s %10s\n" , "f[n][0]  " , "f[n][1]  " ) ;
    for ( n = 1 ; n <= N ; n ++ ) {
      printf ( "%-10.4g %-10.4g %2d\n" , f[n][0] , f[n][1] , n) ; 
    }
    printf ( "\n" ) ; 
    pause_for_return() ; 
  }
  f0 = f[N][0] ;   f1 = f[N][1] ; 
  p0 = 1.0     ;   p1 = 1.0     ;
  for ( n = N ; n >= 1 ; n-- )  {
    if ( n < N ) {
      q0 = q[n+1][0] ;    q1 = q[n+1][1] ;
      p0 = q0 * b0 + q1 * b1 ;     p1 = q1 * b0 + q0 * b1 ; 
    }
    b0 = p0 * bias[n][0] ;       b1 = p1 * bias[n][1] ;  s=b0+b1 ;
    if ( s < BOT || ( s > TOP ) ) { s = 1.0/s ; b0 *= s ; b1 *= s ; }

    /* we can now figure out the most probable state of trellis */
    p_t0 = f0 * p0 ;     p_t1 = f1 * p1 ; 
    to[n] = ( p_t1 > p_t0 ) ;

    /* now to the likelihoods */
    f0 = f[n-1][0] ;    f1 = f[n-1][1] ;

    /*    b[n][0] = b0 ;    b[n][1] = b1 ; */
    r[n][0] = f0*b0 + f1*b1 ; 
    r[n][1] = f1*b0 + f0*b1 ;  /* beware overflow here */
  }

  if ( c->verbose > 3) {
    printf ( "likelihoods:\n%10s %10s\n" , "r[n][0]  " , "r[n][1]  " ) ;
    for ( n = 1 ; n <= N ; n ++ ) {
      printf ( "%-10.4g %-10.4g %2d\n" , r[n][0] , r[n][1] , n ) ; 
    }
    printf ( "\n" ) ; 
    pause_for_return() ; 
  }

  /* pseudo-vertical pass */
  viols = 0 ; 
  for ( m = 1 ; m <= c->K ; m++ ) {   /* run through source bits */
    umax = a->num_mlist[m] ;         /* number of repetitions      */
    mlist = a->mlist[m] ;            /* list of bits               */
    f0 = 1.0 ; f1 = 1.0 ;            /* should insert prior here if source 
					not dense */
    v0 = 0 ; v1 = 0 ; /* count votes for 0 and 1 */
    for ( u = 1 ; u <= umax ; u ++ ) {
      n = mlist[u] ;                  
      f0 *= r[n][0] ;       f1 *= r[n][1] ;   /* beware overflow here */
      if (to[n-1]^to[n]) v1 ++ ; else v0 ++ ;
      /* if ( r[n][0] > r[n][1] ) v0 ++ ; else v1 ++ ; */ /* detect inconsistent state */
    }
    if ( v0 * v1 ) { /* then this bit not yet consistent */