r.c

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

  /* finally, undo the permutation - make mi[l] = mi[perm[l]] etc. */
  undo_cm_perm ( p ) ; 
  return ( done ) ; 
}

void undo_cm_perm ( cm_inversion *p ) {
  int i , pi ;
  
  for ( i = p->N ; i >= p->l ; i-- ) {
    p->mt[i] = p->mi[i] ;                 /* copy all columns */
  } 
  for ( i = p->N ; i >= p->l ; i-- ) {    /* reallocate them */
    pi = p->perm[i] ;
    if ( ( pi >= p->l ) && ( pi <= p->N ) )
      p->mi[i] = p->mt[pi] ;
    else {
      fprintf ( stderr , "error in undo_cm_perm, perm[%d] = %d\n" , i , 
	       pi ) ; 
      exit ( 0 ) ; 
    }
  } 
}

/* invert_ltriangularc takes p->m and p->mi 
   and applies a series of row sums to both 
   such that p->m becomes the identity.
   If p->mi started out as the identity, it becomes the inverse 
   of p->m; otherwise it becomes (that inverse)*(whatever it was before) */
int invert_utriangularc 
( cm_inversion *p ) {
  int i , vj , j  , k ;
   
  for ( i = p->N ; i >= p->l ; i-- ) {
    for ( vj = i + 1 ; vj <= p->N ; vj++ ) { /* virtual j */
      j = p->iperm[vj] ;
      if (!( ( j >= p->l ) && ( j <= p->N ) ))
      {
	fprintf ( stderr , "error in invert_utriangularc, iperm[%d] = %d\n" ,
		 vj ,  j ) ; 
	exit ( 0 ) ; 
      }
      if ( p->m[i][j] ) {
	for ( k = p->l ; k <= p->N ; k ++ ) {
	  p->m[i][k] ^= p->m[vj][k] ; 
	  p->mi[i][k] ^= p->mi[vj][k] ; 
	} 
/*	printf("%d,%d -> \n",i,j);
	print_cm_inversion ( p ) ; */
      }
    }
  }
  return 1 ; /* inverse has been evaluated */
}


int modify_cmatrix_row ( cm_inversion *p ) {
  /* proceed along the current row seeking an available j */
  unsigned char *irow ;
  int j , done = 0 ;

  irow = p->m[p->i] ;
  /* search for a 1 in this row */
  for ( j = p->l ; done == 0 && j <= p->N ; j++ ) {
    if ( p->perm[j] < p->l ) {  /* if j is not already on the list */
      irow[j] = 1 ; 
      p->mo[p->i][j] ^= 1 ; /* modify original matrix */
      done = 1 ;
    }
  }
  return ( done ) ; /* if done == 0 then this 
		       routine has failed to
		       find a j; an unexpected
		       outcome */
}

/* 
   pbm format is 
   P1
   30 90
   (30 rows with 90 cols)

   which gives (under xv) a picture 
   with 90 rows and 30 cols.

   So in fact rows and columns are interchanged
*/
void cmatrix2pbm
(
 unsigned char **m,
 int l1,
 int h1,
 int l2,
 int h2,
 FILE *fp 
)
{
  int i,j;
  
  fprintf( fp , "P1\n%d %d\n" , h1 - l1 + 1 , h2 - l2 + 1 ) ;
  for (i=l1; i<=h1; i++){
    for (j=l2; j<=h2; j++)
      fprintf( fp , "%d ",m[i][j]);
    fprintf( fp , "\n");
  }                   
}

void printoutcmatrix 
(
 unsigned char **m,
 int l1,
 int h1,
 int l2,
 int h2
)
{
  int i,j;
  
  for (i=l1; i<=h1; i++){
    if ( !(i%10) ) printf ( "| " ) ; 
    else printf ( "  " ) ; 
  }
  printf("\n");
  for (i=l1; i<=h1; i++){
    for (j=l2; j<=h2; j++)
      printf("%d ",m[i][j]);
    printf("\n");
  }                   
}

