cmatrix.c

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

    uu = a->num_mlist[mm] ; 
    if ( uu > tbm ) tbm = uu ; 
  }
  a->tot = 0 ; 
  for ( nn  = 1 ; nn <= a->N ; nn++ ) {
    uu = a->num_nlist[nn] ; 
    if ( uu > tbn ) tbn = uu ; 
    a->tot += uu ;
  }
  if ( tbm >= a->biggest_num_m_alloc ) {
    fprintf (stderr , "eek m list overrun %d %d\n" ,
	     tbm , a->biggest_num_m_alloc ) ;
    status -- ; 
  }
  else a->biggest_num_m = tbm ; 
  if ( tbn >= a->biggest_num_n_alloc ) {
    fprintf ( stderr , "eek n list overrun %d %d\n" , 
	     tbn , a->biggest_num_n_alloc ) ;
    status -- ; 
  }
  else a->biggest_num_n = tbn ; 
  return status ; 
}

void report_alist ( alist_matrix *a , int verbose ) {
  int m , n ; 

  printf ( "nlist numbers :\n" ) ;
  printoutivector ( a->num_nlist , 1 , a->N ) ; 
  printf ( "mlist numbers :\n" ) ;
  printoutivector ( a->num_mlist , 1 , a->M ) ; 
  if ( verbose >= 2  )   {
    printf ( "nlist :\n" ) ;
    for ( n = 1 ; n <= a->N ; n++ ) {
      printf ( "%3d : " , n  ) ;
      printoutivector ( a->nlist[n] , 1 , a->num_nlist[n] ) ; 
    }
    printf ( "\n" ) ; 
    printf ( "mlist :\n" ) ;
    for ( m = 1 ; m <= a->M ; m++ ) {
      printf ( "%3d : " , m  ) ;
      printoutivector ( a->mlist[m] , 1 , a->num_mlist[m] ) ; 
    }
    printf ( "\n" ) ; 
  }
}
 

void sparse_random_cmatrix
( unsigned char **A , int per_row , int NN ) {
  int i , j , num_left  ;

  for ( i = 1 ; i <= NN ; i++ )  {
    j = NN ;
    num_left = per_row ;
    for ( ; j >= 1 ; j-- ) {
      if ( num_left > 0 && ( ranf() <=  (double) num_left / (double) j ) ) {
	A[i][j] = 1 ;
	num_left -- ;
      }  else A[i][j] = 0 ;
    }
  }
}

void sparse_random_cmatrix2 /* this ensures that number_per_col equal too */
( unsigned char **A , int per , int NN ) {
  int i , j , num_left  ;
  int *new_num_c , *num_c ; /* number left in column */
  
  new_num_c = ivector ( 1 , NN ) ; 
  num_c = ivector ( 1 , NN ) ; 
  for ( j = 1 ; j <= NN ; j++ )  
    num_c[j] = per ; 
  for ( i = NN ; i >= 1 ; i-- ) {
    do { /* generate the next row at random until the number in the row is
	    right */
      j = NN ;
      num_left = per ;
      for ( ; j >= 1 ; j-- ) {
	if ( num_c[j] > 0 && ( ranf() <=  (double) num_c[j] / (double) i ) ) {
	  A[i][j] = 1 ;
	  num_left -- ;
	  new_num_c[j] = num_c[j] - 1 ; 
	}  else {
	  A[i][j] = 0 ;
	  new_num_c[j] = num_c[j] ;
	}
      }
      if ( num_left == 0 ) { /* move on */
	for ( j = 1 ; j <= NN ; j++ )  
	  num_c[j] = new_num_c[j] ;
      }
/*      printf ( "%d %d; " , i , num_left ) ; fflush(stdout) ;  */
    } while ( num_left != 0 ) ;
  }
  free_ivector ( new_num_c , 1 , NN ) ;
  free_ivector ( num_c , 1 , NN ) ;
}

