r.c

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

/*   r.c         assorted minor subroutines library          release 1.1    

     Copyright   (c) 2002   David J.C. MacKay

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this library; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    GNU licenses are here :
    http://www.gnu.org/licenses/licenses.html

    Author contact details are here :
    http://www.inference.phy.cam.ac.uk/mackay/

    If you find this software useful, please feel free to make a donation to
    support David MacKay's research group.
*/
/* My minor utility routines */

#include "r.h"
#include "rand.h"

/* input */

void inputf( float *pointer )
{
  int junk ;
  scanf( "%f" , pointer ) ;
  junk = getchar( ) ;
}
void inputd( double *pointer )
{
  int junk ;
  scanf( "%lf" , pointer ) ;
  junk = getchar( ) ;
}
void inputc( unsigned char *pointer )
{
  int junk ;
  scanf( "%c" , pointer ) ;
  junk = getchar( ) ;
}
void inputi( int *pointer )
{
  int junk ;
  scanf( "%d" , pointer ) ;
  junk = getchar( ) ;
}

void inputrf( float *pointer )
{
  int junk ;
  scanf( "%f" , pointer ) ;
  do{}while( ( junk = getchar( ) ) != 10 ) ;
}
void inputrd( double *pointer )
{
  int junk ;
  scanf( "%lf" , pointer ) ;
  do{}while( ( junk = getchar( ) ) != 10 ) ;
}
void inputrc( unsigned char *pointer )
{
  int junk ;
  scanf( "%c" , pointer ) ;
  do{}while( ( junk = getchar( ) ) != 10 ) ;
}
void inputri( int *pointer )
{
  int junk ;
  scanf( "%d" , pointer ) ;
  do{}while( ( junk = getchar( ) ) != 10 ) ;
}

void clearscan( )
{
  int junk ;
  do{}while( ( junk = getchar( ) ) != 10 ) ;
}

void typeindvector( double *w , int l , int h )
{
  int i ;
  
  printf( "Please enter %d components \n" , h-l+1 ) ;
  for ( i = l ;i<= h ;i++ )
    scanf( "%lf" , &w[i] ) ;
  clearscan( ) ;
}

void constantdmatrix ( double **w , int l1 , int h1 , int l2, int h2 , 
		      double c )
{
  int i , j ;
  
  for ( i = l1 ; i <= h1 ; i++ )
    for ( j = l2 ; j <= h2 ; j++ )
      w[i][j] = c ;
}

void set_dvector_const( double *w , int l , int h , double c )
{
  int i ;
  
  for ( i = l ; i<= h ;i++ )
    w[i] = c ;
}

void set_dvector_c_dvector( double *w , int l , int h , double c , double *v )
{
  int i ;
  
  for ( i = l ; i<= h ;i++ )
    w[i] = c * v[i] ;
}

void set_ivector_const( int *w , int l , int h , int c )
{
  int i ;
  
  for ( i = l ; i<= h ;i++ )
    w[i] = c ;
}

void typeindmatrix ( double **b , int l1 , int h1 , int l2 , int h2 )
{
  int	i , j ;
  
  printf( "please type matrix\n" ) ;
  for ( i = l1 ; i<= h1-1 ; i++ )
    for ( j = l2 ; j<= h2 ; j++ )
      inputd( &b[i][j] ) ; /* This is all a trick to obtain clearscan at the last read-in */
  for ( j = l2 ; j<= h2-1 ; j++ )
    inputd( &b[i][j] ) ;
  inputrd( &b[i][j] ) ;
}

void typeincmatrix ( unsigned char **b , int l1 , int h1 , int l2 , int h2 )
{
  int	i , j , t ;
  
  printf( "please type matrix\n" ) ;
  for ( i = l1 ; i<= h1-1 ; i++ )
    for ( j = l2 ; j<= h2 ; j++ ) {
      inputi( &t ) ; b[i][j] = (unsigned char) t ; 
    }
  for ( j = l2 ; j<= h2-1 ; j++ ) {
    inputi( &t ) ; b[i][j] = (unsigned char) t ; 
  }
  inputri( &t ) ; b[i][j] = (unsigned char) t ; 
}

