code4.c

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

/* 
   code4.c                                   (c) DJCM 94 07 04
                                                      94 10 26
						      94 10 31
						      95 01 02
						      95 02 23 (from mnc2 to mnc3)
						      95 04 01 continued work on mnc3
						      95 05 01 to mnc4
						      95 05 02 to code4

   - stripped out of mnc4.c to be a special purpose code generator and writer; 
   makes two sparse invertible matrices and writes out the alist matrix, 
   to save everyone making the same inverse and also to track down the 
   strange bug.

   usage: 

   code4 -n N -seed seed -t1 5 -t2 5 -o outfile (Default MNC4/N.t1.t2.seed)

   This program has been used to generate the following default codes: 

   code4 -n 1000 -t1 5 -t2 5 -v 0
   nice +19 code4 -n 2000 -t1 5 -t2 5 -v 0 
   nice +19    code4 -n 4000 -t1 5 -t2 5 -v 0 
   nice +19    code4 -n 2000 -t1 3 -t2 3 -v 0 
   nice +19    code4 -n 4000 -t1 3 -t2 3 -v 0 
   nice +19    code4 -n 2000 -t1 4 -t2 4 -v 0 
   nice +19    code4 -n 4000 -t1 4 -t2 4 -v 0 
   nice +19    code4 -n 2000 -t1 6 -t2 6 -v 0 
   nice +19    code4 -n 4000 -t1 6 -t2 6 -v 0 
   nice +19    code4 -n 1000 -t1 3 -t2 3 -v 0 
   nice +19    code4 -n 1000 -t1 4 -t2 4 -v 0
   nice +19    code4 -n 1000 -t1 6 -t2 6 -v 0
*/

#include "./ansi/r.h"
#include "./ansi/rand.h"
#include "./ansi/mynr.h"
#include "./ansi/cmatrix.h"

/* these three are irrelevant: */
#include "./ansi/macopt.h"
#include "./thing_sort.h"
#include "./fe.h"

#include "./mnc.h"
#include "./code4.h"

static void dc_defaults ( code_creation * ) ; 
static int    process_command ( int , char **  , code_creation *  ) ; 
static void   print_usage ( char ** , FILE * ,   code_creation *  );

static int make_sense ( code_creation * ) ; 
static int make_code ( code_creation * ,  mnc_code * , alist_matrix * ) ;
static int report_code ( code_creation * , mnc_code * , alist_matrix * ) ;
static int make_space ( code_creation * , mnc_code * , alist_matrix * ) ; 
static int check_alist ( alist_matrix *  , code_creation * ) ;
static int free_code ( code_creation * , mnc_code * ) ;

void   main ( int , char ** ) ;

/*
        MAIN
                     */
void main ( int argc, char *argv[] )
{
  FILE   *fp  ; 

  code_creation dc   ;
  mnc_code            code ; 
  alist_matrix a ; 

#ifdef _DEBUG_MALLOC_INC
  unsigned long histid1, histid2 , orig_size , current_size ; 
/*
  double *tmpdv ; 
  unsigned char *tmpuc ; 

  union dbmalloptarg  m;
  m.i = M_HANDLE_CORE | M_HANDLE_DUMP;
  dbmallopt(MALLOC_WARN,&m);

  m.i = M_HANDLE_ABORT;
  dbmallopt(MALLOC_FATAL,&m);

  m.i = 1;
  dbmallopt(MALLOC_CKCHAIN,&m);

  m.str = "log";
  dbmallopt(MALLOC_ERRFILE,&m); */
#endif

#ifdef _DEBUG_MALLOC_INC

#endif

  dc_defaults ( &dc ) ; 

  if ( process_command (argc, argv, &dc ) < 0 ) exit (0) ;
  if ( make_sense ( &dc ) < 0 ) exit (0) ;
  if ( make_space ( &dc , &code , &a ) < 0 ) exit (0) ;
/* here is the space for the alist */

  fprintf(stderr,"code4 MNC=%d, N=%d\n",
	    dc.MNC , dc.N ) ;
  fflush(stderr);

#ifdef _DEBUG_MALLOC_INC
  malloc_chain_check(1) ; 
  malloc_enter ("make_code") ;
  orig_size = malloc_inuse ( &histid1 ) ; 
#endif

  if ( make_code ( &dc , &code , &a ) < 0 ) { 
    report_code ( &dc , &code , &a  ) ; exit ( 0 ) ;     
  }
  if ( check_alist ( &a , &dc ) < 0 ) { 
    report_code ( &dc , &code , &a ) ; exit ( 0 ) ;
  }
  free_code ( &dc , &code ) ;

#ifdef _DEBUG_MALLOC_INC
  malloc_leave ("make_code") ;
  current_size = malloc_inuse ( &histid2 ) ; 
  if ( current_size != orig_size ) {
    fprintf ( stderr , "memory code?\n" ) ; 
    printf ( "mZ:%d\n" , malloc_chain_check(1) ) ; fflush(stdout) ;
    if (dc.verbose)     malloc_list ( 2 , histid1 , histid2 ) ; 
  }
#endif

  if ( dc.out ) {
    fp = fopen ( dc.outfile , "w" ) ; 
    if ( !fp ) {
      fprintf ( stderr , "can't open %s\n" , dc.outfile ) , 
      exit (0) ; 
    }
    write_alist ( fp , &a ) ; 
    fclose ( fp ) ; 
  }
  free_alist ( &a ) ; 
  if ( dc.verbose ) {
#ifdef _DEBUG_MALLOC_INC
    orig_size = malloc_inuse ( &histid1 ) ; 
#endif
    read_allocate_alist ( &a , dc.outfile ) ; 
    if ( dc.out ) {
      fp = fopen ( "tmpout" , "w" ) ; 
      if ( !fp ) {
	fprintf ( stderr , "can't open %s\n" , dc.outfile ) , 
	exit (0) ; 
      }
      write_alist ( fp , &a ) ; 
      fclose ( fp ) ; 
    }
    free_alist ( &a ) ; 
#ifdef _DEBUG_MALLOC_INC
    current_size = malloc_inuse ( &histid2 ) ; 
    if ( current_size != orig_size ) {
      fprintf ( stderr , "memory alist?\n" ) ; 
      printf ( "mZ:%d\n" , malloc_chain_check(1) ) ; fflush(stdout) ;
      if (dc.verbose)     malloc_list ( 2 , histid1 , histid2 ) ; 
    }
#endif
    printf ( "diff tmpout %s\n" , dc.outfile ) ; 
  }
#ifdef _DEBUG_MALLOC_INC
  if (dc.verbose) malloc_dump(1) ; 
#endif

}