void printoutcmatrix1 
(
 unsigned char **m,
 int l1,
 int h1,
 int l2,
 int h2
)
{
  int i,j;
  
  for (i=l1; i<=h1; i++){
    for (j=l2; j<=h2; j++)
      if ( m[i][j] ) printf("1") ; else printf(" ") ;
    printf("\n");
  }                   
}

void printoutivector
(
 int *m,
 int l1,
 int h1
)
{
  int i;
  
  for (i=l1; i<=h1; i++)
    printf( "%d " , m[i] );
  printf("\n");
}

void write_ivector
(
 FILE *fp,
 int *m,
 int l1,
 int h1
)
{
  int i;
  
  for (i=l1; i<=h1; i++)
    fprintf( fp , "%d " , m[i] );
  fprintf( fp , "\n");
}

void write_cvector
(
 FILE *fp ,
 unsigned char *m,
 int l1,
 int h1
)
{
  int i;
  
  for (i=l1; i<=h1; i++)
    fprintf( fp , "%d\n" , m[i] );
/*  fprintf( fp ,"\n"); */
}

void printoutcvector
(
 unsigned char *m,
 int l1,
 int h1
)
{
  int i;
  
  for (i=l1; i<=h1; i++)
    printf( "%d " , m[i] );
  printf("\n");
}

void printoutcvector1
(
 unsigned char *m,
 int l1,
 int h1
)
{
  int i;
  
  printf("|");
  for (i=l1; i<=h1; i++)
    if ( m[i] ) printf( "1" ) ; else printf( " " );
  printf("|\n");
}

void printoutdmatrix 
(
 double **m,
 int l1,
 int h1,
 int l2,
 int h2,
 int style
)
{
	int i,j;
	
	for (i=l1; i<=h1; i++){
		for (j=l2; j<=h2; j++)
        	 	pd(m[i][j],style);
          	printf("\n");
     	}                   
}

void printoutdmatrix3 
(
 double ***m,
 int l1,
 int h1,
 int l2,
 int h2,
 int l3,
 int h3,
 int style
)
{
  int i,j,k;
	
  for (i=l1; i<=h1; i++){
    for (j=l2; j<=h2; j++) {
      for (k=l3; k<=h3; k++)
	pd(m[i][j][k],style);
      printf("\n");
    }                   
    printf("\n");
  }                   
}

void pd
(
 double w,
 int p
) 	/* prints a single double with a specified number of digits */ 
{
	switch(p){
		case(0):		/* Don't print */
			break;
		case(-1):		/* Standard format */
			printf("%g ",w);
			break;
		case(1):
			printf("%1.0f ",w);
			break;
		case(100):
			printf("%1.0f ",w*10);
			break;
		case(11):
			if (w==1.0)printf("%1.0f ",w);
			else printf("  ");
			break;
		case(10): /* 0/1 <- -1/1 */
			printf("%1.0f ",(double)(w+1.0)*0.5);
			break;
		case(2):
			printf("%2.0f ",w);
			break;
		case(29):
			printf("%2f ",w);
			break;
		case(21):
			printf("%2.1f ",w);
			break;
		case(20):
			printf("%2.0f ",w);
			break;
		case(200):
			printf("%2.0f",w*100);
			break;
		case(3):
			printf("%3.1f ",w);
			break;
		case(39):
			printf("%3f ",w);
			break;
		case(30):
			printf("%3.0f ",w);
			break;
		case(300):
			printf("%3.0f ",w*1000);
			break;
		case(31):
			printf("%3.1f ",w);
			break;
		case(32):
			printf("%3.2f ",w);
			break;
		case(4):
			printf("%4.1f ",w);
			break;
		case(49):
			printf("%4f ",w);
			break;
		case(40):
			printf("%4.0f ",w);
			break;
		case(400):
			printf("%4.0f ",w*10000);
			break;
		case(41):
			printf("%4.1f ",w);
			break;
		case(43):
			printf("%4.3f ",w);
			break;
		case(42):
			printf("%4.2f ",w);
			break;
		case(52):
			printf("%5.2f ",w);
			break;
		case(53):
			printf("%5.3f ",w);
			break;
		case(54):
			printf("%5.4f ",w);
			break;
		case(50): 
			printf("%5.0f ",w);
			break;
		case(500): /* for probabilities */
			printf("%5.0f ",w*100000);
			break;
		case(60): 
			printf("%6.0f ",w);
			break;
		case(62):
			printf("%6.2f ",w);
			break;
		case(63):
			printf("%6.3f ",w);
			break;
		case(64):
			printf("%6.4f ",w);
			break;
		case(72):
			printf("%7.2f ",w);
			break;
		case(74):
			printf("%7.4g ",w);
			break;
		case(84):
			printf("%8.4g ",w);
			break;
		case(94):
			printf("%+9.4g ",w);
			break;
		case(76):
			printf("%7.6f ",w);
			break;
		case(5):
			printf("%5g ",w);
			break;
		case(6):
			printf("%6g ",w);
			break;
		case(600):
			printf("%6.0f ",w*1000000);
			break;
		case(7):
			printf("%7g ",w);
			break;
		case(700):
			printf("%7.0f ",w*10000000);
			break;
		default:
			printf("Error in pd rule \n");
			break;
	}
}