int ***imatrix3( int l1 , int h1 , int l2 , int h2 , int l3 , int h3 )
{
  int ***c ;
  int s1 , i ;
  
  s1 = h1-l1+1 ;
  
  c = ( int *** )malloc( ( unsigned ) s1*sizeof( int ** ) ) - l1 ;
  for ( i = l1 ;i<= h1 ;i++ )
    c[i] = imatrix( l2 , h2 , l3 , h3 ) ;
  return c ;
}

double ***dmatrix3( int l1 , int h1 , int l2 , int h2 , int l3 , int h3 )
{
  double ***c;
  int s1,i;
  
  s1=h1-l1+1;
  
  c=(double ***) malloc((unsigned) s1*sizeof(double **)) - l1;
  for (i=l1;i<=h1;i++)
    c[i]=dmatrix(l2,h2,l3,h3);
  return c;
}

long int ***limatrix3( int l1 , int h1 , int l2 , int h2 , int l3 , int h3 )
{
  long int ***c ;
  int s1 , i ;
  
  s1 = h1-l1+1 ;
  
  c = ( long int *** )malloc( ( unsigned ) s1*sizeof( long int ** ) ) - l1 ;
  for ( i = l1 ;i<= h1 ;i++ )
    c[i] = limatrix( l2 , h2 , l3 , h3 ) ;
  return c ;
}

long int **limatrix( int nrl , int nrh , int ncl , int nch )
{
  int i ;
  long int **m ;
  
  m=( long int ** )malloc( ( unsigned ) ( nrh-nrl+1 )*sizeof( long int* ) ) ;
  if ( !m ) nrerror( "allocation failure 1 in limatrix( )" ) ;
  m -= nrl ;
  
  for( i=nrl ;i<=nrh ;i++ ) {
    m[i]=( long int * )malloc( ( unsigned ) ( nch-ncl+1 )*sizeof( long int ) ) ;
    if ( !m[i] ) nrerror( "allocation failure 2 in limatrix( )" ) ;
    m[i] -= ncl ;
  }
  return m ;
}

int ipower( int a , int b )
{
  int     i=1;
  
  for(;b>0;b--)
    i*=a;
  return(i);
}

double readindmatrix 
(
 double **b ,
 int l1,
 int h1,
 int l2,
 int h2,
 char *file
)
{
	int	i, j;
	FILE    *fp;
	double	sumd=0.0;

	fp = fopen( file, "r" );
	if( !fp )   fprintf( stderr, "No such file: %s\n", file ), exit(0);

	printf( "reading in matrix from %s\n",file );
	for (i=l1; i<=h1; i++){
		for (j=l2; j<=h2; j++){
			fscanf(fp,"%lf ",&b[i][j]);
			sumd += b[i][j];
		}
	}
	fclose( fp );
	printf( "matrix in\n" );
	return(sumd);
}

int readinimatrix 
(
 int **b ,
 int l1,
 int h1,
 int l2,
 int h2,
 char *file
)
{
  int	status = 0;
  FILE    *fp;

  fp = fopen( file, "r" );
  if( !fp )   fprintf( stderr, "No such file: %s\n", file ), exit(0);

  printf( "reading in matrix from %s\n",file );
  status = fread_imatrix ( b , l1 , h1 , l2 , h2 , fp ) ; 

  fclose( fp );
  return status ;
}

int fread_imatrix 
(
 int **b ,
 int l1,
 int h1,
 int l2,
 int h2,
 FILE *fp 
)
{
  int	i, j , status = 0;

  for (i=l1; i<=h1; i++){
    /*    fprintf(stderr,"'");  fflush(stderr); */
    for (j=l2; j<=h2; j++){
      if ( fscanf(fp,"%d ",&b[i][j]) == EOF ) {
	status -- ;
	break;
      }
    }
    if ( status < 0 ) break ;
  }
  if ( status == 0 ) fprintf( stderr , "%d * %d int matrix in\n" , h1 , h2 );
  else fprintf( stderr, 
		  "Warning: readinimatrix failed at component %d\n",i);
  return status ;
}


