mnc2.c

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

  dc->true_fs = (double) count_s / (double) dc->N ; 
  if( dc->DEMO || dc->verbose >= 1 ) {
    printf ( "generating string\n" ) ;
    printoutcvector ( v->s , 1 , dc->N ) ; 
  }

  /* Use A to generate t = A s */
  if ( dc->MNC==2 || dc->MNC==4 ) {
    mult_cm_cv ( code->C2 , v->s , v->u , 1 , dc->N ) ;
    mult_cm_cv ( code->C1 , v->u , v->t , 1 , dc->N ) ;
  } else if ( dc->MNC == 1 || dc->MNC == 3 ) {
    mult_cm_cv ( code->C1 , v->s , v->t , 1 , dc->N ) ;
  } else {
    alist_times_cvector_mod2 ( &(p->a) , v->s , v->t ) ; 
  }
  for ( m = 1 ; m <= dc->M ; m++ ) {
    count_t += v->t[m] ; 
  }
  dc->fte = (double) count_t / (double) dc->M ; 

  for ( count_err = 0 , m = 1 ; m <= dc->M ; m++ ) {
    count_err += ( v->noise[m] = ( ranu() < dc->fn ) ? 1 : 0 ) ; 
    v->r[m] = v->noise[m] ^ v->t[m] ;
  }
  dc->true_fn = (double) count_err / (double) dc->M ; 

  if ( dc->MNC==4 ) {
    mult_cm_cv ( code->C1I , v->r , v->cr , 1 , dc->N ) ;
  }

  if( dc->DEMO || dc->verbose >= 1 ) {
    printf ( "generating noise -> received:\n" ) ;
    printoutcvector ( v->noise , 1 , dc->N ) ; 
    printoutcvector ( v->r , 1 , dc->N ) ; 
  }

  return status ; 
}

  
static int evaluate_feasibility ( data_creation_param *dc , fe_min_param *p )    
{
  double fs =   dc->fs ;
  double fn =   dc->fn ;
  int status = 0 ;
  int MNC = dc->MNC ; 
  int tmpi ; 
  double ft , tmpd ; 

/*
   compute Shannon capacity and compare with info content  
*/

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

/* theoretical calculation of ft, using average number of 1's ; 
  only appropriate to C1 and C3 */
  if ( MNC == 1 || MNC == 3 ) {
    tmpi = p->a.tot / p->a.M ;
    for ( ft = 0.0 ; tmpi > 0 ; tmpi -- )
      ft = fs * ( 1- ft ) + (1 - fs ) * ft ; 

    tmpd = fn * (1.0 -ft) + (1-fn) * ft ; 
    printf ( "fs = %f , ft = %f (theo), ft = %f (emp), fn = %f , fz = %f (theo)\n" , fs, ft , dc->fte , fn, tmpd) ; 
    printf ( "H(s) = %f, H(z;t) = %f = %f - %f ; capacity of channel 1 - H(n) = %f. " , 
	    h2(fs) , h2(tmpd) - h2(fn) , h2(tmpd) , h2(fn) , 1 - h2(fn) ) ; 
    if ( (double ) dc->N * h2(fs) > (double ) dc->M * ( h2(ft) - h2(fn) ) ) {
      printf ( "task impossible because output ensemble is suboptimal\n" ) ; 
      status -= 10 ; 
    }
    else printf ( "\n" ) ; 
  }

  return status ; 
}


static void dc_defaults ( data_creation_param *dc ) 
{
  dc->MESSAGE = 1 ; 
  dc->pbm_o = 0 ; 
  dc->pbm_xv = 0 ; 
  dc->mseed = 1234 ; 
  dc->vseed = 9876 ; 
  dc->N = 20 ; 
  dc->M = 0 ; 
  dc->MNC = 0 ; 
  dc->fs = 0.1 ;     /* fraction of signal bits AND/OR assumed prior */
  dc->fpol = 0.1 ; 
  dc->npol = 0 ; 
  dc->fn = 0.15 ;    /* fraction of noise bits AND/OR assumed likelihood */
  dc->nn = 0 ;                        /* if !0, number of noise bits */
  dc->ns = 0 ;                        /* if !0, number of signal bits */
  strcpy ( dc->pbm_ofile , "A.pbm" ) ; 
  dc->DEMO = 0 ; 
  dc->verbose = 0 ;
  dc->pC1.per_row = 5 ; 
  dc->pC2.per_row = 5 ; 
}

