rm.c

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

  /* Write the edges and labels of the trellis to files for gnuplotting */

  for ( l = 0 ; l <= N-1 ; l ++ ) { /* don't do state N cos you can't go from there */
    ii = pint[l+1] ;
    i = 0 ;
    do { /* looping over all states in trellis at l */

      thisv =  v[l][i] ;
      /* look at low bits to find nodes to right */
      thisv = thisv % 64 ;
      zerook = thisv % 8 ;
      oneok = thisv / 8 ;
      if(verbose >= 2)printf("%d%d ",zerook,oneok);
      if(verbose) { printf ( "%2o %4o  " , i , v[l][i] ) ; }

#define WRITEIT(a,b) \
 fprintf ( fpe , "%3d %1d %1d %4d %4d %4d %4d\n",e,l, a ,l,i,l+1,b ) ; \
 fprintf ( fpe , "%3d %1d %1d %4d %4d %4d %4d\n",e,l, a ,l+1,b,l,i ) ; \
 fprintf ( fpe , "\n"); \
 mx = (double)(l)+0.5 ; my = (double)(i + b)*0.5 ; \
  fprintf ((a?fpo:fpz), "%3d %1d %1d %8f %8f\n", e , l , a, mx, my ) 

/* a is the label 0/1 and b is the destination i or i^ii */

      if (zerook) { /* I can go to l+1, i */
	e ++ ; 
	WRITEIT(0,i);
      } 
      if ( oneok ) { /* I can go to l+1 , i^ii */
	e ++ ; 
	WRITEIT(1,(i^ii));
      }

      i = nexti[l][i] ;
    }    while ( i ) ;
    if(verbose){
      printf ("\n");
    }
  }
  fclose (fpe);
  fclose (fpo);
  fclose (fpz);
  fprintf ( stderr , "# gnuplot\nset size 0.4,0.4 ; set nokey ; set noborder; set noxzeroaxis ; set noyzeroaxis ; set noxtics; set noytics\n") ;
  fprintf ( stderr , "plot 'edges' u 4:5 w l,'ones' u 4:5 w p 4 6,'zeros' u 4:5 w p 3 3, 'edges' u 4:5 w p 7 1\n" ) ;

}

void lt_zero_r ( linear_trellis *h ) {
  int l ;
  int N = h->N ; 
  double **r = h->r ;

  for ( l = 1 ; l <= N ; l ++ ) {
      r[l][0] = 0.0 ; 
      r[l][1] = 0.0 ; 
  }
}

int lt_syndrome (  linear_trellis *h ) {
  /*  int verbose = h->verbose ;  */
  int N = h->N ;
  int *pint = h->pint ; 
  int i , l ; 
  unsigned char *xo = h->xo ; 
  /* zips across trellis to compute the syndrome */

  i = 0 ; 
  for ( l = 1 ; l <= N ; l ++ ) { /* going all the way to N+1
					is unnecessary */
    if ( xo[l] ) i ^= pint[l] ;
  }
  return (i) ; /* if this is zero then it was a valid codeword */
}