void readindvector(double *w, int lo, int hi, char *file)
{
	int i, status = 0 ;
	FILE    *fp;
	
	fp = fopen( file, "r" );
	if( !fp )   fprintf( stderr, "No such file: %s\n", file ), exit(0);

	for (i=lo;i<=hi;i++) {
	  if ( fscanf(fp,"%lf ",&w[i]) == EOF ) {
	    status = -1 ; 
	    break ;
	  }
	}
	fclose( fp );
	if ( status < 0 ) 
	  fprintf( stderr, 
		  "Warning: readindvector failed at component %d\n",i);
}

int readdvector(double *w, int lo, int hi, char *file)
{
  int i, status = 0 ;
  FILE    *fp;
  
  fp = fopen( file, "r" );
  if( !fp ) {
    fprintf( stderr, "No such file: %s\n", file ) ;
    status -- ;
  } 
  else {
    for ( i = lo ; i <= hi ; i ++ ) {
      if ( fscanf ( fp , "%lf " , &w[i] ) == EOF ) {
	status = -1 ; 
	break ;
      }
    }
    fclose( fp );
    if ( status < 0 ) 
      fprintf( stderr, 
	      "Warning: readdvector failed at component %d\n",i);
  }
  return status ; 
}

int writedvector ( double *w, int lo, int hi, char *file)
{
  int i, status = 0 ;
  FILE    *fp;
  
  fp = fopen( file, "w" );
  if( !fp ) {
    fprintf( stderr, "writedvector---No such file: %s\n", file ) ;
    status -- ;
  } 
  else {
    for ( i = lo ; i <= hi ; i ++ ) {
      fprintf ( fp , "%g " , w[i] ) ;
    }
    fclose( fp );
  }
  return status ; 
}

int writedmatrix ( 
 double **m,
 int l1,
 int h1,
 int l2,
 int h2,
 char *file
)
{
  int status = 0 ;
  FILE    *fp;
  int i,j;
	
  fp = fopen( file, "w" );
  if( !fp ) {
    fprintf( stderr, "writedmatrix---No such file: %s\n", file ) ;
    status -- ;
  } 
  else {
    for (i=l1; i<=h1; i++){
      for (j=l2; j<=h2; j++)
	fprintf ( fp , "%g " , m[i][j] ) ;

      fprintf ( fp , "\n" ) ;
    }
    fclose( fp );
  }
  return status ; 
}

void readinivector(int *w, int lo, int hi, char *file)
{
  int status = 0 ;
  FILE    *fp;
  
  fp = fopen( file, "r" );
  if( !fp )   fprintf( stderr, "No such file: %s\n", file ), exit(0);

  status = fread_ivector ( w , lo , hi , fp ) ; 

  fclose( fp );
}

int fread_ivector(int *w, int lo, int hi, FILE *fp )
{
  int i, status = 0 ;
  
  for (i=lo;i<=hi;i++) {
    if ( fscanf(fp,"%d ",&w[i]) == EOF ) {
      status = -1 ; 
      break ;
    }
  }

/*  if ( status < 0 ) 
    fprintf( stderr, 
	    "Warning: fread_ivector failed at component %d\n",i);
*/
  return status ; 
}

int fread_dvector( double *w, int lo, int hi, FILE *fp )
{
  int i, status = 0 ;
  
  for (i=lo;i<=hi;i++) {
    if ( fscanf(fp,"%lf ",&w[i]) == EOF ) {
      status = -1 ; 
      break ;
    }
  }

/*  if ( status < 0 ) 
    fprintf( stderr, 
	    "Warning: fread_dvector failed at component %d\n",i);
*/
  return status ; 
}

int fread_cvector(unsigned char *w, int lo, int hi, FILE *fp )
{
  int i, status = 0 , bit ;
  
  for (i=lo;i<=hi;i++) {
    if ( fscanf(fp,"%d", &bit ) == EOF ) {
      status = -1 ; 
      break ;
    }
    w[i] = (unsigned char) bit ;
  }

/*  if ( status < 0 ) 
    fprintf( stderr, 
	    "Warning: fread_cvector failed at component %d\n",i);
*/
  return status ; 
}

