rmtest.c

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

/*
  RMtest.c  April 1997 DJCM (c)

  demonstrates a few uses of RM.c, which 
  contains routines for creating RM trellises and Hamming too,
  and for doing belief propagation on them

  the trellis goes from l=0 to N, with 
  the number of emitted bits being N. 
  from 0 to 1, bit 1 is emitted.
  from N-1 to N, bit N is emitted. 
*/

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

void   main ( int , char ** ) ;
void   show ( linear_trellis *t , int i ) ;

/*
        MAIN
                     */
void main ( int argc, char *argv[] )
{
  double l0 = -2.19722 , l ;
  linear_trellis t ;
  int verbosity = 0 ; 
  int type = 1 ; /* type 1 makes hamming and RM codes; type 0 makes Gallager */
  double fn = 0.1 ;
  double delta = 0.001 ;
  double origS , origS1 , origSrest ; 
  int i ; 
  int flipone = 1 ;
  
  /* these control the convergence- of- the- lagrange- multiplier experiments*/
  int N = 7 ; /*  trellis length */
  int diversity = 1 ;
  double fudgedef = 0.68 ;
  double slope = -0.2 ;
  l0 = 3;
  /* *************** */

  if(0) {
    printf ( "first, let's show a load of trellises\n");
    printf ( "============================\n");
    make_linear_trellis ( &t , 7 , 1, verbosity ) ;
  }
  make_linear_trellis ( &t , N , 0, verbosity ) ;

  t.computeS = 1 ;
  make_energy_trellis ( &t );

  if(verbosity) {
    pause_for_return();
    printf ( "Now see fb happening, when the received vector looks like 0,0,0,0,0\n");
    /* invent a q vector */
  }
  l = l0 ;
  for ( i= 1; i<= t.N ; i++ ) {
    fn = 1.0/(1.0 + exp( - l )) ;
    if(i==1) {
      t.q[i][1] = 1.0 - fn ; t.q[i][0] = fn ;
    } else {      t.q[i][0] = 1.0 - fn ; t.q[i][1] = fn ;
    } 
    l += slope ;
    t.lpob[i] = l ;
  }
  
  t.T = 40 ;
  if ( diversity ) {
    t.fudgeperiod = 3 ;
    t.fudge[3] = 0.1  ;
    t.fudge[4] = 0.25  ;
    t.fudge[0] = 0.5  ;
    t.fudge[1] = 1.0  ;
    t.fudge[2] = 0.68 ;
  } else {
    t.fudgeperiod = 1 ;
    t.fudge[0] = fudgedef  ;
    t.fudge[1] = fudgedef  ;
  }
  lt_match(&t) ;
  
  exit(0);
  lt_fb ( &t ) ;
  show (&t,2 ) ;
  origS = t.S ;
  origS1 = t.S1 ;
  origSrest = origS - origS1  ;
  if(verbosity) {
    pause_for_return();
  }
  l=l0; 
  for ( i=1; i<=t.N ; i++ ) {
    if (verbosity>-1) { printf ("PERTURB %d\n" , i ) ; }
    l += delta ;
    fn = 1.0/(1.0 + exp( - l )) ; 
    t.q[i][0] = 1.0 - fn ; t.q[i][1] = fn ;

    lt_fb ( &t ) ; 
    show (&t,1 ) ;
    if(verbosity) {
      pause_for_return();
    }
    printf(" EMPIRICAL DERIVS %d %9.7g ( %9.7g %9.7g)\n" , i ,  (t.S - origS) / delta ,  (t.S1 - origS1) / delta ,  ((t.S-t.S1) - origSrest) / delta ); 
    l -= delta ;
    fn = 1.0/(1.0 + exp( - l )) ; 
    t.q[i][0] = 1.0 - fn ; t.q[i][1] = fn ;
    l +=0.05 ;
  }
    
  if ( flipone ) { 
    printf ( "Now the received vector is (1,0,0,0,0..) \n");
    fn =  t.q[1][1] ;  t.q[1][1] = t.q[1][0] ; t.q[1][0]  = fn ;
    lt_fb ( &t ) ; 
    pause_for_return();
  }
  printf ( "Now see fb happening at a sequence of different noise levels\n");
  pause_for_return();
  for ( fn = 0.02 ; fn <= 0.5 ; fn += 0.04 ) {
    printf ("fn = %9.5g ========== \n",fn) ; 
    for ( i= 1; i<= t.N ; i++ ) {
      t.q[i][0] = 1.0 - fn ; t.q[i][1] = fn ;
    }
    if ( flipone ) { t.q[1][0] = fn ; t.q[1][1] =  1.0 - fn ; }
    lt_fb ( &t ) ; 
    pause_for_return();
  }


  free_energy_trellis ( &t );

}

void   show ( linear_trellis *h , int localverbose ) {
  int verbose = h->verbose ; 
  int N = h->N ;
  double **Scbit = h->Scbit ; 
  double *dSdpo = h->dSdpo ; 
  double *dSdl = h->dSdl ; 
  double *Srest = h->Srest ; 
  double **po = h->po ; 
  double **r = h->r ; 
  double **q = h->q ; 
  int l  ;

  for ( l = 1 ; l <= N ; l ++ ) {
    printf("%9.5g, ", q[l][1] ) ;
  }
  printf(" -> Entropy %9.5g (%9.5g bits)\n", h->S , h->S/log(2.0) );
  if(localverbose) {
    printf(" Entropy of whole trellis : %9.5g\n" , h->S ) ;
    printf(" Entropy of trellis step 1: %9.5g\n" , h->S1 ) ;
    printf(" Entropy of rest          : %9.5g\n" , h->S - h->S1 ) ;
  }

  if(localverbose>1) {
    printf("\n DERIV (w r t posterior) \n dS/dp(i) \t\t  w.r.t. logit\n" );
    for ( l = 1 ; l <= N ; l ++ ) {

       printf("%d %9.5g (%9.5g %9.5g)\t %9.5g (%9.5g %9.5g)\n",
	       l , dSdpo[l] , log ( po[l][0] / po[l][1] ) ,
	       (r[l][1] * (Scbit[l][1]-Srest[l]) - r[l][0] * ( Scbit[l][0] - Srest[l] ))/h->ML,
	       dSdl[l] , log ( po[l][0] / po[l][1] )  * po[l][0] * po[l][1] ,
	       ( r[l][1] * (Scbit[l][1]-Srest[l]) - r[l][0] * ( Scbit[l][0] - Srest[l] )  )  * q[l][0]*q[l][1]/h->ML  ) ;

      
    }
  }

}

/*
<!-- hhmts start -->
Last modified: Tue Jul 15 17:14:44 1997
<!-- hhmts end -->
*/