gallagers.c

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

  }

  return status ; 
}

/* needs personalising for fe or depersonalising for bnd */
int make_vectors_quick ( data_creation_param *dc , mnc_vectors *v , mnc_all *all )
{
  int status = 0 ; 
  int count_s = 0 , count_n = 0 ; 

  /* the two main options used in 1998 are
     dc->gc = 1 (ordinary gaussian channel )
     and 
     dc->gc = 0 (BSC)

     other minor options are: dc->nsn  fixed weight noise ; dc->fgc  fixed gaussian noise power.
     */
  if ( dc->nsn ) { /* a fixed weight vector for all s and n is required */
    fixed_wt_cvector ( v->x , dc->nsn , v->nsfrom , v->nnto ) ;
    v->count_s = dc->nsn ; 
    dc->gcx =  (double)  dc->nsn / (double) v->N ;
    dc->gcxact = dc->gcx ;
  } else { /* other data generation rules */
    if ( !dc->ns ) { /* in NMN codes there are source bits to make */ 
      /* normally ns=0 , and nsfrom=1 and nsto=1 */
      count_s = random_cvector ( v->x , dc->fs , v->nsfrom , v->nsto ) ; 
    } else {
      fixed_wt_cvector (  v->x , dc->ns , v->nsfrom , v->nsto ) ;
      count_s = dc->ns ; 
    }
    
    if ( dc->fgc ) { /* fixed noise level first, assumed noise level second */
      count_n = make_fgaussian_noise_bits_and_fix_biases 
	( v->x , v->bias , v->n , dc->fgcx, dc->fgcnf , dc->gcx , v->nnfrom , v->nnto ) ;
    }
    else if ( dc->gc ) { /* Here is assumed that using bnd */
      count_n = make_gaussian_noise_bits_and_fix_biases 
	( v->x , v->bias , dc->gcx , v->nnfrom , v->nnto ) ;
      dc->fgcx = dc->gcx ; 
      dc->gcxact = dc->gcx ;  /* ideally this should be set to the 
				 actual noise level each time.... */
    } else { /* ordinary BSC stuff */
      if ( !dc->nn ) {
	count_n = random_cvector ( v->x , dc->fn , v->nnfrom , v->nnto ) ; 
	dc->gcx = dc->fn ; /* write fn into gcx */
      } else { /* this makes the noise fixed weight, independent of what the {s} would be (for MN code) */
	fixed_wt_cvector (  v->x , dc->nn , v->nnfrom , v->nnto ) ;
	count_n = dc->nn ; /* the desired weight */
	dc->gcx =  (double)  dc->nn / (double) v->N ;
      }
      dc->gcxact =  (double) count_n / (double) v->N ;
    }
    v->count_s = count_s + count_n ;
  }
    
  dc->true_fs = 0.0 ;
  dc->true_fn = (double) count_n / (double) dc->M ; 
    
  alist_times_cvector_mod2 ( all->a , v->x , v->z ) ; 

  return status ; 
}
int make_vectors_quickER ( data_creation_param *dc , mnc_vectors *v , mnc_all *all ) /* just makes the calls to the rand number generator */
{
  int status = 0 ; 
  int count_s = 0 , count_n = 0 ; 

  if ( !dc->ns ) {
    count_s = random_cvector ( v->x , dc->fs , v->nsfrom , v->nsto ) ; 
  } else {
    fixed_wt_cvector (  v->x , dc->ns , v->nsfrom , v->nsto ) ;
  }

  if ( dc->gc ) { /* Here is assumed that using bnd */
    count_n = make_gaussian_noise_bits_ONLY 
      ( v->x , v->bias , dc->gcx , v->nnfrom , v->nnto ) ;
  }
  else {
    if ( !dc->nn ) {
      count_n = random_cvector ( v->x , dc->fn , v->nnfrom , v->nnto ) ; 
    } else {
      fixed_wt_cvector (  v->x , dc->nn , v->nnfrom , v->nnto ) ;
    }
  }

  return status ; 
}

static void process_vector_fe ( mnc_vectors *v , fe_min_param *p ) {
  int m ; 
  p->true_s = 1 ; 
  for ( m = 1 ; m <= p->M ; m++ ) {
    p->g[m] = ( v->z[m] ) ? 1.0 : -1.0 ; 
  }
}

