bnd.c

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

  double **pdpf = p->dpf ;
  double *dpr ; 
  double **pdpr = p->dpr ;
  int *lupto  = p->a->l_up_to ;         /* keeps track of what l we are up to */
  int *nlist ; 
  alist_matrix *a = p->a ;
  unsigned char *z = p->z ; 
  double **pc0 = p->pc0 ;
  double **pc1 = p->pc1 ;
  int atanhwarn=0 ; /* whether we have spat a warning this time */
                                        /* * * * * * FORWARD PASS  * * * * * */
  for ( m = 1 ; m <= p->M ; m++ )   
    lupto[m] = 0 ;                      /* check all columns are at start    */

  for ( n = 1 ; n <= p->N ; n++ ) {     /* Run through bits x[n] left->right */
    umax = a->num_nlist[n] ;            /* number of bits in this column     */
    nlist = a->nlist[n] ;               /* list of bits in this column       */
    dqc = p->dqc[n] ;                   /* contributions to be made          */

    for ( u = 1 ; u <= umax ; u ++ ) {  /* Run down column              */
      m = nlist[u] ;                    /* what row are we on ?         */
      dpf = pdpf[m] ;                   /* -- for efficiency --         */

      l = lupto[m] ;                    /* what value of l = left-right */
      lupto[m] ++ ;                     /*   are we at on this row ?    */

      dpf[l+1] = dqc[u] * dpf[l] ;      /* the forward step             */
    }
  }
/*
  for ( m = 1 ; m <= p->M ; m++ )   {
    if ( lupto[m] != a->num_mlist[m] ) {
      fprintf ( stderr , "alist visitation problem %d %d %d\n" ,
	       m , lupto[m] , a->num_mlist[m] ) ;
      pause_for_return () ; 
    }
  }
*/
                                        /* * * * * * BACKWARD PASS * * * * * */

  for ( n = p->N ; n >= 1 ; n-- ) {     /* Run back through bits n           */
    umax = a->num_nlist[n] ;            /* number of bits in this column     */
    nlist = a->nlist[n] ;               /* list of bits in this column       */
    dqc = p->dqc[n] ;                   /* contributions to be made          */
    

    for ( u = 1 ; u <= umax ; u ++ ) {  /* Run down column              */
      m = nlist[u] ;                    /* what row are we on ?         */
      dpr = pdpr[m] ;                   /* -- for efficiency --         */
      dpf = pdpf[m] ;                   

      l = lupto[m] ;                    /* what value of l = left-right */
      lupto[m] -- ;                     /*   are we at on this row ?    */

      dpr[l] = dqc[u] * dpr[l+1] ;      /* the backward step             */
      dpc = dpf[l-1] * dpr[l+1] ;
      if ( c->vpl ) {  /* vertical pass wants the logarithm */
	if ( dpc >= 1.0 ) { 
	  if ( !atanhwarn ) { fprintf ( stderr , "A" ) ; fflush ( stderr ) ; atanhwarn=1; }
	  dpc=0.5; /* hack it back to an ad hoc value */
	} else if ( dpc <= -1.0 ) {
	  if ( !atanhwarn ) { fprintf ( stderr , "A" ) ; fflush ( stderr ) ; atanhwarn=1; }
	  dpc=-0.5; /* hack it back */
	}	
	if ( !z[m] ) { pc1[n][u] = - 2.0 * atanh ( dpc )  ;
	} else { pc1[n][u] = 2.0 * atanh( dpc ) ; }
      } else { /* standard algm */
	if ( !z[m] ) {
	  pc0[n][u] = 0.5 * (1.0 + dpc) ;         /* messages to bit n */
	  pc1[n][u] = 0.5 * (1.0 - dpc) ;
	} else {
	  pc1[n][u] = 0.5 * (1.0 + dpc) ; 
	  pc0[n][u] = 0.5 * (1.0 - dpc) ; /* if space is short then only 
					     a vector is really needed here,
					     if the vertical computation
					     done immediately */
	}
      }
    }
  }
}