void sparse_rectangular_cmatrix ( unsigned char **A , int per , int M , int N , int hd ){
/* makes a sparse rectangular matrix with per per column */
/* and number per row as equal as possible 
   if hd > 0 then all columns must be more than hd distance from each other.
   if a column is hd or closer to another column, then a "h" is printed
   and we restart that column */
  int i , j , num_left  ;
  int *num_r ; /* number left in row */
  int *mylist ; 
  int per_row , extras , hdfail = 0 ;
  int dist = 0  , col , m ;
  double di ;

  /* N columns, M rows */

  per_row = ( per * N ) / M ;  
  extras = per * N - per_row * M ;
  if ( extras < 0 ) fprintf ( stderr , "eek unexpected extras < 0 %d\n" , extras )  ;

  num_r = ivector ( 1 , M ) ; 
  mylist =  ivector ( -1 , per ) ; 

  printf ( "Start:" ) ; fflush(stdout) ;

  for ( j = 1 ; j <= M ; j++ )  
    num_r[j] = per_row ; 
  if ( extras > 0 ) {
    fprintf ( stderr , "Adding %d extras because inexact division\n" , extras ) ; 
    for ( j = M ; j >= 1 && ( extras > 0 ); j-- )  {
      if ( ranf() <=  (double) extras / (double) j ) {
	num_r[j] ++ ;
      }
    }
  }

  for ( i = N ; i >= 1 ; i-- ) {
    di = 1.0 / (double) i ; 
    printf ( "%d;" , i ) ; fflush(stdout) ; 
    do { /* generate the next col at random until the number in the col is
	    right 
	    if numleft has got down to -1, then we can certainly stop
	    and start again */
      j = M ;
      num_left = per ;
      for ( ; j >= 1 && (num_left > -1)  ; j-- ) {
	if ( num_r[j] > 0 && ( ranf() <=  (double) num_r[j] * di ) ) {
	  mylist[num_left] = j ; 
	  num_left -- ;
	} 
      }
      if ( num_left == 0 ) { /* check hd then move on */
	if ( hd > 0 ) {
	  hdfail = 0 ; 
	  for ( col = N ; (col > i) && ( hdfail == 0 ) ; col-- ) {
	    dist = 0 ; 
	    for (  j = 1 ; j <= per && ( dist <= hd ) ; j ++ ) {
	      m = mylist[j] ;
	      dist += A[m][col] ? 0 : 2 ; 
	    }
	    if ( dist <= hd ) { hdfail = 1 ; }
	  }
	}
	if ( hd > 0 && hdfail == 1 ) {
	  printf ( "h" ) ; fflush(stdout) ;
	  if ( i == 1 ) { /* tough luck, have to accept the failure */
	    fprintf ( stderr, "\nwarning, final hamming check failure %d\n" , dist ) ;
	    hdfail = 0 ; 
	  }
	} 
	if ( hdfail == 0 ) {
	  /* final acts before moving on */
	  for ( j = 1 ; j <= per ; j ++ ) {
	    m = mylist[j] ;
	    A[m][i] = 1 ;
	    num_r[m] -- ; 
	  }
	}

      } else {
	printf ( "." ) ; /* fflush(stdout) ;   */
      }
    } while ( num_left != 0 || hdfail != 0 ) ;
  }
  free_ivector ( mylist , -1 , per ) ;
  free_ivector ( num_r , 1 , M ) ;
}

int sparse_rectangular_cmatrix2 ( unsigned char **A , int c0 , int per , int M , int N , int hd ){
/* Two differences from 1st routine.
   The matrix is constructed by ROW, from bottom to top.
   And a start can be made a distance c0 from the bottom if desired 
*/
/*
   if hd > 0 then all rows must be more than hd distance from each other.
   if a row is hd or closer to another row, then a "h" is printed
   and we restart that row 

   This version starts from row N - c0, and compares with all previous
   rows right up to N
*/
  int i , j , num_left  ;
  int *num_r ; /* number left in col */
  int *mylist ; 
  int per_col , extras , hdfail = 0 ;
  int dist = 0 , row , m ;
  double di ;
  int num_rows_to_do = N - c0 ; 
  int status = 0 ; 

  /* N rows, M cols */

  per_col = ( per * num_rows_to_do ) / M ;  
  extras = per * num_rows_to_do - per_col * M ;
  if ( extras < 0 ) fprintf ( stderr , "eek unexpected extras < 0 %d\n" , extras )  ;

  num_r = ivector ( 1 , M ) ; 
  mylist =  ivector ( -1 , per ) ; 

  printf ( "Start:" ) ; fflush(stdout) ;

  for ( j = 1 ; j <= M ; j++ )  
    num_r[j] = per_col ; 
  if ( extras > 0 ) {
    fprintf ( stderr , "Adding %d extras because inexact division\n" , extras ) ; 
    for ( j = M ; j >= 1 && ( extras > 0 ); j-- )  
      if ( ranf() <=  (double) extras / (double) j ) 
	num_r[j] ++ ;
  }

  for ( i = num_rows_to_do ; i >= 1 ; i-- ) {
    di = 1.0 / (double) i ; 
    printf ( "%d;" , i ) ; fflush(stdout) ; 
    do { /* generate the next row at random until the number in the row is
	    right */
      j = M ;
      num_left = per ;
      for ( ; j >= 1 && (num_left > -1)  ; j-- ) {
	if ( num_r[j] > 0 && ( ranf() <=  (double) num_r[j] * di ) ) {
	  mylist[num_left] = j ; 
	  num_left -- ;
	} 
      }
      if ( num_left == 0 ) { /* check hd then move on */
	if ( hd > 0 ) {
	  hdfail = 0 ; 
	  for ( row = N ; (row > i) && ( hdfail == 0 ) ; row-- ) {
	    dist = 0 ; 
	    for (  j = 1 ; j <= per && ( dist <= hd ) ; j ++ ) {
	      m = mylist[j] ;
	      dist += A[row][m] ? 0 : 2 ; /* worst case - assume any 1
					     implies another difference
					     somewhere; correct if both
					     rows have same weight */
	    }
	    if ( dist <= hd ) { hdfail = 1 ; }
	  }
	}
	if ( hd > 0 && hdfail == 1 ) {
	  printf ( "h" ) ; fflush(stdout) ;
	  if ( i == 1 ) { /* tough luck, have to accept the failure */
	    fprintf ( stderr, "\nwarning, final hamming check failure %d\n" , dist ) ;
	    hdfail = 0 ; 
	    status -- ; 
	  }
	} 
	if ( hdfail == 0 ) {
	  /* final acts before moving on */
	  for ( j = 1 ; j <= per ; j ++ ) {
	    m = mylist[j] ;
	    A[i][m] = 1 ;
	    num_r[m] -- ; 
	  }
	}

      } else {
	printf ( "." ) ; /* fflush(stdout) ;  */
      }
    } while ( num_left != 0 || hdfail != 0 ) ;
  }
  free_ivector ( mylist , -1 , per ) ;
  free_ivector ( num_r , 1 , M ) ;
  return status ; 
}