void printoutimatrix 
(
 int **m,
 int l1,
 int h1,
 int l2,
 int h2
)
{
  write_imatrix ( stdout , m , l1 , h1 , l2 , h2 ) ; 
}

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

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


void allocate_cm_inversion 
( unsigned char **mo,
 int l,
 int N,
 unsigned char **mi,
 cm_inversion *p
)
{
  int i,j ; 

  p->l = l ; p->N = N ; p->mo = mo ; p->mi = mi ; 
  p->m = cmatrix ( l , N , l , N ) ; 
  p->mt=(unsigned char **)malloc((unsigned) (N-l+1)*sizeof(unsigned char*));
  p->mt -= l ;
  p->perm = ivector ( l , N ) ; 
  p->iperm = ivector ( l , N ) ; 
  p->i = l ; 
  for (i=l; i<=N; i++){
    p->perm[i] = l-1 ;  /* deliberately off scale */
    p->iperm[i] = l-1 ;
    for (j=l; j<=N; j++) {
      mi[i][j] = 0 ; 
      p->m[i][j] = mo[i][j] ;
    }
    mi[i][i] = 1 ;
  }
}

void print_cm_inversion 
( cm_inversion *p
)
{
/*  printf( "mo:\n" ) ; 
  printoutcmatrix (  p->mo , p->l , p->N , p->l , p->N ) ; */
  printf( "m:\n" ) ; 
  printoutcmatrix (  p->m , p->l , p->N , p->l , p->N ) ; 
  printf( "mi:\n" ) ; 
  printoutcmatrix (  p->mi , p->l , p->N , p->l , p->N ) ; 
  printf( "i = %d\n" , p->i ) ;
  printf( "perm:\n" ) ; 
  printoutivector ( p->perm , p->l , p->N ) ; 
  printf( "iperm:\n" ) ; 
  printoutivector ( p->iperm , p->l , p->N ) ; 
  pause_for_return () ;
}

void free_cm_inversion 
( cm_inversion *p
)
{
  free_ivector ( p->perm , p->l , p->N ) ; 
  free_ivector ( p->iperm , p->l , p->N ) ; 
  free_cmatrix ( p->m , p->l  , p->N , p->l , p->N ) ;
  free( ( char * ) ( p->mt + p->l ) ) ;
}

/* invert_cmatrix attempts to invert p->m, or to continue to invert it, 
   starting from the state indicated by i,j,perm,etc. 
   if invert_cmatrix succeeds it returns 1, else 0, with p->i 
   recording the row at which non-full-rank was detected . 
*/
int invert_cmatrix
(
 cm_inversion *p
)
{
  int j , done , i2 , k ; 
  unsigned char *irow ;

  for ( ; p->i <= p->N ; p->i ++ ) {
    irow = p->m[p->i] ;
    /* search for a 1 in this row */
    for ( j = p->l , done = 0 ; j <= p->N ; j++ ) {
      if ( irow[j] &&            /* if we find a 1              and */
	  p->perm[j] < p->l ) {  /* if j is not already on the list */
	p->perm[j]  = p->i ; 
	p->iperm[p->i] = j ; 
	done = 1 ;
	break ; /* break out of j loop */
      }
    }
    if ( !done ) {
      return ( 0 ) ;
      break ;
    }
    else { /* search for 1s in subsequent rows. */
      for ( i2 = p->i + 1 ; i2 <= p->N ; i2 ++ ) {
	if ( p->m[i2][j] ) {
	  for ( k = p->l ; k <= p->N ; k ++ ) {
	    p->m[i2][k] ^= irow[k] ;
 	    p->mi[i2][k] ^= p->mi[p->i][k] ; 
	  }
	}
      }
    }
  }      
  /* at this point, if done==1 then m contains a (permuted) upper
     triangular matrix and mi contains a (permuted) lower triangular
     matrix. 
     */
/*  printf ("lu decomposition achieved\n" ) ; */
/*  print_cm_inversion ( p ) ; */
/*  printf ("now working back\n" ) ;  */
  done = invert_utriangularc ( p ) ;