rm.c

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

/*
  rm.c derived from  gh.c
  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 zero_linear_trellis ( linear_trellis *h ) { /* this is a crude one,
						  zeroing every entry */
  double **f = h->f ; 
  double **b = h->b ; 
  double **r = h->r ; 
  int N = h->N ; 
  int I = h->I ; 
  int i ;
  int verbose = h->verbose ; 

  for ( i = 1 ; i <= N ; i++ ) {
    set_dvector_const( r[i] , 0 , 1 , 0.0 ) ;
  }
  for ( i = 0 ; i <= N ; i++ ) {
    if ( verbose >= 12 ) printf ( "%d %9.6g %9.6g\n"  , i , f[i][0]  , f[i][1] ) ; 
    set_dvector_const( f[i] , 0 , I , 0.0 ) ;
    set_dvector_const( b[i] , 0 , I , 0.0 ) ;
  }
  f[0][0] = 1.0 ;
  b[N][0] = 1.0 ;
}

void free_linear_trellis ( linear_trellis *h ) {
  int N = h->N ; 
  int M = h->M ; 
  int I = h->I ; 

  free_ivector ( h->pint , 1 , N ) ;    
  free_cmatrix ( h->p , 1 , M , 1 , N ) ;  
  free_imatrix ( h->v , 0 , N+1 , 0, I ) ;      
  free_imatrix ( h->nexti , 0 , N+1 , 0, I ) ;      
  free_dmatrix ( h->q , 1 , N , 0, 1 ) ;    
  free_dmatrix ( h->r , 1 , N , 0, 1 ) ;    
  free_dmatrix ( h->f , 0 , N+1 , 0, I ) ;    
  free_dmatrix ( h->b , 0 , N+1 , 0, I ) ;    
  free_dmatrix ( h->po , 0 , N+1 , 0, 1 ) ;   /* Tue 24/4/01 */ 
}

void make_energy_trellis (  linear_trellis *h ) {
  int N = h->N ; 
  int I = h->I ;
  
  h->Ef = dmatrix ( 0 , N+1 , 0 , I ) ;      
  h->Eb = dmatrix ( 0 , N+1 , 0 , I ) ;      
  h->Ec = dmatrix ( 0 , N+1 , 0 , I ) ;      
  h->Sc = dmatrix ( 0 , N+1 , 0 , I ) ;      
  h->Scbit = dmatrix ( 0 , N+1 , 0 , 1 ) ;      
  h->Sf = dmatrix ( 0 , N+1 , 0 , I ) ;      
  h->Sb = dmatrix ( 0 , N+1 , 0 , I ) ;      
  h->Srest = dvector ( 0 , N ) ;      
  h->dSdpo = dvector ( 0 , N ) ;      
  h->dSdl  = dvector ( 0 , N ) ;      

  h->lq   = dvector ( 0 , N ) ;      
  h->pob  = dvector ( 0 , N ) ;      
  h->lpob = dvector ( 0 , N ) ;      
  h->lpo  = dvector ( 0 , N ) ;      
}

void zero_energy_trellis ( linear_trellis *h ) { /* this is a crude one,
						  zeroing every entry */
  double **Ef = h->Ef ; 
  double **Ec = h->Ec ; 
  double **Eb = h->Eb ; 
  double **Sc = h->Sc ; 
  double **Sf = h->Sf ; 
  double **Sb = h->Sb ; 
  int N = h->N ; 
  int I = h->I ; 
  int i ;
  int verbose = h->verbose ; 

  for ( i = 0 ; i <= N ; i++ ) {
    set_dvector_const( Ef[i] , 0 , I , 0.0 ) ;
    set_dvector_const( Ec[i] , 0 , I , 0.0 ) ;
    set_dvector_const( Eb[i] , 0 , I , 0.0 ) ;
    set_dvector_const( Sc[i] , 0 , I , 0.0 ) ;
    set_dvector_const( Sf[i] , 0 , I , 0.0 ) ;
    set_dvector_const( Sb[i] , 0 , I , 0.0 ) ;
  }
  set_dvector_const( h->Srest , 0 , N , 0.0 ) ;
  set_dvector_const( h->dSdpo , 0 , N , 0.0 ) ;
  set_dvector_const( h->dSdl  , 0 , N , 0.0 ) ;
}