int old_sparse_rectangular_alist ( alist_matrix *a , int c0 , int c1 , int per , int Cols , int Rows , int hd ){
/* 
   Constructs a matrix by row, from bottom to top.
   A start can be made a distance c0 from the bottom if desired 
   A stop can be called after only c1 rows from the bottom have been 
   finished also.  

   if hd > 0 then all rows must be more than hd distance from each other.
   if a row is hd or closer to another row, then a "h" is printed
   and we restart that row 

   This version starts from row Rows - c0, and compares with all previous
   rows right up to Rows
*/
  int status = 0 ; 
  int nn , u , latest ; 
  int i , j , num_left  ;
  int *num_r ; /* number left in col */
  int *mylist ; 
  int per_col , tot_per_col , extras , hdfail = 0 ;
  int m ;
  double di ;
  int num_rows_to_do = c1 - c0 ; 
  int *relatives ; 
  int relmax , rr , ss , prevrr ;
  int *num_nlist = a->num_nlist ; 
  int **nlist = a->nlist ; 

  /* N rows, M cols */

  per_col = ( per * num_rows_to_do ) / Cols ;  
  extras = per * num_rows_to_do - per_col * Cols ;
  if ( extras < 0 ) fprintf ( stderr , "eek unexpected extras < 0 %d\n" , extras )  ;

  num_r = ivector ( 1 , Cols ) ; 
  mylist =  ivector ( -1 , per ) ; 

  tot_per_col = ( per * Rows ) / Cols ;  
  relmax = per * ( tot_per_col + 2 ) ; /* was tot_per_col + 1 */
  relatives = ivector ( 1 ,  relmax + 10 ) ; /* was relmax */
  
  printf ( "Start:" ) ; fflush(stdout) ;

  for ( j = 1 ; j <= Cols ; j++ )  
    num_r[j] = per_col ; 
  if ( extras > 0 ) {
    fprintf ( stderr , "Adding %d extras because inexact division\n" , extras ) ; 
    for ( j = Cols ; j >= 1 && ( extras > 0 ); j-- )  
      if ( ranf() <=  (double) extras / (double) j ) 
	num_r[j] ++ ;
  }

  for ( i = num_rows_to_do , m = Rows - c0 ; i >= 1 ; i-- , m-- ) {
    di = 1.0 / (double) i ; 
    printf ( "%d;" , i ) ; fflush(stdout) ; 
    do { 
      j = Cols ;
      num_left = per ;
      for ( ; j >= 1 && (num_left > -1)  ; j-- ) {
	if ( num_r[j] > 0 && ( ranf() <=  (double) num_r[j] * di ) ) {
	  mylist[num_left] = j ; 
	  num_left -- ;
	} 
      }
      if ( num_left == 0 ) { /* check hd then move on */
	if ( hd > 0 ) {
	  hdfail = 0 ; 
/* 
   for each j on the list
   go and find any other rows that share this j
   check if we have already met this row today. If so, then that's 
   a failure for the traditional setting of hd
*/
	  rr = 0 ; 
	  prevrr = 0 ; 
	  for ( j = 1 ; j <= per && ( hdfail == 0 ) ; j ++ ) {
	    nn = mylist[j] ;
	    for ( u = num_nlist[nn] ; u >= 1 && ( hdfail == 0 ) ; u-- ) {
	      rr ++ ;
	      latest = relatives[rr] = nlist[nn][u] ; 
	      for ( ss = prevrr ; ss >= 1 ; ss -- ) {
		if ( relatives[ss] == latest ) {
		  hdfail = 1 ; 
		  break ; 
		}
	      }
	    }
	    prevrr = rr ; 
	    if ( rr > relmax ) {
	      fprintf ( stderr, "\nwarning, relatives overran vector %d %d\n" , rr , relmax ) ;
	      status -- ; 
	    }
	  }
	}
	if ( hd > 0 && hdfail == 1 ) {
	  printf ( "h" ) ; fflush(stdout) ;
	  if ( i == 1 ) { /* tough luck, have to accept the failure */
	    fprintf ( stderr, "\nwarning, final hamming check failure at relative %d %d out of %d\n" , rr , ss , relmax ) ;
	    hdfail = 0 ; 
	  }
	} 
	if ( hdfail == 0 ) {
	  /* final acts before moving on */
	  for ( j = 1 ; j <= per ; j ++ ) {
	    nn = mylist[j] ;
	    add_to_alist ( a , nn , m ) ; 
	    num_r[nn] -- ; 
	  }
	}
      } else {
	printf ( "." ) ; /* fflush(stdout) ;  */
      }
    } while ( num_left != 0 || hdfail != 0 ) ;
  }
  free_ivector ( mylist , -1 , per ) ;
  free_ivector ( num_r , 1 , Cols ) ;
  free_ivector ( relatives , 1 , relmax ) ;
  
  return status ; 
}

