gallagers.c

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

      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 " , v->count ) ; /* bits wrong this block (have option to put source bits here) */
    fprintf ( fp , "%4d " , v->viols ) ; /* number of detected problems this block */
    fprintf ( fp , "%4d " , v->count ) ; /* bits wrong this block (transmitted) */
  }

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

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

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

  if ( totals ) {

    fprintf ( fp , "%4d " , (v->N)-(dc->trueM) ) ; /* K source bits  */
    fprintf ( fp , "%4d " , (v->N) ) ; /* N transmitted bits, NO, this is the number of bits in the graph. The transmitted number is "nn"; this has been taken into account in the EbNo computation */

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

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

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

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

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

  fprintf ( fp , "\n" ) ; 
}

static int make_sense ( data_creation_param *dc , mnc_all *all ) 
{ /*  correct silly control parameters */
/* 
   first the data creation 
                             */
  int status = 0 ; 
  char junk[200] ; 
  int N ; 

  if ( dc->M == 0 ) {
    dc->M = dc->N ; /* default is a square MNC */
  }
                                 /* there is potential confusion 
				    here:
				    the "n,k" of the error correcting code
				    are called "M,N" in the dc-> structure.
				    But when it comes to decoding, 
				    the matrix is A is M * N
				    where M = dc->M (traditionally known as n)
				    and   N = dc->N + dc->M = "k + n" 
				    
				    this comment needs to go to check_alist also */
  /* here are the standard settings if Gallanger code */  
  all->vec->M = dc->M ; 
  N = dc->M + dc->N ; 
  all->vec->N = N ; 
  all->vec->NS = 0 ; 
  all->vec->NN = N ; 
  all->vec->nsfrom = 1 ; 
  all->vec->nsto   = 1 ; 
  all->vec->nnfrom = 1 ; 
  all->vec->nnto   = N ;
  /* but here are the modifications caused by dc->ns */
  if ( dc->ns ) {
    all->vec->NS = dc->ns ; 
    all->vec->NN = N - dc->ns ;
    if (dc->nsfrom == 1) {
      all->vec->nsfrom = 1 ;
      all->vec->nsto   = dc->ns ;
      all->vec->nnfrom = dc->ns + 1 ; 
      all->vec->nnto   = N ;
    } else {
      all->vec->nsfrom = N - dc->ns + 1 ;
      all->vec->nnfrom = 1 ; 
      all->vec->nnto   = N - dc->ns ;
      all->vec->nsto   = N ;
    }
  }
  dc->nn = all->vec->NN ; /* NB this overwrites anything in nn */
  
/* default name for afile */

  sprintf ( junk , "MNC%d/%d.%d.%ld" , 
	   dc->MNC , dc->N , dc->tc , dc->mseed ) ; 
  if ( dc->reada==1 )  /* if this is 2 then filename
			  's been put there already */
    sprintf ( dc->afile , "%s" ,  junk ) ; 
  if ( dc->reada==0 )  
    fprintf ( stderr , "oi, how can reada be zero?\n" ) ; 

/* 
   and now the free energy control
                                       */
  if ( dc->use_fe ) {
    status += make_sense_c ( dc , all->c ) ; 
    if ( all->c->second ) 
      status += make_sense_c ( dc , all->c2 ) ; 
    if ( all->c->third ) 
      status += make_sense_c ( dc , all->c3 ) ; 
  }

  return status ; 
}


static int make_sense_c (  data_creation_param *dc , fe_min_control *c ) {
  int status = 0 ; 

  if ( ( dc->MNC ==3 || dc->MNC == 4 || c->NL > 1 )  
      && ( !(c->opt==5 || c->opt==6) )
      && ( c->betastyle == 0 ) )
    fprintf ( stderr , "Warning: no annealing schedule given, MNC = %d, NL = %d\n" , 
	     dc->MNC , c->NL ) ; 
  if ( ( c->NL > 1 || c->betastyle != 0 ) && ( (c->opt==5) ) )
    fprintf ( stderr , "Warning: this optimizer doesn't need annealing or multiple loops, opt = %d, NL = %d\n" , 
	     c->opt , c->NL ) ; 

  if ( c->betastyle == 1 && c->NL > 1 ) 	
    c->betaf = ( c->beta1 - c->beta0 ) / ((double) ( c->NL - 1) ) ; 
  else if ( ( c->betastyle == 2 || c->betastyle == 22 ) && c->NL > 1 )
    c->betaf = exp ( log ( c->beta1 / c->beta0 ) / ((double) ( c->NL - 1) )) ; 

  if ( c->opt > 6 || c->opt < 0 ) {
      fprintf ( stderr , " opt %d not existo\n" , c->opt ) ;   status-- ; }

  if ( c->opt == 5 || c->opt == 6 ) {
    c->zero_temperature = 1 ; /* changed from mnc5 on */
  } 

  c->dgamma = (double) (c->gamma) ;

/* any optimizer-specific stuff? */
  switch ( c->opt ) {
  case ( 0 ) : case ( 1 ) : case(2): case(3): case (4): default: 
    break ; 
      
  case(5):    case(6):
    c->macarg.itmax = c->itmax ; 
    break ;
  }

  return status ; 
}