void vertical_pass ( bnd_param *p , bnd_control *c  ) 
{
  int u , umax ; 
  int n ;
  double  *dqc , *qt0 = p->qt0 , *qt1 = p->qt1 ;
  double         *qb0 = p->qb0 , *qb1 = p->qb1 ; 
  double **pc0 = p->pc0 ;  double **pc1 = p->pc1 ;
  double *qc0 = p->qc0 ;   double *qc1 = p->qc1 ;
  double sum , dif ; 
  alist_matrix *a = (p->a) ;
  int clipwarn=0 ; /* whether we have spat a warning this time */
  int fwarn=0 ; /* whether we have spat a warning this time */
  int swarn=0 ; /* whether we have spat a warning this time */

  for ( n = 1 ; n <= p->N ; n++ ) {     /* Run through bits x[n] left->right */
    umax = a->num_nlist[n] ;              /* number of bits in this column   */
    dqc = p->dqc[n] ;                     /* contributions to be computed    */

    qt0[0] = 1.0 - p->bias[n] ;
    qt1[0] = p->bias[n] ;
    qb0[umax+1] = 1.0 ;
    qb1[umax+1] = 1.0 ;
                                        /* * * * * * DOWNWARD PASS * * * * * */
    for ( u = 1 ; u <= umax ; u ++ ) {  /* Run down column                   */
      qt0[u] = qt0[u-1] * pc0[n][u] ;   /*     downward step                 */
      qt1[u] = qt1[u-1] * pc1[n][u] ;   /*     downward step                 */
    }
    sum = ( qt0[umax] + qt1[umax] ) ;
    if ( sum > c->tinydiv ) {
      p->q1[n] = qt1[umax] / sum ;      /* pseudoposterior prob, bit n */
    } else { 
/*      p->q1[n] = 0.49 ;  */ /* leave it as it was */
      if ( !swarn ) { fprintf ( stderr , "*" ) ; fflush ( stderr ) ; swarn=1; }
    }

                                        /* * * * * * UPWARD PASS * * * * * * */
    for ( u = umax ; u >= 1 ; u -- ) {  /* Run up column                     */
      qb0[u] = qb0[u+1] * pc0[n][u] ;   /*    upward step                    */
      qb1[u] = qb1[u+1] * pc1[n][u] ;   /*    upward step                    */

      qc0[u] = qt0[u-1] * qb0[u+1] ;    /*    everyone else's messages       */
      qc1[u] = qt1[u-1] * qb1[u+1] ;    /*    everyone else's messages       */

      sum = qc0[u] + qc1[u] ;
      dif = qc0[u] - qc1[u] ;
      if ( sum > c->tinydiv ) {
	dqc[u] = dif / sum ;              /*   difference for the next iteration */
	if ( c->dofudge ) { /* research option - does stretching or shrinking the qs help
			       the decoder? */
	  if ( fudge( &dqc[u] , c->fudge ) < 0 ) {
	    if ( !fwarn ) { fprintf ( stderr , "F" ) ;  fflush ( stderr ) ; fwarn=1; } 
	  }
	}
	if ( c->doclip ) { /* pragmatic clipping procedure to stop overflow */
	  if ( dqc[u] > c->clip ) { 
	    dqc[u] = c->clip ;
	    if ( !clipwarn ) {
	      fprintf ( stderr , "c" ) ; fflush ( stderr ) ; clipwarn=1;
	    }
	  }
	  else if ( dqc[u] < - c->clip ) { 
	    dqc[u] = - c->clip ; 
	    if ( !clipwarn ) { 
	      fprintf ( stderr , "c" ) ; fflush ( stderr ) ; clipwarn=1; 
	    }
	  }
	}
	
      } else { 
	dqc[u] = 0.0 ;
	if ( !swarn ) { fprintf ( stderr , "*" ) ; fflush ( stderr ) ; swarn=1; }
      }
    }
  }
}

/* Mon 10/6/02
   New version that uses a log scale for the vertical pass, to try to reduce roundoff errors.
   Use the identical horizontal pass (modified a bit internally),
   but take log(1+dr)/(1-dr) = log pc1[n][u] - log pc0[n][u]
   at first opportunity.

   We store the values of these quantities in    pc1[n][u]

   (actually that is done in horizontal pass now
 */