static int make_sense ( code_creation *dc ) 
{ /* routine to correct silly control parameters */
/* 
   first the data creation 
                             */
  int status = 0 ; 
  char junk[100] ; 

  switch ( dc->MNC ) {
  case(0) :
    if ( dc->M == 0 ) {
/* if MNC == 0, M and N better had got set on command line */
      fprintf ( stderr , " No M specified \n" ) ;  status-- ; }
    break ; 
  case ( 1) : 
  case ( 2) : 
  case ( 3) : 
  case ( 4) : 
    dc->M = dc->N ; 
    break ; 
  default : 
    fprintf ( stderr , "invalid MNC = %d\n" , dc->MNC ) ; 
    break ; 
  }

  sprintf ( junk , "MNC%d/%d.%d.%ld" , 
	   dc->MNC , dc->N , dc->pC1.per_row ,/*  dc->pC2.per_row , */
	   dc->seed 
	   ) ; 
  if ( dc->out==1 )  
    sprintf ( dc->outfile , "%s" ,  junk ) ; 

  return status ; 
}

static int make_space ( code_creation *dc  , 
		       mnc_code *code , alist_matrix *a) {
/* see also mnc_free */
  int status = 0 ; 
  int tr , tc , M , N ; 

  malloc_enter ("make_space") ;
  if ( dc->verbose ) printf ( "Making space\n" ) ; 

/* 
   Now estimate size of alist matrix 

     MNC1: Number per row and column is almost fixed to dc.pC1.per_row
     MNC2: Per row is the max of C1 and C2. Per col is the sum plus a bit.
     MNC3: As MNC1, plus 1. 
     MNC4: converse of MNC2.
     MNC0: identity matrix followed by a load of sparse junk 
*/
  switch ( dc->MNC ) {
  case ( 1 ) : 
    N = dc->N ; M = dc->M ; 
    tr = dc->pC1.per_row + 2 ; 
    tc = tr + 10 ; break ;
  case ( 2 ) : 
    N = dc->N ; M = 2 * dc->M ; 
    tr = MAX ( dc->pC1.per_row , dc->pC2.per_row ) + 2 ; 
    tc = 2 * tr + 10 ; break ;
  case ( 3 ) : 
    N = 2 * dc->N ; M = dc->M ; 
    tr = dc->pC1.per_row + 3 ; 
    tc = tr + 10 ; break ;
  case ( 4 ) : 
    N = 2 * dc->N ; M = dc->M ; 
    tc = MAX ( dc->pC1.per_row , dc->pC2.per_row ) + 5 ; 
    tr = 2 * tc ; break ;
  default : 
  case ( 0 ) :
    N = dc->N ; M = dc->M ; 
    tr = dc->pC1.per_row + 1 ;
    tc = 1 + ( dc->M - dc->N ) * tr * 3 / ( 2 * dc->N )  ; /* 3/2 for safety */
    break ;
  }

  /* Set up memory , etc. */

  initialize_alist ( a , N , M , tc , tr ) ;

  malloc_leave ("make_space");
  return status ; 
}

static int make_code ( code_creation *dc , mnc_code *code , alist_matrix *a ) {
  int status = 0 ; 
  FILE   *fp  ; 
  char 	junk[100] ; 

  if ( dc->verbose > 1 ) printf ( "s.%6ld." , dc->seed ) ; 
  fflush (stdout ) ; 
  if ( dc->verbose ) printf ( "Making code\n" ) ; 
  ran_seed ( dc->seed ) ; 

  switch ( dc->MNC ) {
  case ( 2) :
    code->C1 = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 
    code->C1I = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 
    sparse_invertible_cmatrix ( code->C1 , code->C1I , dc->pC1.per_row , dc->N ) ;
  
    code->C2 = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 
    code->C2I = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 
    sparse_invertible_cmatrix ( code->C2I , code->C2 , dc->pC2.per_row , dc->N ) ;
  
    status += cmatrices_2_alist ( code->C1 , code->C2I , a , dc->N , 1 ) ;