void pdv
(
 double *w, int m, int n, int p
 )	/* prints dvector without newline, with specified number of digits */
{
	int i;
	for(i=m;i<=n;i++)
		pd(w[i],p);
}
 

void 	pause_for_return()
{
	int c;

	printf( "press return to continue" );
	do{}while((c=getchar())!=10);
}

double gammln(	      double xx) /* log of gamma function */
{
  double x,tmp,ser;
  static double cof[6]={76.18009173,-86.50532033,24.01409822,
                -1.231739516,0.120858003e-2,-0.536382e-5};
  int j;
  
  x=xx-1.0;
  tmp=x+5.5;
  tmp -= (x+0.5)*log(tmp);
  ser=1.0;
  for (j=0;j<=5;j++) {
    x += 1.0;
    ser += cof[j]/x;
  }
  return -tmp+log(2.50662827465*ser);
}

double digamma_dif ( double s , double c , int exact ) {
/* returns Psi(s) - Psi(c) [approx, good for s>=c+1] */
  int F = (int)(s+0.01 - c) ; /* recover the integer (round-down assumed) */ 
  double ans = 0.0 ; 
  if ( F==0 || s==c || fabs(s-c) < 0.01 ) return ( 0.0 ) ; 
  else if ( !exact )
    return ( ( 1.0 / c ) + log ( ( s - 0.5 ) / ( c + 0.5 ) ) ) ; 
  else {
    for ( s -- ; F >= 1 ; F -- , s-- )   /* assume that F really is an integer */ {
      ans += 1.0 / s ;  /* sum f=0..F-1   1/(f+u)  */
    }
    return ans ; 
  }
}

double digamma_dif_deriv ( double s , double c  ) {
  int F = (int)(s+0.01 - c) ; /* recover the integer (round-down assumed) */ 
  double ans = 0.0 ; 
  if ( F==0 || s==c || fabs(s-c) < 0.01 ) return ( 0.0 ) ; 
  else {
    for ( s -- ; F >= 1 ; F -- , s-- )   /* assume that F really is an integer */ {
      ans -= 1.0 / ( s * s ) ;  /* sum f=0..F-1   1/(f+u)  */
    }
    return ans ; 
  }
}