static int make_space ( mnc_all *all , mnc_vectors *vec ) {
/* see also mnc_free */
  int status = 0 ; 

  vec->xo = cvector ( 1 , vec->N ) ; 
  vec->t = cvector ( 1 , vec->M ) ;

  vec->x = cvector ( 1 , vec->N ) ; 
  vec->z = cvector ( 1 , vec->M ) ;
  vec->y = cvector ( 1 , vec->M ) ; 

  vec->bias = dvector ( 1 , vec->N ) ;
  vec->n = dvector ( 1 , vec->N ) ;
  vec->nsum = ivector ( -1 , vec->N ) ;
  set_ivector_const ( vec->nsum ,  -1 , vec->N , 0 ) ;

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

  return status ; 
}

static int hook_mnc_vec_to_fe ( fe_min_param *p , mnc_vectors *vec ) {
  int status = 0 ; 
  p->M = vec->M ;
  p->N = vec->N ;
  p->z = vec->z ; 
  p->s = vec->x ; 
  /* where the true vector lives */
  p->so = vec->xo ; 
  /* where the output goes */
  return status ; 
}

static int hook_mnc_vec_to_bnd ( bnd_param *p , mnc_vectors *vec ) {
  int status = 0 ; 
  p->M = vec->M ;
  p->N = vec->N ;
  p->z = vec->z ; 
  /* where the true vector lives */
  p->x = vec->x ; 
  p->bias = vec->bias ;
  /* where the output lives */
  p->xo = vec->xo ; 
  return status ; 
}

int check_alist_MN ( alist_matrix *a , mnc_vectors *v , data_creation_param *dc ) {
  int status = 0 ; 
  if (     v->N != a->N 
      ||   v->M != a->M      ) {
    fprintf ( stderr , "setting N and M from alist: %d %d %d %d\n" , v->N , a->N , v->M , a->M ) ; 
    v->N = a->N ;
    v->M = a->M ; 
  }
  if ( !(dc->trueM) ) dc->trueM = a->M ; /* assume that the number of constraints equals M. */
  /* rate is number of free bits, which is total length of matrix minus trueM,
     divided by the number of transmitted bits, which is NN */
  dc->Rate = (double) (a->N - dc->trueM) / (double) (dc->nn) ; 
  if ( dc->gc ) {
    printf ( "R = %9.4g, x = %9.4g\n" , dc->Rate , dc->gcx ) ; 
    dc->ebno = 10.0 * log ( dc->gcx * dc->gcx / ( 2.0 * dc->Rate )  ) / log ( 10.0 ) ; 
    printf ( "Eb/No = %9.4g\n" , dc->ebno ) ; 
  } else {
    dc->ebno = -10.0 ; 
  }
  if (dc->lowweight <= 0 ) {
    dc->lowweight = a->N / 10 ;
  }

  return status ; 
}