void vertical_pass2 ( bnd_param *p , bnd_control *c  ) 
{
  int u , umax ; 
  int n ;
  double  *dqc , *qt1 = p->qt1 ;
  double  *qb1 = p->qb1 ; 
  double  **pc1 = p->pc1 ;
  double  *qc1 = p->qc1 ;
  alist_matrix *a = (p->a) ;
  int clipwarn=0 ; /* whether we have spat a warning this time */
  /*  int slwarn=0 ;  whether we have spat a warning this time */

  for ( n = 1 ; n <= p->N ; n++ ) {     /* Run through bits x[n] left->right */
    umax = a->num_nlist[n] ;              /* number of bits in this column   */
    dqc = p->dqc[n] ;                     /* contributions to be computed    */

    qt1[0] =  p->lbias[n] ; 
    qb1[umax+1] = 0.0 ;

                                        /* * * * * * DOWNWARD PASS * * * * * */
    for ( u = 1 ; u <= umax ; u ++ ) {  /* Run down column                   */
      qt1[u] = qt1[u-1] + pc1[n][u] ;   /*     downward step                 */
    }
    /* make pseudoposterior from the result at bottom */
    p->q1[n] = tanh ( 0.5 * qt1[umax] )  ;      /* pseudoposterior prob, bit n */

                                        /* * * * * * UPWARD PASS * * * * * * */
    for ( u = umax ; u >= 1 ; u -- ) {  /* Run up column                     */
      qb1[u] = qb1[u+1] + pc1[n][u] ;   /*    upward step                    */

      qc1[u] = qt1[u-1] + qb1[u+1] ;    /*    everyone else's messages       */
      if ( c->dofudge ) { /* research option - does stretching or shrinking the qs help
			       the decoder? */
	dqc[u] = -tanh (  0.5 * qc1[u] * c->fudge ) ;    /*   difference for the next iteration */
      } else {
	dqc[u] = -tanh (  0.5 * qc1[u] ) ;    /*   difference for the next iteration */
      }

      if ( c->doclip ) { /* pragmatic clipping procedure to stop overflow */
	if ( dqc[u] > c->clip ) { 
	  dqc[u] = c->clip ;
	  if ( !clipwarn ) {
	    fprintf ( stderr , "c" ) ; fflush ( stderr ) ; clipwarn=1;
	  }
	}
	else if ( dqc[u] < - c->clip ) { 
	  dqc[u] = - c->clip ; 
	  if ( !clipwarn ) { 
	    fprintf ( stderr , "c" ) ; fflush ( stderr ) ; clipwarn=1; 
	  }
	}
      }
	
    }
  }
}

void bnd_tweak_prior ( bnd_param *p , bnd_control *c  ) 
{
  int n , N=p->N ; 
  int since =  c->loop - c->dotweak ;
  if ( !(since % c->tweakinc) ) { /* time for a tweak */ 
    for ( n = 1 ; n <= N ; n++ )  {
      tweak ( &(p->bias[n]) , c->tweak ) ;
    }
    fprintf ( stderr , "t" ) ; fflush ( stderr ) ;
    if ( c->tweakreset >=2 || (c->tweakreset && (since==0)))
      bnd_load_dqc ( p , c ) ;
  }
/* else nothing */
}

int tweak ( double *p , double f ) {
  double l=0.0 , q0 , q1 ;
  int status = 0 ; 
  /* turn p into a logit, scale it by f, then turn it back */
  q0 =       *p ; 
  q1 = 1.0 - *p ; 
  if ( q0 > 0.0 ) { l = - log ( q0 ) ; }
  else { status -- ; }
  if ( q1 > 0.0 ) { l += log ( q1 ) ; }
  else { status -- ; }
  if ( !status ) 
    *p = 1.0 / ( 1.0 + exp ( f * l ) ) ; 
  /* else *p gets left alone */
  return status ; 
}

int fudge ( double *d , double f ) {
  double l=0.0 , q0 , q1 ;
  int status = 0 ; 
  /* turn d into a logit, scale it by f, then turn it back */
  q0 = 1.0 + *d ; 
  q1 = 1.0 - *d ;  /* PROBLEMS occur here because often *d is equal to 1.0 */
  if ( q0 > 0.0 ) { l = - log ( q0 ) ; }
  else { status -- ; }
  if ( q1 > 0.0 ) { l += log ( q1 ) ; }
  else { status -- ; }
  if ( !status ) 
    *d = - tanh ( 0.5 * f * l ) ; 
  /* else *d gets left alone - is that wise? */
  /*  else
      fprintf ( stderr , "%8.5g", *d ) ; */
  return status ; 
}

void fudges ( double *d , double f ) { /*   NOT USED */
  double l ;
  /* turn d into a logit, scale it by f, then turn it back */
  l = d_2_logit ( *d ) ;
  *d = logit_2_d ( f * l ) ; 
}

double d_2_logit ( double d ) { /*   NOT USED */
  double l = 0.0 , q0 , q1 ;
  double BIG = 10.0 ;
  q0 = 1.0 + d ; 
  q1 = 1.0 - d ; 
  if ( q0 > 0.0 ) { l = - log ( q0 ) ; } /* if  q0=0.000000001 then -log(q0) is BIG */
  else { fprintf ( stderr , "L" ) ; fflush ( stderr ) ; l = BIG ; }
  if ( q1 > 0.0 ) { l += log ( q1 ) ; }
  else { fprintf ( stderr , "l" ) ; fflush ( stderr ) ; l += - BIG ; }
  return l ;
}

double logit_2_d ( double l ) { /*   NOT USED */
  double d ; 
  d = - tanh ( l * 0.5 ) ;
  return d ;
}

/*
<!-- hhmts start -->
Last modified: Wed Jun 11 12:14:15 1997
<!-- hhmts end -->
*/