int find_rank ( double act , double *v , int lo , int hi ) {
  /* given an ordered list v[lo] < v[lo+1] < ... < v[hi], 
     and assuming that act > v[lo], the task is to return 
     the ranking of act, i.e. the integer place where it will push 
     in. For example, if we find act < v[lo+1], the answer is lo. 
     If act > v[hi], answer is hi. 
     Note no comparison with v[lo] should occur, we already know it's
     bigger. 
     */

  int query ;

  if ( lo > hi ) { 
    fprintf ( stderr , "Eek, lo > hi ???\n" ) ; 
    exit ( 0 ) ; 
  }
  if ( lo == hi ) { /* we are done */
    return ( lo ) ; 
  } 
  else { 
    query = lo + ( hi + 1 - lo ) / 2 ; 
    if ( query == lo ) { 
      fprintf ( stderr , "Eek, query == lo ???\n" ) ; 
      exit ( 0 ) ; 
    }
    else if ( act > v[query] ) { 
      return ( find_rank ( act , v , query , hi ) ) ;
    }
    else {
      return ( find_rank ( act , v , lo , query - 1 ) ) ;
    }
  }
} 

int dotprod_mod2 ( int *a , int *b , int lo , int hi ) {
  int i , ans = 0 ; 
  for ( i = lo ; i <= hi ; i ++ ) {
    ans ^= a[i] * b[i] ;
  }
  return ans ;
}
int idotprod ( int *a , int *b , int lo , int hi ) {
  int i , ans = 0 ; 
  for ( i = lo ; i <= hi ; i ++ ) {
    ans += a[i] * b[i] ;
  }
  return ans ;
}

unsigned char cdotprod_mod2 ( unsigned char *a , unsigned char *b , int lo , int hi ) {
  unsigned char ans = 0 ; 
  int i ; 
  for ( i = lo ; i <= hi ; i ++ ) {
    ans ^= a[i] & b[i] ;
  }
  return ans ;
}

void mult_cms 
( unsigned char **A , unsigned char **B , unsigned char **C , int l , int N ) {
  int i , j , k ; 
  unsigned char a ;

  for ( i = l ; i <= N ; i ++ ) {
    for ( k = l ; k <= N ; k ++ ) {
      a = 0 ; 
      for ( j = l ; j <= N ; j ++ ) {
	a ^= A[i][j] & B[j][k] ;
      }
      C[i][k] = a ;
    }
  }
}

void mult_cm_cv 
( unsigned char **A , unsigned char *b , unsigned char *c , int l , int N ) {
  int i , j ; 
  unsigned char a ;

  for ( i = l ; i <= N ; i ++ ) {
    a = 0 ; 
    for ( j = l ; j <= N ; j ++ ) {
      a ^= A[i][j] & b[j];
    }
    c[i] = a ;
  }
}

int cdotprod ( unsigned char *a , unsigned char *b , int lo , int hi ) {
  int i , ans = 0 ; 
  for ( i = lo ; i <= hi ; i ++ ) {
    ans += a[i] & b[i] ;
  }
  return ans ;
}

/* For useful routines for dealing with triangular char matrices and 
   with lists of matrices of the form
   U = 10
       11
   (just one off diagonal 1)
   see
   ~/code/old/mnc.c.940106
*/


#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)

float ran3(
	   int *idum
)
{
	static int inext,inextp;
	static long ma[56];
	static int iff=0;
	long mj,mk;
	int i,ii,k;

	if (*idum < 0 || iff == 0) {
		iff=1;
		mj=MSEED-(*idum < 0 ? -*idum : *idum);
		mj %= MBIG;
		ma[55]=mj;
		mk=1;
		for (i=1;i<=54;i++) {
			ii=(21*i) % 55;
			ma[ii]=mk;
			mk=mj-mk;
			if (mk < MZ) mk += MBIG;
			mj=ma[ii];
		}
		for (k=1;k<=4;k++)
			for (i=1;i<=55;i++) {
				ma[i] -= ma[1+(i+30) % 55];
				if (ma[i] < MZ) ma[i] += MBIG;
			}
		inext=0;
		inextp=31;
		*idum=1;
	}
	if (++inext == 56) inext=1;
	if (++inextp == 56) inextp=1;
	mj=ma[inext]-ma[inextp];
	if (mj < MZ) mj += MBIG;
	ma[inext]=mj;
	return mj*FAC;
}

#undef MBIG
#undef MSEED
#undef MZ
#undef FAC