static int score ( data_creation_param *dc , 
       		   mnc_vectors  *v , mnc_all *all 
) /* compares so with s and works out quality of solution if wrong. */
{
  int status = 0 ; 
  int n , count = 0 , m , M = v->M , *nsum = v->nsum ; 
  FILE *fp ; 

  if ( dc->verbose ) printf ( "Scoring answer --- " ) ; 
  v->count = 0 ; v->count_high = 0 ; v->viols = 0  ;
  dc->blocks ++ ;
  for ( n = 1 ; n <= v->N ; n++ ) { /* this compares all bits */
    if ( v->xo[n] != v->x[n] ) {    
      v->count ++ ; 
    }
    v->count_high += v->xo[n] ;
  }
  if ( v->count > 0 ) {
    alist_times_cvector_mod2 ( (all->bndp->a) , v->xo , v->t ) ; 

    for ( m = 1 ; m <= M ; m ++ ) {
      if ( v->t[m] != v->z[m] ) v->viols ++ ; 
    }
    if ( dc->verbose ) printf ( "fail\n" ) ; 
/*    else fprintf ( stderr , " X" ) ;  */
    dc->failcount ++ ; 

    dc->bit_errs += v->count ;
    if ( v->viols ) {
      dc->block_det ++ ;
      dc->bit_det += v->count ;
      /* check for near-codeword */
      if ( v->count < dc->lowweight ) {
	dc->block_detlw ++ ;
      }
    } else {
      printf ( "UNDETECTED ERROR\n" ) ; 
      dc->block_undet ++ ;
      dc->bit_undet += v->count ;
    } 

    /*    if ( v->count ) { strange 'if' in mncN */
    if ( dc->maxcount != 0 && dc->error_log ) {
      if ( dc->verbose > 3 )
	fprintf ( stdout , "logging error to %s\n" , dc->error_logfile ) ; 
      fp = fopen ( dc->error_logfile , "a" ) ;
      if ( !fp ) {
	fprintf ( stderr , " couldn't open logfile %s\n" , dc->error_logfile ) ; 
	dc->error_log = 0 ; 
      } else {
	/* tell me more about this failure */ 	    /* write to a file */
	finalline ( fp , all , dc , 0 ) ;
	fprintf ( fp , "# " ) ;
	for ( n = 1 ; n <= v->N && ( ( count <= dc->maxcount ) 
				     || dc->maxcount < 0 ) ; n++ ) { 
	  if ( v->xo[n] != v->x[n] ) {    
	    count ++ ;
	    if ( v->xo[n] ) {
	      fprintf ( fp , "+%4d " , n ) ; /* + indicates extra bit given in so */
	    } else {
	      fprintf ( fp , "-%4d " , n ) ; /* - indicates bit lacking     */
	    }
	  }
	}
	fprintf ( fp , "\n" ) ;
	fclose (fp ) ; 
      }
      if ( dc->noise_sum ) { 
/* Keep track of the noise vectors that caused the errors (mean only) */
	fp = fopen ( dc->noise_sumfile , "w" ) ;
	if ( !fp ) {
	  fprintf ( stderr , " couldn't open logfile %s\n" , 
		   dc->noise_sumfile ) ; 
	  dc->noise_sum = 0 ; 
	} else { /* pop the total number of failures and number of bits wrong
		     at the top in nsum[0], nsum[-1] */
	  nsum[0] ++ ;
	  for ( n = 1 ; n <= v->N ; n++ ) { 
	    if ( v->x[n] ) {    
	      nsum[n] ++ ; nsum[-1] ++ ;
	    }
	  }
	  fprintf ( fp , "#noise %d blocks %d bits %d block_error_rate %d / %d \n" , nsum[0] , nsum[-1] , v->N  , dc->failcount, dc->message ) ;
	  for ( n = 1 ; n <= v->N ; n++ ) { 
	    fprintf ( fp , "%d" , nsum[n] ) ;
	    fprintf ( fp , "\n" ) ;
	  }
	  fclose (fp ) ; 
	}
      }
    }
  }  else {
    if ( dc->verbose ) printf ( "SUCCESS\n" ) ; 
  }
  if ( !(v->viols) ) {
    dc->block_valid ++ ;
    dc->totloops += all->bndc->loop ; 
    v->histo[all->bndc->loop] ++ ;
  }

  /* could also work out the typicality of the solution */
  return status ; 
}

static int make_norder ( alist_matrix *alist , fe_min_control *c )
{ /*
     Once this is changed to include other orders, 
     then the command must be called every time there is a 
     switch of c (or some other solution to this problem found)
*/
  int status = 0 ; 
  int n , N=alist->N ; 

    /* set up the norder */
  switch ( c->NORDER ) {
  case(0):
    
    for ( n = 1 ; n <= N ; n++ ) {
      alist->norder[n] = n ; 
    }
    break;
  default:
    status -- ;
    fprintf ( stderr , "invalid NORDER %d \n" , c->NORDER ) ; 
    break ;