static int make_gaussian_noise_bits_and_fix_biases ( unsigned char *x ,
						    double *b,
						    double gcx ,
						    int lo , int hi ) {
  int i , c = 0 ;
  double z , p ;

  for ( i = lo ; i <= hi ; i ++ ) {
    /* make a random normal variate with s.d. 1.0 and mean gcx */
    z = gcx + rann() ; 
    p = 1.0 / ( 1.0 + exp ( - 2.0 * z * gcx ) ) ; 
    if ( z <= 0.0 ) { /* a noise bit is high */
      x[i] = 1 ; c ++ ;
      b[i] = p ; 
      if ( b[i] > 0.5 ) fprintf ( stderr , "something wrong %f\n" , b[i] ) ;
    } else {
      x[i] = 0 ; 
      b[i] = 1.0 - p ; 
      if ( b[i] > 0.5 ) fprintf ( stderr , "something wrong %d %f\n" , x[i] , b[i] ) ;
    }
  }
  return c ; 
}

/* fixed noise - by making n a unit-like vector */ 
static int make_fgaussian_noise_bits_and_fix_biases ( unsigned char *x ,
						    double *b,double *n,
						    double fgcx ,
						    double fgcnf ,
						    double gcx ,
						    int lo , int hi ) {
  int i , c = 0 ;
  double z , p , power=0.0 , factor , count = 0.0 ;

  for ( i = lo ; i <= hi ; i ++ ) {
    /* make a random normal variate with s.d. 1.0 and mean gcx */
    n[i] = rann() ; 
    power += n[i] * n[i] ; 
    count += 1.0 ; 
  }
  factor = sqrt( count / power ) ; 
  fprintf ( stdout , "f=%8.6g -> %5g : " , factor , fgcnf ) ; 
  fflush(stdout) ; 
  factor *= fgcnf ; /* boost or reduce noise for importance sampling */
  for ( i = lo ; i <= hi ; i ++ ) {
    /* make a random normal variate with s.d. 1.0 and mean gcx */
    z = fgcx + n[i] * factor ; /* fixed signal to noise level */
    p = 1.0 / ( 1.0 + exp ( - 2.0 * z * gcx ) ) ; /* gcx - assumed value */ 
    if ( z <= 0.0 ) { /* a noise bit is high */
      x[i] = 1 ; c ++ ;
      b[i] = p ; 
      if ( b[i] > 0.5 ) fprintf ( stderr , "something wrong %f\n" , b[i] ) ;
    } else {
      x[i] = 0 ; 
      b[i] = 1.0 - p ; 
      if ( b[i] > 0.5 ) fprintf ( stderr , "something wrong %d %f\n" , x[i] , b[i] ) ;
    }
  }
  return c ; 
}

static int make_gaussian_noise_bits_ONLY ( unsigned char *x ,
						    double *b,
						    double gcx ,
						    int lo , int hi ) {
  int i , c = 0 ;
  double z ;

  for ( i = lo ; i <= hi ; i ++ ) {
    /* make a random normal variate with s.d. 1.0 and mean gcx */
    z = rann() ; 
  }
  return c ; 
}

static void set_up_biases ( double *b , mnc_vectors *v , data_creation_param *dc ) {
  int n ;
  double lfn1fn = log ( dc->fn / ( 1 - dc->fn ) ) ;
  double lfs1fs = log ( dc->fs / ( 1 - dc->fs ) ) ;

  for ( n = v->nsfrom ; n <= v->nsto ; n++ ) {
    b[n] = lfs1fs ;
  }
  for ( n = v->nnfrom ; n <= v->nnto ; n++ ) {
    b[n] = lfn1fn ;
  }
}
  
static void set_up_priors ( double *b , mnc_vectors *v , data_creation_param *dc ) {
  int n ;

  for ( n = v->nsfrom ; n <= v->nsto ; n++ ) {
    b[n] = dc->fs ;
  }
  for ( n = v->nnfrom ; n <= v->nnto ; n++ ) {
    b[n] = dc->fn ;
  }
}
  