static int process_command ( int argc , char **argv , data_creation_param *dc , fe_min_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 --
#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], "-DV") == 0 )        {
      dc->verbose = 1;
    }
    else  if ( strcmp (argv[i], "-FV") == 0 )        {
      c->verbose = 1;
    }
    else  if ( strcmp (argv[i], "-FVV") == 0 )        {
      c->verbose = 2 ;
    }
    else  if ( strcmp (argv[i], "-VV") == 0 )        {
      c->verbose = 2;
      dc->verbose = 2;
    }
    else  if ( strcmp (argv[i], "-SV") == 0 )        {
      c->seqverbose = 1;
    }
    else  if ( strcmp (argv[i], "-SVV") == 0 )        {
      c->seqverbose = 2;
    }
    else  if ( strcmp (argv[i], "-CG") == 0 )        {
      c->CG = 1;
    }
    else  if ( strcmp (argv[i], "-NL") == 0 )        {
      c->NL = 0; 
    }
    else  if ( strcmp (argv[i], "-DEMO") == 0 )        {
      c->DEMO = 1; 
    }
    else  if ( strcmp (argv[i], "-nl") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(c->NL)); 
    }
    else  if ( strcmp (argv[i], "-mess") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(dc->MESSAGE)); 
    }
    else  if ( strcmp (argv[i], "-n") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(dc->N)); 
    }
    else  if ( strcmp (argv[i], "-mnc") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(dc->MNC) ); 
    }
    else  if ( strcmp (argv[i], "-opt") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(c->opt)); 
    }
    else  if ( strcmp (argv[i], "-norder") == 0 ) {
      if ( i + 1 == argc ) { ERROR1; }
      else cs *= sscanf(argv[++i], "%d", &(c->NORDER) ) ; 
    }
    else  if ( strcmp (argv[i], "-itmax") == 0 ) {
      if ( i + 1 == argc ) { ERROR1; }
      else cs *= sscanf(argv[++i], "%d", &(c->itmax) ) ; 
    }
    else   if ( strcmp (argv[i], "-seed") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%ld", &(c->seed)); 
    }  
    else   if ( strcmp (argv[i], "-mseed") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%ld", &(dc->mseed)); 
    }  
    else   if ( strcmp (argv[i], "-vseed") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%ld", &(dc->vseed)); 
    }  
    else  if ( strcmp (argv[i], "-seqp") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(c->seq_period)); 
    }
    else if ( strcmp (argv[i], "-e") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &(c->epsilon)); 
    }
    else if ( strcmp (argv[i], "-fn") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &(dc->fn) ); 
    }
    else if ( strcmp (argv[i], "-sdx") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &(c->sdx)); 
    }
    else if ( strcmp (argv[i], "-betap") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &(c->betap)); 
    }
    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 if ( strcmp (argv[i], "-ftol") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &(c->ftol)); 
    }
    else if ( strcmp (argv[i], "-fs") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &(dc->fs)); 
    }
    else if ( strcmp (argv[i], "-fpol") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &(dc->fpol)); 
    }
    else if ( strcmp (argv[i], "-t1") == 0 ) {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(dc->pC1.per_row ) ); 
    }
    else if ( strcmp (argv[i], "-t2") == 0 ) {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(dc->pC2.per_row ) ); 
    }
    else if ( strcmp (argv[i], "-pheading") == 0 ) {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(c->pheading_period)); 
    }
    else  if ( strcmp (argv[i], "-o") == 0 )	{
      if ( i + 1 == argc ) 	    { ERROR1;      }
      else {
	strcpy(c->outfile, argv[++i]);
	if ( !(c->printout) ) c->printout = 1 ;
      }
    }
    else  if ( strcmp (argv[i], "-xv") == 0 )	{
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(dc->pbm_xv )); 
    }
    else  if ( strcmp (argv[i], "-pbm") == 0 )	{
      if ( i + 1 == argc ) 	    { ERROR1;      }
      else {
	strcpy(dc->pbm_ofile, argv[++i]);
	if ( !dc->pbm_o ) dc->pbm_o = 1 ;
      }
    }
    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 , dc , c ) ;
  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 ,
			 data_creation_param *dc   ,
			 fe_min_control      *c    
)
{
  fprintf( fp, "Usage: %s ",argv[0]);
  fprintf( fp, "datafile [optional arguments]");

  NLNE; fprintf( fp, " Verbosity of run-time output: <defaults>");
  DNT;  fprintf( fp, "-V | -VV                 (verbose or very verbose)"); 
  DNT;  fprintf( fp, "-FV | -FVV               (verbose in fe_min)"); 
  DNT;  fprintf( fp, "-DV | -DVV               (verbose in data creation)"); 
  DNT;  fprintf( fp, "-SV | -SVV               (verbose in seq_opt)"); 
  DNT;  fprintf( fp, "-DEMO                    (only do demo of MNC) "); 
  DNT;  fprintf( fp, "-CG                      (check gradient) "); 
  DNT;  fprintf( fp, "[-e epsilon]     <%g> (step size for gradient)", c->epsilon); 
  DNT;  fprintf( fp, "-pheading period <%d>" , c->pheading_period ); 

  NLNE; fprintf( fp, " Data creation:" ) ; 
  DNT;  fprintf( fp, "-mnc MNC         <%d> (which code to use [0;1,2,3,4]) " , dc->MNC ); 
  DNT;  fprintf( fp, "-n N             <%d> (number of bits in string) " , dc->N ); 
  DNT;  fprintf( fp, "-m M             <%d> (number of bits in transmission - only applicable if MNC=0) " , dc->M ); 
  DNT;  fprintf( fp, "-t1 t1           <%d> (C1 number of 1s per row) " , dc->pC1.per_row ); 
  DNT;  fprintf( fp, "-t2 t2           <%d> (C2 [if appropriate] number of 1s per row) " , dc->pC2.per_row ); 
  DNT;  fprintf( fp, "-fn fn           <%g> (error probability) " , dc->fn ); 
  DNT;  fprintf( fp, "-fs fs           <%g>  (fraction in string)", dc->fs); 
  DNT;  fprintf( fp, "-pbm pbm_ofile   <%s>  (write A.pbm)", dc->pbm_ofile);  
 DNT;  fprintf( fp, "-xv 0/1          <%d>  (xv it too)", dc->pbm_xv ); 

  NLNE; fprintf( fp, " Inference:" ) ; 
  DNT;  fprintf( fp, "-sdx sdx         <%3g> (sd of initial state)" , c->sdx ); 
  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" ); 
  DNT;  fprintf( fp, "-betap betap     <%3g> (prior's beta)" , c->betap ); 
  DNT;  fprintf( fp, "-opt optimizer   <%d>   (0=frp,1=macopt,2=jumpopt,3=seqopt,4=betaopt)", c->opt); 
  DNT;  fprintf( fp, "-itmax itmax     <%d>   (num of linmins or jumps per optimization)", c->itmax); 
  DNT;  fprintf( fp, "-norder n        <%d>   (for seqopt - 0/1/2 = 1..N / random / suspicion / suspicion/num )", c->NORDER); 
  DNT;  fprintf( fp, "-mess MESSAGES   <%d>   (number of messages to send)", dc->MESSAGE); 
  DNT;  fprintf( fp, "-nl loops        <%d>   (number of cg runs)", c->NL); 
  DNT;  fprintf( fp, "-ftol ftol       <%g> ", c->ftol ); 
  DNT;  fprintf( fp, "-seed seed "); 
  DNT;  fprintf( fp, "-mseed mseed "); 
  DNT;  fprintf( fp, "-vseed vseed "); 
  DNT;  fprintf( fp, "-seqp seq_period   <%d>  (how many loops per check of delta F)" , c->seq_period ); 
  NLNE; fprintf( fp, " Files:");
  DNT;  fprintf( fp, "-o outfile               "); 
  fprintf( fp, "\n");
  return ;
}
#undef DNT
#undef NLNE


int cmp_thing ( const void *a , const void *b ) {
  const thing *ma = (const thing*)a ;
  const thing *mb = (const thing*)b ;
  return ( ( ma->val - mb->val ) ) ; 
}

double h2 ( double x ) {
  double tmp ; 
  tmp = x * log ( x ) + ( 1.0 - x ) * log ( 1.0 - x ) ;
  return - tmp / log ( 2.0) ; 
}