ra.c

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

      viols ++ ; 
      if ( c->verbose > 3) {
	printf ( "%dv" , m ) ; 
      }
    }
    if ( f1 > f0 ) { so[m] = 1 ;  }
    else { so[m] = 0 ; }
    if ( so[m] ^ ( v1 > v0 ) ) { /* disagreement between mp bit and votes */
      viols ++ ;
      if ( c->verbose > 5 ) {
	printf ( "so[m] = %d but v0 = %d and v1 = %d. (%f,%f)\n" ,
		 so[m] , v0 , v1 , f0 , f1 ) ; 
      }
    }
    for ( u = umax ; u >= 1 ; u -- ) {
      n = mlist[u] ;                  
      p0 = f0 / r[n][0] ;       p1 = f1 / r[n][1] ;  s = 1.0/(p0 + p1) ;  
      q[n][0] = p0 * s ;       q[n][1] = p1 * s ; 
      if ( c->doclip ) { /* pragmatic clipping procedure to stop overflow */
	if ( q[n][0] > c->clip ) { 
	  q[n][0] = c->clip ;
	  q[n][1] = 1.0 - c->clip ;
	  if ( !clipwarn ) { fprintf ( stderr , "c" ) ; fflush ( stderr ) ; clipwarn=1; }
	}
	else           if ( q[n][1] > c->clip ) { 
	  q[n][1] = c->clip ;
	  q[n][0] = 1.0 - c->clip ;
	  if ( !clipwarn ) { fprintf ( stderr , "c" ) ; fflush ( stderr ) ; clipwarn=1; }
	}
      }

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

static double bern ( int errors , int trials , double *p_point , double *p_upper
		   ,  double *p_lower , double zscore ) {
  double factor ; double lowerbound = 1e-6 ;
  if ( errors == 0  ) {
    *p_point = lowerbound ;
    *p_lower = lowerbound ;
    *p_upper = 1.0 - exp( - zscore / (double)(trials) ) ; 
    factor = 100.0 ; 
  } else {
    *p_point =  (double)(errors)/(double)(trials) ;
    /* factor corresponding z sd's in log p */
    factor = exp ( zscore * sqrt( (double)(trials - errors ) / ((double)(errors * trials ) ) )) ;
    *p_upper = *p_point * factor ; 
    *p_lower = *p_point / factor ; 
  }
  return factor ;
}

static void finalline ( FILE *fp , RA_control *c , int totals )
{
  double factor ; 
  if ( totals ) {
    factor = bern ( c->block_errs , c->blocks , &(c->bep) , &(c->bep_u) , &(c->bep_l) , c->zscore  ) ;
    c->bitep = (double) c->bit_errs / (double)( c->K * c->blocks ) ;
    c->bitep_u = c->bitep * factor ; 
    c->bitep_l = c->bitep / factor ; 
    c->bep_det   = (double) c->block_det   / (double)( c->blocks ) ;
    c->bep_detlw = (double) c->block_detlw / (double)( c->blocks ) ;
    c->bep_undet = (double) c->block_undet / (double)( c->blocks ) ;
  }

  if ( (!(c->outappend)) || (totals) ||
       ( c->pheading_period &&  
	 ! ( ( ++ c->pheading_count ) % c->pheading_period ))  ) {
    if ( !totals ) {
      fprintf ( fp , "#es  viols t_es ") ; 
    }
    fprintf ( fp , "#blks es un det valid " ) ; 
    fprintf ( fp , "Eb/No(dB)    x    xass ") ; 
    if ( !totals ) {
      fprintf ( fp , "   xact flipd ") ; 
    }
    if ( totals ) {
      fprintf ( fp , "  K    N ") ; 
      fprintf ( fp , "biters undet det ") ; 
      fprintf ( fp , "block_ep u     l    ") ; 
      fprintf ( fp , "bit_ep   u     l    ") ; 
      fprintf ( fp , "blep_det undet ") ; 
      fprintf ( fp , "detlw blep_detlw " ) ; 
      fprintf ( fp , "totloops lmax " ) ; 
    }
    if ( !totals ) {
      fprintf ( fp , "loops ") ; 
    }
    fprintf ( fp , "\n" ) ; 
  }

  if ( !totals ) {
    fprintf ( fp , "%4d " , c->count ) ; /* bits wrong this block */
    fprintf ( fp , "%4d " , c->viols ) ; /* number of detected problems this block */
    fprintf ( fp , "%4d " , c->tcount ) ; /* bits wrong this block */
  }

  fprintf ( fp , "%4d " , c->blocks ) ;  /* blocks transmitted */
  fprintf ( fp , "%3d " , c->block_errs ) ;  /* block errors */
  fprintf ( fp , "%1d " , c->block_undet ) ; /* undetected */
  fprintf ( fp , "%3d " , c->block_det ) ;   /* detected */
  fprintf ( fp , "%4d " , c->block_valid ) ;  /* blocks that reached a valid decode state */

  fprintf ( fp , "%7.4g " , c->ebno ) ;  /* gaussian channel EbNo */

  fprintf ( fp , "%7.4g " , c->gcx ) ;    /* gaussian channel mean x */
  fprintf ( fp , "%7.4g " , c->gcxass ) ; /* assumed x */
  if ( !totals ) {
    fprintf ( fp , "%7.4g " , c->gcxact ) ; /* actual x this time 
					      (including fluctuation) */
    fprintf ( fp , "%4d " , c->flipped ) ; /* number of t bits "flipped" (a measure of noise level) */
  }

  if ( totals ) {

    fprintf ( fp , "%4d " , c->K ) ; /* K source bits  */
    fprintf ( fp , "%4d " , c->N ) ; /* N transmitted bits */

    fprintf ( fp , "%4d " , c->bit_errs ) ;  /* bit errors */
    fprintf ( fp , "%4d " , c->bit_undet ) ; /* undetected */
    fprintf ( fp , "%4d " , c->bit_det ) ;   /* detected */

    fprintf ( fp , "%8.4g " , c->bep ) ;  /* block error probability */
    fprintf ( fp , "%7.2g " , c->bep_u ) ;  /* block error probability (upper) */
    fprintf ( fp , "%7.2g " , c->bep_l ) ;  /* block error probability (lower) */

    fprintf ( fp , "%8.4g " , c->bitep ) ;  /* bit error probability */
    fprintf ( fp , "%7.2g " , c->bitep_u ) ;  /*  error probability (upper) */
    fprintf ( fp , "%7.2g " , c->bitep_l ) ;  /*  error probability (lower) */

    fprintf ( fp , "%9.3g " , c->bep_det ) ;  /* detected block error probability */
    fprintf ( fp , "%9.3g " , c->bep_undet ) ;  /* undetected block error probability  */

    fprintf ( fp , "%3d " , c->block_detlw ) ;   /* detected errors involving low weight error */
    fprintf ( fp , "%9.3g " , c->bep_detlw ) ;  /* detected block error probability involving low weight error words */
    fprintf ( fp , "%4d " , c->totloops ) ; /* totloops used in finding valid decodings */
    fprintf ( fp , "%4d " , c->loops ) ;  /* when we would have stopped */
  } else {
    fprintf ( fp , "%4d " , c->loop ) ;   /* iterations this time */
  }

  fprintf ( fp , "\n" ) ; 

}

static int make_sense ( RA_control *c ) 
{ /*  correct silly RA_control parameters */
  int status = 0 ; 
  
  /* if no assumed noise level stated, set assumed equal to true */
  if ( c->gcxass < 0.0 ) { c->gcxass = c->gcx ; c->gcxact = c->gcx ; }
  if ( c->fnass < 0.0 ) { c->fnass = c->fn ; }

  c->R = (double) (c->K) / (double) (c->N) ; 
  c->ebno = 10.0 * log ( c->gcx * c->gcx / ( 2.0 * c->R )  ) / log ( 10.0 ) ; 

  return status ; 
}

static int make_space ( RA_control *c ) {
/* see also mnc_free */
  int status = 0 ; 

  c->so = cvector ( 1 , c->K ) ;  
  c->s = cvector ( 1 , c->K ) ;  
  c->t = cvector ( 1 , c->N ) ;
  c->to = cvector ( 0 , c->N ) ;
  c->to[0] = 0 ;

  c->histo = ivector ( 0 , c->loops+1 ) ; /* histogram of number of loops taken  */
  set_ivector_const ( c->histo ,  0 , c->loops+1 , 0 ) ; 

  /*  vec->bias = dvector ( 1 , vec->N ) ; */

  c->q    = dmatrix ( 1 , c->N , 0 , 1 ) ;
  c->bias = dmatrix ( 1 , c->N , 0 , 1 ) ;
  c->f    = dmatrix ( 0 , c->N , 0 , 1 ) ;
  c->r    = c->f ; /* they can share space ! */
  c->f[0][0] = 1 ; c->f[0][1] = 0 ; 

  c->block_valid = 0 ; 
  c->blocks=0 ; 
  c->block_det = 0 ; c->block_undet = 0 ; c->block_errs = 0 ; 
  c->bit_det = 0 ; c->bit_undet = 0 ; c->bit_errs = 0 ; 
  c->block_detlw = 0 ; 
  c->lowweight = c->N / 10 ;
  /* AD HOC definition of a 'low weight detected error event' */

  return status ; 
}

int check_alist_MN ( alist_matrix *a , RA_control *c ) {
  int status = 0 ; 
  if (     c->N != a->N 
      ||   c->K != a->M      ) {
    fprintf ( stderr , "setting N and K from alist: %d %d %d %d\n" , c->N , a->N , c->K , a->M ) ; 
    c->N = a->N ;
    c->K = a->M ; 
  }
  return status ; 
}

static int score ( RA_control *c ) 
/* compares so with s and works out quality of solution if wrong. */
{
  int status = 0 ; 
  int n , k , count = 0 ; 
  FILE *fp ; 

  if ( c->verbose ) printf ( "Scoring answer --- \n" ) ; 
  if ( c->verbose > 2) {
    printf ( "original source vector:\n" ) ;
    for ( k = 1 ; k <= c->K ; k ++ ) {
      if ( c->s[k] ) printf ("1 ") ; else printf ( "0 " ); 
    }
    printf ( "\n" ) ; 
  }
  if ( c->verbose > 2) {
    printf ( "reconstructed source vector:\n" ) ;
    for ( k = 1 ; k <= c->K ; k ++ ) {
      if ( c->so[k] ) printf ("1 ") ; else printf ( "0 " ); 
    }
    printf ( "\n" ) ; 
  }

  c->blocks ++ ; 
  if ( !(c->viols) ) {
    c->block_valid ++ ;
    c->totloops += c->loop ; 
    c->histo[c->loop] ++ ;
  }
  c->count = 0 ; 
  for ( n = 1 ; n <= c->K ; n++ ) { /* this compares all bits */
    if ( c->so[n] != c->s[n] ) {    
      c->count ++ ;  /* count is the errors */
    }
  }
  c->tcount = 0 ;
  /* count the weight of the difference in codeword land */
  for ( n = 1 ; n <= c->N ; n++ ) { 
    if ( c->to[n] != c->t[n] ) {    
      c->tcount ++ ;  
    }
  }
  if ( c->count > 0 ) { 
    c->block_errs ++ ; c->failcount ++ ; /* why two registers!? */
    c->bit_errs += c->count ;
    if ( c->viols ) {
      c->block_det ++ ;
      c->bit_det += c->count ;
      /* check for near-codeword */
      if ( c->tcount < c->lowweight ) {
	c->block_detlw ++ ;
      }
    } else {
      printf ( "UNDETECTED ERROR\n" ) ; 
      c->block_undet ++ ;
      c->bit_undet += c->count ;
    } 
    if ( c->count && c->maxcount != 0 && c->error_log ) {
      fp = fopen ( c->error_logfile , "a" ) ;
      if ( !fp ) {
	fprintf ( stderr , " couldn't open logfile %s\n" , c->error_logfile ) ; 
	c->error_log = 0 ; 
      } else {
	/* tell me more about this failure */ 	    /* write to a file */
	finalline ( fp , c , 0 ) ;
	fprintf ( fp , "# " ) ;
	for ( n = 1 ; n <= c->K && ( ( count <= c->maxcount ) 
				    || c->maxcount < 0 ) ; n++ ) { 
	  if ( c->so[n] != c->s[n] ) {    
	    count ++ ;
	    if ( c->so[n] ) {
	      fprintf ( fp , "+%4d " , n ) ; /* + indicates extra bit given in so */
	    } else {
	      fprintf ( fp , "-%4d " , n ) ; /* - indicates bit lacking     */
	    }
	  }
	}
	fprintf ( fp , "\n" ) ;
	fclose (fp ) ; 
      }
    }
  }  else {
    if ( c->verbose ) printf ( "SUCCESS\n" ) ; 
  }
  return status ; 
}

static int t_to_b (  unsigned char *x ,	RA_control *c ) {
  int i , co = 0 ;
  int lo = 1 ; int hi = c->N ; 
  double z , p ;
  double gcx = c->gcx ; 
  double gcxass = c->gcxass ; 
  double **b = c->bias ; 

  if ( c->gc ) {
    for ( i = lo ; i <= hi ; i ++ ) {
      /* make a random normal variate with s.d. 1.0 and mean +/- gcx */
      z = ( x[i] ? gcx : - gcx ) + rann() ; 
      p = 1.0 / ( 1.0 + exp ( - 2.0 * z * gcxass ) ) ; /* likely to be large if x = 0 */
      /*   was   b[i] = 1.0 - p ;  */
      b[i][1] =  p ;      b[i][0] = 1.0 - p ;
      if ( !(( p < 0.5 )^(x[i])) ) co++ ; /* counts as a flipped bit */
    }

  } else { /* BSC  
    co = random_cvector ( c->x , c->fn , lo, hi ) ; 
    for ( i = lo ; i <= hi ; i ++) {
      if ( (c->x[i]) ^ (t[i]) ) {
	c->bias[i] = 1.0 - c->fnass ;
      } else {
	c->bias[i] =  c->fnass ;
      }
    } */
  }
    
  return co ;
}

static void dc_defaults ( RA_control *c ) 
{
#include "RAdc_var_def.c"
}

static int process_command ( int argc , char **argv , RA_control *c ) {
  int p_usage = 0 ;
  int status = 0 ;
  int cs , i ;

  if ( argc < 1 )     {
    p_usage = 1 ; 
    status -- ;
  }

#define ERROR1 fprintf ( stderr , "arg to `%s' missing\n" , argv[i] ) ; \
               status --
#define ERROR2 fprintf ( stderr , "args to `%s' missing\n" , argv[i] ) ; \
               status --
  for (i = 1 ; i < argc; i++)    {
    cs = 1 ;
    if ( strcmp (argv[i], "-V") == 0 )        {
      c->verbose = 1;
    }
    else  if ( strcmp (argv[i], "-VV") == 0 )        {
      c->verbose = 2;
    }
#include "RAdc_var_clr.c"
    else {
      fprintf ( stderr , "arg `%s' not recognised\n" , argv[i] ) ; 
      p_usage = 1 ;
      status -- ;
    }
    if ( cs == 0 ) {
      fprintf ( stderr , "arg at or before `%s' has incorrect format\n" , 
	       argv[i] ) ;
      p_usage = 1 ;
      status -- ;
    }
  }
  if ( p_usage ) print_usage ( argv , stderr , c ) ;
  return ( status ) ;
}
#undef ERROR1
#undef ERROR2

#define DNT fprintf( fp, "\n        ")
#define NLNE  fprintf( fp, "\n")

static void print_usage ( char **argv , FILE * fp , RA_control *c )
{
  fprintf( fp, "Usage: %s ",argv[0]);
  fprintf( fp, " [optional arguments]");

  DNT;  fprintf( fp, "-V | -VV                 (verbose or very verbose)"); 
  NLNE; fprintf( fp, " Data creation:" ) ; 
#include "RAdc_var_usg.c"
  pause_for_return();
  fprintf( fp, "\n");
  return ;
}

#undef DNT
#undef NLNE

static void RA_free ( RA_control *vec ) {

  free_cvector ( vec->so, 1 , vec->K ) ; 
  free_cvector ( vec->t , 1 , vec->N ) ; 

  free_cvector ( vec->s , 1 , vec->K ) ; 

  free_dmatrix ( vec->bias , 1 , vec->N , 0 , 1 ) ; 
  free_dmatrix ( vec->q , 1 , vec->N , 0 , 1 ) ; 
  free_dmatrix ( vec->f , 0 , vec->N , 0 , 1  ) ; 

}

double h2 ( double x ) {
  double tmp ; 
  if ( x <= 0.0 || x>= 1.0 ) tmp = 0.0 ;
  tmp = x * log ( x ) + ( 1.0 - x ) * log ( 1.0 - x ) ;
  
  return - tmp / log ( 2.0) ; 
}

/*
<!-- hhmts start -->
<!-- hhmts end -->
*/