static int evaluate_feasibility ( data_creation_param *dc )    
{
  int status = 0 ;
/*
   compute Shannon capacity and compare with info content  
*/

  dc->h2fs =  h2(dc->fs) ;
  dc->h2fn =  h2(dc->fn) ;
  dc->rho =  (double ) dc->N / (double ) dc->M  ;
  dc->rate = dc->rho * dc->h2fs ;
  dc->capacity = 1.0 - dc->h2fn ;
  if ( !dc->gc ) {
    if ( dc->rate >   dc->capacity ) {
      printf ( "task impossible, signed Shannon\n" ) ; 
      if ( dc->notabovecap ) {status -- ; }
    }
    if ( status < 0 || dc->verbose ) 
      printf ( "H(s) = %f, capacity of channel = %f.\n" , 
	      dc->h2fs , dc->capacity ) ; 
  }
  return status ; 
}


static void dc_defaults ( data_creation_param *dc ) 
{
#include "dc_var6_def.c"
}

static int process_command ( int argc , char **argv , mnc_all *all ) {
  data_creation_param *dc = all->dc ; 
  fe_min_control *c = all->c ; 
  fe_min_control *c2 = all->c2 ; 
  fe_min_control *c3 = all->c3 ; 
  bnd_control *bndc = all->bndc ; 

  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 --
#define ERRORREG fprintf ( stderr , "regtype must be defined before `%s'\n" , argv[i] ) ; \
               status --
  for (i = 1 ; i < argc; i++)    {
    cs = 1 ;
    if ( strcmp (argv[i], "-V") == 0 )        {
      c->verbose = 1;
      dc->verbose = 1;
    }
    else  if ( strcmp (argv[i], "-VV") == 0 )        {
      c->verbose = 2;
      dc->verbose = 2;
    }
#include "fe_var6_clr.c"
#include "fed2_var_clr.c"
#include "dc_var6_clr.c"
#include "bnd_var_clr.c"
    else if ( strcmp (argv[i], "-b") == 0 ) {
      if ( i + 3 >= argc ) { ERROR2;      }
      else {
	cs *= sscanf(argv[++i], "%d", &(c->betastyle)); 
	cs *= sscanf(argv[++i], "%lf", &(c->beta0)); 
	cs *= sscanf(argv[++i], "%lf", &(c->beta1)); 
      }
    }
    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 , all ) ;
  return ( status ) ;
}
#undef ERROR1
#undef ERROR2
#undef ERRORREG

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

static void print_usage ( char **argv , FILE * fp ,
			 mnc_all *all 
)
{
  data_creation_param *dc = all->dc ; 
  
  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 "dc_var6_usg.c"
  pause_for_return();
  NLNE; fprintf( fp, " Inference:" ) ; 
  DNT;  fprintf( fp, "-b betastyle beta0 beta1 (what to do with beta)" ); 
  DNT;  fprintf( fp, "       betastyle 0: const; 1: linear; 2: multiply; 22: multiply and go wild on last loop" ); 
  fprintf( fp, "\n");
  fe_print_usage ( argv , fp , all ) ; 
  pause_for_return();
  bnd_print_usage ( argv , fp , all ) ; 
  return ;
}

static void fe_print_usage ( char **argv , FILE * fp ,
			    mnc_all *all 
)
{
  fe_min_control *c = all->c ; 
  fe_min_control *c2 = all->c2 ; 
  fe_min_control *c3 = all->c3 ; 

  NLNE; fprintf( fp, " Further free energy minimization stuff: <defaults>");
#include "fe_var6_usg.c"
#include "fed2_var_usg.c"
  fprintf( fp, "\n");
  return ;
}

static void bnd_print_usage ( char **argv , FILE * fp ,
			    mnc_all *all 
)
{
  bnd_control *bndc = all->bndc ; 

  NLNE; fprintf( fp, " Belief Net decoder:     <defaults>");
#include "bnd_var_usg.c"
  fprintf( fp, "\n");
  return ;
}

#undef DNT
#undef NLNE

static void mnc_free ( mnc_all *all ) {
/*  data_creation_param *dc = all->dc ;*/
  mnc_vectors *vec = all->vec ; 

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

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

  free_dvector ( vec->bias , 1 , vec->N ) ; 
  free_dvector ( vec->n , 1 , vec->N ) ; 

}

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 -->
Last modified: Sat Aug 23 16:46:51 1997
<!-- hhmts end -->
*/