void free_energy_trellis (  linear_trellis *h ) {
  int N = h->N ; 
  int I = h->I ; 

  free_dmatrix ( h->Ef , 0 , N+1 , 0, I ) ;    
  free_dmatrix ( h->Eb , 0 , N+1 , 0, I ) ;    
  free_dmatrix ( h->Ec , 0 , N+1 , 0, I ) ;    
  free_dmatrix ( h->Sc , 0 , N+1 , 0, I ) ;    
  free_dmatrix ( h->Scbit , 0 , N+1 , 0, 1 ) ;    
  free_dmatrix ( h->Sf , 0 , N+1 , 0, I ) ;    
  free_dmatrix ( h->Sb , 0 , N+1 , 0, I ) ;    
  free_dvector ( h->Srest , 0 , N ) ;    
  free_dvector ( h->dSdpo , 0 , N ) ;    
  free_dvector ( h->dSdl , 0 , N  ) ;    

  free_dvector ( h->lq    , 0 , N  ) ;    
  free_dvector ( h->pob   , 0 , N  ) ;    
  free_dvector ( h->lpob  , 0 , N  ) ;    
  free_dvector ( h->lpo   , 0 , N  ) ;    
}

void make_linear_trellis ( linear_trellis *h , int N , int type , int verbose ) {
  int **ff , **bb ; /* temporary objects used to create v */
  int M=-1 , K , I ; 
  int c , i , d , w , newi , lostbit ; /* generic integers */
  unsigned char **p ; /* will = h->p ;  */
  int *pint ; /* = h->pint ; */
  int  ** v ; 
  int **nexti ;
  double **q, **r ; 
  int thisf , previ , l ;
  int aug = 0 ; /* whether it is a RM (1) or a Hamming (0) */
  int H237[23] = {1,3,7,14,12,10,8,16,20,22,18,
		  50,54,52,48,	    66,98,102,70
				,100,  96,68,64} ;
  int H246[24] = {1,3,7,5,13,9,11,15, 26,30,22,   18,28,20,24,16
		, 42,46,38,34,44,40, 36,32} ;
/*  int H143[14] = { 1,2,3,4,5,6,7,  1,2,3,4,5,6,7 } ; */
  int H143[14] = { 1,1,2,2,3,3,4,4,5,5,6,6,7,7 } ;  /* type 4 */
  int H144[14] = { 1,2,3,4,5,6,7,  1,2,3,14,15,6,7 } ; /* type 3 */
  int starti = 1 ;
 
  h->computeS = 0 ; /* default settings */
  
  if ( type == 0 ) { /* gallager constraint */
    M = 1 ; aug = 0 ; 
  } else if( type==1) { /* RM or Hamming */
    switch ( N ) {
    case( 3 ) : aug = 0 ; M = 2 ; break ; 
    case( 4 ) : aug = 1 ; M = 3 ; break ; 
    case( 7 ) : aug = 0 ; M = 3 ; break ; 
    case( 8 ) : aug = 1 ; M = 4 ; break ; 
    case( 15 ) : aug = 0 ; M = 4 ; break ; 
    case( 16 ) : aug = 1 ; M = 5 ; break ; 
    case( 31 ) : aug = 0 ; M = 5 ; break ; 
    case( 32 ) : aug = 1 ; M = 6 ; break ; 
    case( 63 ) : aug = 0 ; M = 6 ; break ; 
    case( 64 ) : aug = 1 ; M = 7 ; break ; 
    default: fprintf(stderr," sorry, I don't know this code %d - fatal error\n",N);
    }
  } else if (type==2 ) { /* shortened hamming code, expect N=12 */
    if ( N <= 14 && N >= 9 ) {
      aug = 0 ; M = 4 ;
      starti = 1 + 15 - N ;
    } 
    else
      fprintf(stderr," sorry, I don't know this code %d - fatal error\n",N);
  } else if (( type==3 ) && ( N == 14 ) ) {
      M = 4 ;
  } else if (type==3 ) { /* a rate 5/8 scheme */
    if ( N == 24 ) 
      M = 6 ;
    else if ( N == 23 ) 
      M = 7 ;
    else
      fprintf(stderr," sorry, I don't know this code %d - fatal error\n",N);
  } else if (( type==4 ) && ( N == 14 ) ) {
      M = 3 ;
  }
  K =  N - M ;   /* 11 */
  h->M = M ;   /* eg 4 , or 5 if aug'd */
  h->N = N ; 
  h->K = K ;
  I = ipower ( 2 , M ) - 1 ; /* largest number of states needed in trellis */
  h->I = I ;
  h->verbose = verbose ; 

  if ( verbose ) {
    printf ( "The (%d,%d[m=%d]) (%d) code will give me a trellis of max width %d+1\n" , N,K,M,aug,I ) ;
  }
  
  /* make space */
  pint = ivector ( 1 , N ) ;        h->pint = pint ; 
  p = cmatrix ( 1 , M , 1 , N ) ;  h->p = p ;  
  h->xo = cvector ( 1 , N ) ; 
  ff = imatrix ( 0 , N+1 , 0, I ) ;         
  bb = imatrix ( 0 , N+1 , 0, I ) ;       
  v = imatrix ( 0 , N+1 , 0, I ) ;        h->v = v ; 
  nexti = imatrix ( 0 , N+1 , 0, I ) ;      h->nexti = nexti ; 
  for ( i = 0 ; i <= N+1 ; i++ ) {
    set_ivector_const( ff[i] , 0 ,I,0 ) ;
    set_ivector_const( bb[i] , 0 ,I,0 ) ;
    set_ivector_const( v[i] , 0 ,I,0 ) ;
    set_ivector_const( nexti[i] , 0 ,I,0 ) ;
  }
  ff[0][0] = 1 + 8 ; /*these assignments assert that both transitions are possible from the first and last states - not necessarily true, but true for RM and Hmming */
  bb[N][0] = 8*8*8 + 8*8 ;
  q = dmatrix ( 1 , N , 0, 1 ) ;        h->q = q ; 
  r = dmatrix ( 1 , N , 0, 1 ) ;        h->r = r ; 
  h->po = dmatrix ( 1 , N , 0, 1 ) ;        /* not actually needed, added Tue 24/4/01 */
  if ( verbose >= 2 ) printf ("N = %d, I = %d\n" , N , I  ) ;
  h->f = dmatrix ( 0 , N+1 , 0 , I ) ;      
  h->b = dmatrix ( 0 , N+1 , 0 , I ) ;      
  zero_linear_trellis ( h ) ; 

  if ( type == 0 ) { /* plain gallager constraints */
    for ( c = 1 , i = 1 ; c <= N ; c++  ) { /* c runs across cols */
      newi = i ; 
      pint[c] = i ; 
      for ( w = 0 , d = 1 ; d <= M-aug ; d ++ ) { /* d runs down rows */
	lostbit = newi % 2 ; 
	p[d][c] = lostbit ;
	w += lostbit ;
	newi = newi / 2 ;
      }
    }    
  } else if ( type == 1 ) {
  /* here we make the parity check matrix
     of a RM or hamming 
     code by creating its columns, which are the binary integers
     plus 16 or whatever */

    for ( c = 1 , i = (aug)?0:1 ;
	  c <= N ;
	  c++ , i ++ ) { /* c runs across cols */
      newi = i ; 
      if (aug) { 
	p[M][c] = 1 ; /* aug adds an extra 1 at the bottom */
	pint[c] = i + N ;
      }    else {
	pint[c] = i ; 
      }
      for ( w = 0 , d = 1 ; d <= M-aug ; d ++ ) { /* d runs down rows */
	lostbit = newi % 2 ; 
	p[d][c] = lostbit ;
	w += lostbit ;
	newi = newi / 2 ;
      }
    }
  } else if ( type == 2 ) {
  /* shortened   RM or hamming */
    for ( c = 1 , i = starti ;
	  c <= N ;
	  c++ , i ++ ) { /* c runs across cols */
      newi = i ; 
      if (aug) { 
	p[M][c] = 1 ; /* aug adds an extra 1 at the bottom */
	pint[c] = i + N ;
      }    else {
	pint[c] = i ; 
      }
      for ( w = 0 , d = 1 ; d <= M-aug ; d ++ ) { /* d runs down rows */
	lostbit = newi % 2 ; 
	p[d][c] = lostbit ;
	w += lostbit ;
	newi = newi / 2 ;
      }
    }
  } else if ( type == 3 ) {
    if ( N == 23 ) { /* Golay codes */
      for ( c = 1 ; c <= N ;  c++ ) { /* c runs across cols */
	pint[c] = H237[c-1] ;
      }
    } else if ( N == 24 ) { /* Golay codes */
      for ( c = 1 ; c <= N ;  c++ ) { /* c runs across cols */
	pint[c] = H246[c-1] ;
      }
    } else if ( N == 14 ) {
      for ( c = 1 ; c <= N ;  c++ ) { /* c runs across cols */
	pint[c] = H144[c-1] ;
      }
    } else {
      fprintf ( stderr , "sorry can't make type %d %d - fatal\n"  , type , N ) ;
    }
    for ( c = 1 ; c <= N ;  c++ ) { /* c runs across cols */
      newi =       pint[c] ; 

      for ( w = 0 , d = 1 ; d <= M-aug ; d ++ ) { /* d runs down rows */
	lostbit = newi % 2 ; 
	p[d][c] = lostbit ;
	w += lostbit ;
	newi = newi / 2 ;
      }
    }
  } else if ( type == 4 ) {
    if ( N == 14 ) {
      for ( c = 1 ; c <= N ;  c++ ) { /* c runs across cols */
	pint[c] = H143[c-1] ;
      }
    } else {
      fprintf ( stderr , "sorry can't make type %d %d - fatal\n"  , type , N ) ;
    }
  } else {
    fprintf ( stderr , "sorry can't make type %d - fatal\n"  , type ) ;
  }
  if ( verbose  ) {  
    for ( d = 1 ; d <= M ; d ++ ) { 
      for ( c = 1 ; c <= N ; c ++ ) { 
	printf( "%d ", p[d][c] ) ; 
      }
      printf ( "\n" ) ; 
    }
    printf ( "integer names (octal)\n" ) ; 
    for ( c = 1 ; c <= N ; c ++ ) { 
      printf( "%2o ", pint[c] ) ; 
    }
    printf ( "integer names (digital)\n" ) ; 
    for ( c = 1 ; c <= N ; c ++ ) { 
      printf( "%2d ", pint[c] ) ; 
    }
    printf ( "\n" ) ; 
  }

  /* propagate forward */
  if ( verbose >=2 ) printf (" propagating all possible words (state, number)\n" ) ; 
  for ( l = 1 ; l <= N ; l ++ ) {
    i = pint[l] ;
    for ( previ = 0 ; previ <= I ; previ ++ ) {
      if ( ff[l-1][previ] ) {
	thisf = ff[l-1][previ];
	if(verbose >=2) { printf ( "%2o %-3o  " , previ , thisf ) ; }
	ff[l][previ] |= 8*8 ;
	ff[l][previ^i] |= 8*8*8; /* this node can go back to somewhere */
      }
    }
    if(verbose >=2){
      printf ("\n");
    }
  }
  
  if(verbose >=2)printf ("propagate backward\n");
  for ( l = N ; l >= 1 ; l -- ) {
    i = pint[l] ;
    for ( previ = 0 ; previ <= I ; previ ++ ) {
      thisf = bb[l][previ];
      if ( thisf ) {
	if(verbose >=2) { printf ( "%2o %-3o  " , previ , thisf ) ; }
	bb[l-1][previ] |= 1 ; /* this indicates that this node can 
				 go forward by adding zero */
	bb[l-1][previ^i] |= 8;
      }
    }
    if(verbose >=2){
      printf ("\n");
    }
  }
  
  if(verbose>=3)printf ("== list valid nodes in trellis\n" ) ; 
  if(verbose>=3)printf (" l: (i, num)\n" ) ; 
  for ( l = 0 ; l <= N ; l ++ ) {
    if(verbose>=3)printf ( "%2d: ",l) ; 
    previ = 0 ; 
    for ( i = 0 ; i <= I ; i ++ ) {
      if ( bb[l][i] && ff[l][i] ) { /* then it's valid */
	v[l][i] = bb[l][i] | ff[l][i] ;
	if(verbose>=3) { printf ( "%2o %4o  " , i , v[l][i] ) ; }
	/* go back to the last i and make it point to this */
	nexti[l][previ] = i ;
	previ = i ;
      }
    }
    /* go back to the last i and make it say end */
    nexti[l][previ] = 0 ;
    if(verbose>=3){
      printf ("\n");
    }
  }

  if(verbose)printf ("== list valid nodes in trellis\n" ) ; 
  if(verbose)printf (" l: (i, valid jumps back+forw)\n" ) ; 
  for ( l = 0 ; l <= N ; l ++ ) {
    if(verbose)printf ( "%2d: ",l) ; 
    i = 0 ; do { /* looping over all states in trellis at l */
      if(verbose) { printf ( "%2o %4o  " , i , v[l][i] ) ; }
       /* end loop */
      i = nexti[l][i] ;
    }    while ( i ) ;
    if(verbose){
      printf ("\n");
    }
  }

  /* free ff and bb */
  free_imatrix ( ff ,  0 , N+1, 0 , I ) ;
  free_imatrix ( bb ,  0 , N+1, 0 , I ) ;
}

void lt_gnu (  linear_trellis *h , const char path[] ) {
  int l ;
  int N = h->N ; 
  int  ** v = h->v ; 
  int  ** nexti = h->nexti ; 
  int verbose = 1 ; 
  int i , e = 0 ; 
  FILE *fpe , *fpo , *fpz ; 
  char edges[1000] ; 
  char ones[1000] ; 
  char zeros[1000] ; 
  int *pint = h->pint ; 
  int ii , thisv , zerook, oneok ;
  double mx,my ; 

  sprintf ( edges , "%s%s" , path , "edges" ) ; 
  sprintf ( ones  , "%s%s" , path , "ones"  ) ; 
  sprintf ( zeros , "%s%s" , path , "zeros" ) ; 
  fpe = fopen ( edges , "w" ) ; 
  fpo = fopen ( ones , "w" ) ; 
  fpz = fopen ( zeros , "w" ) ; 
  fprintf ( fpe , "# edges file written by RM.c lt_gnu\n" ) ; 
  fprintf ( fpo , "# ones  file written by RM.c lt_gnu\n" ) ; 
  fprintf ( fpz , "# zeros file written by RM.c lt_gnu\n" ) ;