int sparse_rectangular_alist ( alist_matrix *a , int c0 , int c1 , int per , int Cols , int Rows , int hd ){
/* 
   Constructs a matrix by row, from bottom to top.
   A start can be made a distance c0 from the bottom if desired 
   A stop can be called after only c1 rows from the bottom have been 
   finished also.  

   if hd > 0 then all rows must be more than hd distance from each other.
   if a row is hd or closer to another row, then a "h" is printed
   and we restart that row 

   This version starts from row Rows - c0, and compares with all previous
   rows right up to Rows

   added hfailmax 9809, but then noticed that uneven_sparse....
 has this attempts feature in it already. Use code5?
*/
  int status = 0 ; 
  int hfailmax = 200 ; /* max permitted attempts  per row (set to zero to disable this rule) */
  int hfails = 0 ; 
  int nn , u , latest ; 
  int i , j , num_left  ;
  int *num_r ; /* number left in col */
  int *pikd_r ;   /* whether this col pikd this time */
  int *cum ;   /* cumulative probability */
  int *mylist ; 
  int per_col , tot_per_col , extras , hdfail = 0 ;
  int m ;
  double di ;
  int num_rows_to_do = c1 - c0 ; 
  int *relatives ; 
  int relmax , rr , ss , prevrr ;
  int *num_nlist = a->num_nlist ; 
  int **nlist = a->nlist ; 
  double un ; /* random 0,1 */

  /* N rows, M cols */

  per_col = ( per * num_rows_to_do ) / Cols ;  
  extras = per * num_rows_to_do - per_col * Cols ;
  if ( extras < 0 ) fprintf ( stderr , "eek unexpected extras < 0 %d\n" , extras )  ;

  num_r = ivector ( 1 , Cols ) ; 
  pikd_r = ivector ( 1 , Cols ) ; 
  cum = ivector ( 0 , Cols ) ; cum[0] = 0 ; 
  mylist =  ivector ( -1 , per ) ; 

  tot_per_col = ( per * Rows ) / Cols ;  
  relmax = per * ( tot_per_col + 2 ) ; /* was + 1 */
  relatives = ivector ( 1 , relmax + 10 ) ; 
  
  printf ( "Start:" ) ; fflush(stdout) ;

  for ( j = 1 ; j <= Cols ; j++ )  
    num_r[j] = per_col ; 
  if ( extras > 0 ) {
    fprintf ( stderr , "Adding %d extras because inexact division\n" , extras ) ; 
    for ( j = Cols ; j >= 1 && ( extras > 0 ); j-- )  
      if ( ranf() <=  (double) extras / (double) j ) 
	num_r[j] ++ ;
  }

  for ( i = num_rows_to_do , m = Rows - c0 ; i >= 1 ; i-- , m-- ) {
    di = 1.0 / (double) i ; 
    if ( !(i%1000) ) { printf ( "%d;" , i ) ; fflush(stdout) ; }

    /* compute cumulatives */
    for ( j = 1 ; j <= Cols ; j ++ ) 
      cum[j] = cum[j-1] + num_r[j] ; 

    hfails = 0 ; 
    do { /* generate a row */