/* runs forward-backward on trellis. Mon 23/4/01 added entropy computer. */
void lt_fb ( linear_trellis *h ) {
  int verbose = h->verbose ; 
  int oneok , zerook ;
  int thisv ; 
  double thisfd, thisbd , thisbd1 ;
  int N = h->N ;
  double **q = h->q ; 
  double **r = h->r ; 
  double **Ef = h->Ef ; 
  double **Eb = h->Eb ; 
  double **Ec = h->Ec ; 
  double **Scbit = h->Scbit ; 
  double **Sc = h->Sc ; 
  double **Sf = h->Sf ; 
  double **Sb = h->Sb ; 
  double *dSdpo = h->dSdpo ; 
  double *Srest = h->Srest ; 
  double *dSdl = h->dSdl ; 
  double **fd = h->f ; 
  double **bd = h->b ; 
  double **po = h->po ; 
  int **v = h->v ; 
  int **nexti = h->nexti ; 
  int *pint = h->pint ; 
  int i , previ , l , newi ;
  double thisS   ; /* conditional entropy */
  double pthis , pthisnorm ; /* normalized probability of being in this state */
  int gobackto ;
  double condprob0 , condprob1 , prevE , pastSc ,       dSdpoA ; 
  
  /* forward pass works by going to each *next* node and adding up 
     appropriate stuff from all its precursors */

  if(verbose >= 2) { printf ( "== doing fb\n" ) ; }
  if(verbose >= 2) {
    for ( l = 1 ; l <= N ; l ++ ) 
      printf ( "q[%d] = %9.5g,%9.5g\n" , l , q[l][0] , q[l][1] ) ; 
  }
  if(verbose >= 2) printf("l   i  prob  choices(0,1)?\n");

  /* forward pass */
  for ( l = 1 ; l <= N ; l ++ ) { /* going all the way to N+1
					is unnecessary, but we do it anyway */
    i = pint[l] ;
    newi = 0 ; do { /* looping over all states in trellis at l */

      thisv = v[l][newi] ; thisfd = 0.0 ; /* v records the valid transitions */

      thisv =  thisv / 64 ;	/* look at high bits to find precursors */
      zerook = thisv % 8 ;
      oneok = thisv / 8 ;
      if(verbose >= 2)printf("%d%d",zerook,oneok);
      if (zerook) {
	/* get the juice and copy it over */
	thisfd += fd[l-1][newi] * q[l][0] ;
	if(verbose >= 2)printf(" 0<- %9.4g * %9.4g " , fd[l-1][newi] ,
			       q[l][0] ) ;
      }
      if ( oneok ) {
	previ = newi^i ; 
	thisfd += fd[l-1][previ] * q[l][1] ; 
	if(verbose >= 2)printf(" 1<- %9.4g * %9.4g " , fd[l-1][previ] ,
			       q[l][1] ) ;
      }
      fd[l][newi] = thisfd ;
      if(verbose >= 2)printf ("%2o %9.4g\n",newi,thisfd);

      /* end loop */
      newi = nexti[l][newi] ; 
    }    while ( newi ) ;


    if(verbose >= 2){
      printf ("\n");
    }
  }
  if ( h->computeS ) {
    /* compute the entropy of the posterior distribution over paths */
    /* Mon 23/4/01  correct only for Gallager type trellis */
    h->S = 0.0 ; 
    /* Make a note of the marginal likelihood, i.e. the probability of reaching
	 the (presumed unique) final state of trellis */
    h->ML = thisfd ;
    zero_energy_trellis(h) ; 
  }

  if(verbose >= 2)printf("backpass\n");
  if(verbose >= 2)printf("l   i  prob  choices\n");
  gobackto = (h->computeS) ? 0 : 1 ;
  for ( l = N-1 ; l >= gobackto ; l -- ) {
    i = pint[l+1] ;
    newi = 0 ; do { /* looping over all states in trellis at l */

      thisv =  v[l][newi] ;
      thisfd = 0.0 ; 
      /* look at low bits to find nodes to right */
      thisv = thisv % 64 ;
      zerook = thisv % 8 ;
      oneok = thisv / 8 ;
      if(verbose >= 2)printf("%d%d ",zerook,oneok);
      if (zerook) {
	thisfd += bd[l+1][newi] * q[l+1][0] ;
      } 
      if ( oneok ) {
	thisfd += bd[l+1][newi^i] * q[l+1][1] ; 
      }
      bd[l][newi] = thisfd ;
      if(verbose >= 2)printf ("%2o %9.4g\n",newi,thisfd);

      if( h->computeS ) {
	/* ALSO compute mean energy, conditioned on being in state  */
	if (zerook) { 
	  Eb[l][newi] += ( Eb[l+1][newi]-log(q[l+1][0])) * q[l+1][0] * bd[l+1][newi] / bd[l][newi] ;
	}
	if ( oneok ) {
	  Eb[l][newi] += ( Eb[l+1][newi^i] -log(q[l+1][1])) * q[l+1][1] * bd[l+1][newi^i] / bd[l][newi] ;
	}
	Sb[l][newi] = Eb[l][newi] + log( bd[l][newi] ) ;
	if(verbose>=2) {
	  printf(" mean back-energies conditional on state%d=%d %9.5g\n", l , newi, Eb[l][newi] );
	  printf(" entropy of future condl on state%d=%d: %9.5g =%9.5g %9.5g\n", l , newi, Sb[l][newi] , Eb[l][newi] , log( bd[l][newi] ) );
	}
      } 

       /* end loop */
      newi = nexti[l][newi] ;
    }    while ( newi ) ;
    if(verbose >= 2){
      printf ("\n");
    }
  }
  lt_zero_r(h) ; 
  if(verbose >= 2)printf("readout r\n");
  if(verbose >= 2)printf("col state\n");
  if(verbose >= 2)printf("l   i  prob  choices\n");
  for ( l = 1 ; l <= N ; l ++ ) {
    i = pint[l] ;
    previ = 0 ; do { /* looping over all states in trellis at l-1 */

      thisv = v[l-1][previ] ; 
      thisfd = fd[l-1][previ]; /* get the juice */
      if(verbose >= 2) { printf ( "%2d %2o %9.4g  " , l-1 , previ , thisfd ) ; }

      thisv =  thisv % 64 ;	/* look at low bits */
      zerook = thisv % 8 ;
      oneok = thisv / 8 ;
      if(verbose >= 2)printf("%d%d",zerook,oneok);
      if (zerook) { /* whether there is a zero transition from the *previous* state */
	newi = previ ; /* consider transition zero */
	thisbd = bd[l][newi] ;
	r[l][0] += thisfd * thisbd ;
	if(verbose >= 2)printf ( " meet %9.4g ",thisbd);
      }
      if ( oneok ) {
	newi = previ^i ; /* where in the current state a 1 takes us to */
	thisbd1 = bd[l][newi] ;
	    
	r[l][1] += thisfd * thisbd1 ;
	if(verbose >= 2)printf ( " meet %9.4g ",thisbd1);
      }
      if(verbose >= 2)printf ("\n");
      
      if ( h->computeS ) { 
	/* compute conditional entropy of 2nd state given first state */
	/* find probability of being in first state */
	if( ( zerook && oneok ) ) {
	  pthis = thisfd * bd[l-1][previ];
	  pthisnorm = pthis / h->ML ;
	  condprob0 =  thisbd  * q[l][0] ;
	  condprob1 =  thisbd1 * q[l][1] ;
	  thisS = H2( condprob0 / ( condprob0 + condprob1 ) );
	  if( l==1 ) { h->S1 = thisS ; } /* special debugging */
	  if (verbose >= 2) 
	    printf ("P(state(%d)=%d) = %9.4g /%3.2g = %9.4g ; H(S%d|s=%d)=%9.3g (%6.3g)\n",
		    l-1, previ , pthis , h->ML , pthisnorm , l , previ , thisS  , thisS/log(2.0) );
	  h->S += thisS * pthisnorm ; 
	}
	/* ALSO compute mean energy, conditioned on being in state  */
	prevE = Ef[l-1][previ] ; 
	if (zerook) { /* whether there is a zero transition from the *previous* state */
	  Ef[l][previ] += (prevE-log(q[l][0])) * q[l][0]  *  thisfd  / fd[l][previ] ;
	}
	if ( oneok ) {
	  Ef[l][newi] +=  (prevE-log(q[l][1])) * q[l][1]  *  thisfd / fd[l][newi] ;
	}
      } 
      /* end loop */
      previ = nexti[l-1][previ] ; 
    }    while ( previ ) ;
    

    if ( h->computeS ) {
      /* compute the entropy of the posterior distribution over paths */
      if(verbose >= 2){
	printf("%d : entropy so far %9.5g (%9.5g bits)\n", l , h->S , h->S/log(2.0) );
	printf(" mean f-energies conditional on state%d=0/1 %9.5g %9.5g\n", l , Ef[l][0] , Ef[l][1] );
      }
      previ = 0 ; do { /* looping over all  current states in trellis  */
	pastSc =  Ef[l][previ] + log( fd[l][previ] );
	Sc[l][previ] = Sb[l][previ] + pastSc ;
	Sf[l][previ] = pastSc ;
	if(verbose >= 2){
	  printf(" entropy of past condl on state%d=%d: %9.5g =%9.5g %9.5g\n", l , previ , pastSc , Ef[l][previ] , log( fd[l][previ] ) );
	}
	previ = nexti[l][previ] ; 
      }    while ( previ ) ;
	/* and find the conditional entropy given this bit l is
	   either 0 or 1 */
      previ = 0 ; do { /* looping over all  current states in trellis  */
	thisv = v[l-1][previ] ; /* repeat procedures used to find r1 and r0 */
	thisfd = fd[l-1][previ]; /* get the juice */
	thisv =  thisv % 64 ;	/* look at low bits */
	zerook = thisv % 8 ;
	oneok = thisv / 8 ;
	if (zerook) { /* whether there is a zero transition from the *previous* state */
	  newi = previ ; /* consider transition zero */
	  thisbd = bd[l][newi] ;
	  Ec[l][0] += thisfd * thisbd *  ( Ef[l-1][previ] + Eb[l][newi] ) / r[l][0] ;
	}
	if ( oneok ) {
	  newi = previ^i ; /* where in the current state a 1 takes us to */
	  thisbd1 = bd[l][newi] ;
	  Ec[l][1] += thisfd * thisbd1 * ( Ef[l-1][previ] + Eb[l][newi] ) / r[l][1] ;
	}
	/* end loop */
	previ = nexti[l-1][previ] ; 
      }    while ( previ ) ;
      for ( previ = 0 ; previ <= 1 ; previ ++ ) { /* run over bit values now */
	Scbit[l][previ] = Ec[l][previ] + log( r[l][previ] );
	if(verbose >= 2){
	  printf(" entropy of all condl on bit%d=%d: %9.5g =%9.5g %9.5g\n", l , previ , Scbit[l][previ] , Ec[l][previ] , log( r[l][previ] ) );
	}
	/* and finally, compute the derivative wrt the parameter q of
	   the JOINT entropy, using these conditional entropies */
	po[l][previ] =  r[l][previ] * q[l][previ]/( r[l][0]* q[l][0] + r[l][1]* q[l][1] ) ;
      } /* derivative wrt posterior is: */
      /* DERIV */
      /* DERIV */
      /* DERIV */
      Srest[l] = po[l][1] *  Scbit[l][1] + po[l][0] * Scbit[l][0] ;
      dSdpoA =  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 ;
      dSdl[l] = dSdpoA ;
      dSdpo[l] = dSdpoA ;
      if(verbose >= 2){ /* see also RMtest.c: show */
	printf(" DERIV (w r t posterior) dS/dp%d = %9.5g (%9.5g %9.5g);\n\t\t\t  w.r.t. logit %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  ) ;
      }
      
    }
  }

  if ( verbose ) {
    printf ("== outcome (r) [fn[2]=%9.6g]\n", q[2][1] );
    printf ("    ri0       ri1       norm0    norm1   |  post0  post1\n--------------------------------------\n" );
    for ( i = 1 ; i <= N ; i++ ) {
      printf(" %9.4g %9.4g %9.4g %9.4g" , r[i][0] , r[i][1] , r[i][0]/( r[i][0] + r[i][1] ) , r[i][1]/( r[i][0] + r[i][1] )) ;
      printf(" | %9.4g %9.4g\n" , q[i][0] * r[i][0]/( q[i][0] * r[i][0] + q[i][1] * r[i][1] ) , q[i][1] * r[i][1]/( q[i][0] * r[i][0] + q[i][1] * r[i][1] )) ;
    }
  }
  if(verbose >= 3) { pause_for_return () ; }

  /* since r's are about to multiplied by priors and things there is no point in 
     our normalizing them. We leave the r's where they are and return. */

}

/* Adjusts the incoming q's (via lq) #
   such that the required posteriors (as defined by pob) are matched */
void lt_match ( linear_trellis *h ) {
  int verbose = h->verbose ; 
  double **q = h->q ; 
  double **po = h->po ;
  double *lpo = h->lpo ;
  double *lpob = h->lpob ;
  double *lq = h->lq ;
  double *fudge = h->fudge ;
  int fudgeperiod = h->fudgeperiod ;
  int N = h->N ; 
  int i ,  t , T = h->T ;

  /* at start, user should have supplied
     lpob (and pob maybe) (target values)

     Also an initial condition for lq, could be provided, but to start with
     we will set it from lpob
  */
  /* set initial condn for lq */
  for ( i=1;i<=N;i++ ) {
    lq[i] = lpob[i] ;
  }
  for ( t=1; t<=T ; t++ ) {
    /* transfer lq to q */
    for ( i=1;i<=N;i++ ) {
      q[i][1] = 1.0/(1.0+exp(-lq[i])) ;
      q[i][0] = 1.0/(1.0+exp(lq[i])) ;
    }
    /* run fb */
    lt_fb ( h ) ;
    /* extract lpo */
    for ( i=1;i<=N;i++ ) {
      lpo[i] = log( po[i][1] / po[i][0] ) ; 
    }
    /* report */
    if(verbose>2) {
      printf("Iteration %d\ni      lq       lpob       lpo       delta\n====================================\n",t) ;
      for ( i=1;i<=N;i++ ) {
	printf("%d %9.6g %9.6g %9.6g %9.6g\n" , i, lq[i] , lpob[i] , lpo[i], lpob[i] - lpo[i]  ) ;
      }
    }    else if(verbose>-1) {
      if (t==1) { /* headings */
	fprintf(stderr,"# gnuplot \n plot 'dat' u 1:5 w linesp,  'dat' u 1:9 w linesp\n" ) ;
      }
      if (!((t+9)%10)) { /* headings */
	printf("# Its ");
	for ( i=1;i<=N;i++ ) 
	  printf("   lq[%d]   lpob[%d]    lpo[%d]    delta[%d]",i,i,i,i) ;
	printf( "\n" );
      }
      printf( "%d " ,t);
      for ( i=1;i<=N;i++ ) {
	printf("%9.5g %9.5g %9.5g %9.5g " , lq[i] , lpob[i] , lpo[i], lpob[i] - lpo[i]  ) ;
      }
      printf( "\n" );
    }
    /* update lq */
    for ( i=1;i<=N;i++ ) {
      lq[i] += fudge[t%fudgeperiod] * (lpob[i] - lpo[i])  ;
    }
  }
}
  

double H2 ( double x0 ) {
  double tmpd , l0 , l1 , x1 ;
  x1 = 1.0 - x0 ; 
  if ( (x1 > 0.0) && (x0 > 0.0 )  ) {
    l1 = - log ( x1 ) ; 
    l0 = - log ( x0 ) ;

    tmpd =  ( x0 * l0 + x1 * l1 ) ; 
    return tmpd ;
  }
  else return 0.0 ; 
}