sfconvert.c

来自「speech signal process tools」· C语言 代码 · 共 1,201 行 · 第 1/3 页

/*
 * This material contains unpublished, proprietary software of 
 * Entropic Research Laboratory, Inc. Any reproduction, distribution, 
 * or publication of this work must be authorized in writing by Entropic 
 * Research Laboratory, Inc., and must bear the notice: 
 *
 *    "Copyright (c) 1990-1997 Entropic Research Laboratory, Inc. 
 *                   All rights reserved"
 *
 * The copyright notice above does not evidence any actual or intended 
 * publication of this source code.     
 *
 * Written by:  Derek Lin
 * Checked by:
 * Revised by:
 *
 * Brief description:
 *        This Program converts sampling frequency of FEA_SD data.  Sfconvert
 *        consists of three main modules.  gen_filter() designs the lowpass
 *        filter for interpolation.  gen_index() creates indexing array for
 *        both filter and sampled data.  convert() only computes points
 *        of interests for output by using the arrays created by last
 *        two modules.
 *
 */

static char *sccs_id = "@(#)sfconvert.c	1.14	24 Mar 1997	ERL";

#include <math.h>
#include <esps/esps.h>
#include <esps/constants.h>
#include <esps/fea.h>
#include <esps/feasd.h>
#include <esps/feafilt.h>
#include <esps/window.h>

#define DEBUG 0
#define BUF_SIZ 1000
#define MAX_NO_CHANNELS 512

#define SYNTAX USAGE("sfconvert [-P param_file] [-s new_sample_freq] [-r range][-v deviation] [-c corner_freq] [[-R rej_db][-t trans_band]] [-w sfwin_type -l sfwin_len] [-e channels] [-o output_type] [-f] [-F filt_file] [-d] [-x debug_level] [intput.SD] output.SD")

#define ERROR_EXIT(text) {(void) fprintf(stderr, "%s: %s - exiting\n", \
					 argv[0], text); exit(1);}

/*
 * Global declaration
 */

long *grange_switch();
long get_fea_siz();
short get_fea_type();
int debug_level = 0;
char *arr_alloc();

union data                /* pointer for diffrent types of input &   */
{                                /* output data array.                      */
  float *fsr;                    /*   f/d: float or double                  */
  float_cplx *fsc;               /*   s/m: single or multi channels         */
  float **fmr;                   /*   r/c: real or complex                  */
  float_cplx **fmc;
  double *dsr;
  double_cplx *dsc;
  double **dmr;
  double_cplx **dmc;
};

union filt                /* pointer to filter data, float or double */
{
  float *f;
  double *d;
};

/*
 * main program
 */

main(argc, argv)
     int argc;
     char **argv;
{
  char *Version = "1.0";
  char *Date = "9/10/92";
  char *param_file = NULL;           /* parameter file name */
  FILE *isdfile = stdin,             /* input and output file streams */
       *osdfile = stdout;
  char *iname = NULL,                /* input and output file names */
       *oname = NULL,
       *fname = NULL;
  struct header *ihd,                /* input and output file headers */
                *ohd;
  struct feasd *isd_rec,             /* record for input and output data */
               *osd_rec;

  extern int optind;
  extern char *optarg;
  extern char *window_types[];       /* window type code words, window.h */
  extern char *type_codes[];         /* data type code words. */
  int ch;
  
  int gen_index();                   /* prepares indexes for filtering*/
  long **kernel,                     /* filter kernels */
       *jumpsd;                      /* # of sample to jump for each conv. */
  long dimk[2];                      /* kernel dimension for arr_alloc argum.*/
  long ibuffer_len,                  /* input buffer length */
       obuffer_len,                  /* output buffer length */
       patch = 0;                    /* for some patch work for get_feasd_orecs
					in returning number of data points
					read at the end of file */
  long step1,                        /* 1st step argument for get_feasd_orecs*/
       stepn;                        /* subsequent step in get_feasd_orecs */
  long step;
  long actsize;                      /* # of valid samples read in on buffer */

  void gen_filter();                 /* lowpass design by Kaiser window */
  void write_filter_out();
  short win_type = WT_KAISER;        /* window type */
  float trans = 200;                 /* transistion band in the low pass */
  float db = 60;                     /* number of db down for alias */
  float len_s = 0.5;                 /* filter length in seconds */
  long len;                          /* filter length in samples*/
  union filt filter;
  double cor_freq = 0;               /* corner frequency, defualt is Nyquist*/
  double omega_c;                    /* corner frequency in radian */

  void convert();                    /* core for convolution computation */
  static union data 
    obuffer,                         /* output sd buffer */
    ibuffer;                         /* input sd buffer */
  long up = 1;                       /* upsampling factor */
  long down = 1;                     /* downsampling factor */
  
  long no_samples = 0;               /* number of samples in output file */
  double *istart_time,               /* input and output start_time */
         *ostart_time;            
  double *max_val;                   /* maximum value of output data file */
  double samp_rate;                  /* input sampling rate */
  double new_samp_rate = 0;          /* new sampling rate */
  void find_ratio();                 /* finds conversion factors up, down
					from 2 sampling freqs.      */
  float dev = 0;                     /* maximum tolerable deviation from      
					new sampling rate in percentage */
  int sflag = 0;                     /* flag for new_samp_rate */
  int rflag = 0;                     /* flag for selecting a range */
  int cflag = 0;                     /* flag for corner frequency */
  int dflag = 0;                     /* flag for double precision compute */
  int oflag = 0;                     /* flag for output data type */
  int vflag = 0;                     /* flag of maximum allowable deviation 
					from desired sampling frequency */
  int Rflag = 0;                     /* flag of rejection ratio in db from 
					pass to stop band */
  int tflag = 0;                     /* flag of transition bandwith in Hertz */
  int eflag = 0;                     /* flag of channels selection */
  int wflag = 0;                     /* flag of window selection  */
  int lflag = 0;                     /* flag of window length */
  int Kaiserflag = 1;                /* flag for Kaiser window in params file*/
  int fflag = 0;                     /* flag for saving filter coeff */
  long *channels = NULL;             /* array for multichannels */
  void get_range();                  /* get range function */
  char *range=NULL;                  /* string from -r option */
  long s_rec, e_rec;                 /* start and ending position of range*/
  long skip_rec = 0;                 /* number of records to skip */
  long no_samp_rf_left;              /* number of output samples left
					remaining to be written if range
					is specified */
  long num_channels = 0;             /* number of channels */
  int common_file = 0;               /* common processing flag */
  int input_type, output_type;       /* input and output data types */
  int code = 0;                      /* code word for data type and structure
				        default is 'fsr': float,single channels
					and real */
  char *sym;                         /* array to hold data from paramter file*/
  int i, type, size, check = 1;      /* utility variables */

/*
 * process command line options
 */

  while((ch = getopt(argc, argv, "P:s:r:v:c:R:t:w:l:e:o:F:x:df")) != EOF)
    switch (ch) {
    case 'd':
      dflag++;                             
      break;
    case 'f':
      fflag++;
      break;
    case 'F':
      fname = optarg;
      break;
    case 'P':
      param_file = optarg;
      break;
    case 's':
      new_samp_rate = atof(optarg);
      sflag++;
      break;
    case 'r':
      range = optarg;
      rflag++;
      break;
    case 'v':
      dev = atof(optarg) / 100.0;
      vflag++;
      break;
    case 'c':
      cor_freq = atof(optarg);
      cflag++;
      break;
    case 'R':
      db = atof(optarg);
      Rflag++;
      break;
    case 't':
      trans = atof(optarg);
      tflag++;
      break;
    case 'w':
      win_type = win_type_from_name(optarg);
      wflag++;
      break;
    case 'l':
      len_s = atof(optarg);
      lflag++;
      break;
    case 'e':
      channels = grange_switch(optarg, &num_channels);
      eflag++;
      break;
    case 'o':
      output_type = lin_search(type_codes, optarg);
      oflag++;
      break;
    case 'x':
      debug_level = atoi(optarg);
      break;
    default:
      SYNTAX;
      break;
    }

/*
 * Determine and open input/output file
 */

  if (argc - optind < 2) {
    Fprintf(stderr, "%s: no output file specified.\n",argv[0]);
    SYNTAX;
  }
  iname = eopen(argv[0],argv[argc-2],"r", FT_FEA, FEA_SD, &ihd, &isdfile);
  oname = eopen(argv[0], argv[argc - 1], "w", NONE, NONE, &ohd, &osdfile);

  if (debug_level){
    Fprintf(stderr, "%s: output file is %s\n", argv[0], oname);
    Fprintf(stderr, "%s: input file is %s\n", argv[0], iname);
  }

/*
 * Get parameter values
 */

  if (oflag == 1 && output_type == -1)
    ERROR_EXIT(" Unrecognized output data type with -o option.");
  if (wflag == 1 && win_type == WT_NONE)
    ERROR_EXIT(" Unrecognized window type with -w option.");

  (void) read_params(param_file, SC_CHECK_FILE, iname);

  if (!sflag && symtype("new_sample_freq") != ST_UNDEF)
    new_samp_rate = getsym_i("new_sample_freq");
  get_range( &s_rec, &e_rec, range, rflag, 0, ihd);
  if (symtype("start") != ST_UNDEF || symtype("nan") != ST_UNDEF)
    rflag = 1;
  if (!vflag && symtype("deviation") != ST_UNDEF)
    dev = getsym_d("deviation") / 100.0;
  if (!cflag && symtype("corner_freq") != ST_UNDEF)
    cor_freq = getsym_d("corner_freq");
  if (!Rflag && symtype("rej_db") != ST_UNDEF)
    db = getsym_d("rej_db");
  if (!tflag && symtype("trans_band") != ST_UNDEF)
    trans = getsym_d("trans_band");
  if (!wflag && symtype("sfwin_type") != ST_UNDEF)
  {
    sym = getsym_s("sfwin_type");
    win_type = win_type_from_name(sym);
    if (win_type == WT_NONE)
	ERROR_EXIT(" Unrecognized window type in param file.");
  }
  if( !lflag && symtype("sfwin_len") != ST_UNDEF)
    len_s = getsym_d("sfwin_len");
  if( !eflag && symtype("channels") != ST_UNDEF)
  {
    channels = grange_switch(getsym_s("channels"), &num_channels);
  }
  if( !oflag && symtype("output_type") != ST_UNDEF)
  {
    sym = getsym_s("output_type");
    output_type = lin_search(type_codes, sym);
  }
  if( !dflag)
    if( symtype("dflag") != ST_UNDEF)
      if ( getsym_i("dflag") == 1) dflag = 1;
  if ( symtype("Kaiserflag") != ST_UNDEF )
    Kaiserflag = getsym_i("Kaiserflag");

/*
 *  check for inconsistency 
 */

  /* up to now, all values of -l/-w/-R/-t are all defined either by command
     line, paramter files, or defaults.  To resolve conflicts of window type,
     win_type == non Kaiser only if -w exists || (sfwin_type exists && 
     Kaiserflag == 0 ).  All other case, win_type = kaiser */

  if( ! (wflag || symtype("sfwin_type") != ST_UNDEF && Kaiserflag == 0 ) )
    win_type = WT_KAISER;

  if (win_type == WT_ARB)
    ERROR_EXIT("Sorry, window type ARB not supported");
    
  if (win_type != WT_KAISER && len_s <= 0) 
    ERROR_EXIT("sfwin_len must be > zero.");
  if ( new_samp_rate < 0 ) ERROR_EXIT("new_sample_freq must >= zero.");
  if ( dev < 0 ) ERROR_EXIT("deviation must be >= zero.");
  if ( cflag && cor_freq < 0 ) ERROR_EXIT("corner_freq must be >= zero.");
  if ( db <= 0 ) ERROR_EXIT("rej_db must be > zero.");
  if ( trans <= 0 ) ERROR_EXIT("trans_band must be > zero.");

/*
 *   determine real/complex, single/multichannel, or float/double 
 *   input data kind for computation.  code 0 corresponds to float
 *   single channel, real input data
 */

  if (( input_type = get_fea_type( "samples", ihd )) == UNDEF )
    ERROR_EXIT(" Input file has no field named samples!");
  if ( input_type == FLOAT_CPLX || input_type == DOUBLE_CPLX ||
      input_type == LONG_CPLX || input_type == SHORT_CPLX ||
      input_type == BYTE_CPLX ) code += 1;
  if ( input_type == DOUBLE || input_type == DOUBLE_CPLX || dflag == 1) 
  {
    code +=100;
    dflag = 1;
  }
  if ( (size = get_fea_siz("samples", ihd, (short *)NULL,(long **)NULL)) > 
      MAX_NO_CHANNELS && num_channels > MAX_NO_CHANNELS )
    ERROR_EXIT(" Input file exceeds 512 maximum allowable channels");
  if ( size > 1 ) code +=10;
  if ( num_channels > size )
    ERROR_EXIT(" Input file has less channels than specified");
  if ( num_channels >= 1)
    if ( channels[num_channels-1] >= size )
      ERROR_EXIT("  Non-existent channel specified");
  if ( num_channels == 0)
  {
    num_channels = size;
    channels = (long *) malloc( (unsigned)num_channels * sizeof(int));
    spsassert(channels,"malloc failed in main");
    for( i= 0; i<num_channels; i++) channels[i] = i;
  }

  if(debug_level)
  {
    if (1 <= code%10) Fprintf(stderr, "\nComplex input data.\n");
    if (1 > code%10 ) Fprintf(stderr, "\nReal input data.\n");
    if (10 > code%100 ) Fprintf(stderr, "Single channel input data.\n");
    if ( code < 100  ) Fprintf(stderr, "Floating point computation.\n\n");
    if ( code >=100  ) Fprintf(stderr, "Double precision computation.\n\n");
    if (10 <= code%100) 
    {
      Fprintf(stderr, "Multichannel input data.\n");
      Fprintf(stderr, "Number of output channels: %i\n",num_channels);
    }
  }
  
/*
 *   find up/down ratio, change new sampling rate accordingly within
 *   deviation level to get the smallest conversion factor 
 */

  if( (type = genhd_type("record_freq",&size,ihd)) == HD_UNDEF)
    ERROR_EXIT("record_freq header item is missing from input file");
  samp_rate = *get_genhd_d("record_freq", ihd);
  if( new_samp_rate ==0 ) new_samp_rate = samp_rate;
  find_ratio( &new_samp_rate, samp_rate, dev, &up, &down);
  len = (long) (len_s * samp_rate * MAX(up,down));
  if (win_type != WT_KAISER && len > 5000 && debug_level)
    Fprintf(stderr,"%s: A window of %d data points is used.\n", argv[0], len);
  if (len%2 == 0) len++;              /* make sure odd window length */
      

/*
 *   Compute low pass filter using Kaiser or other window design method
 */

  if( cor_freq > (samp_rate) / 2){
    Fprintf(stderr, "%s: Corner frequency violates the Nyquist rate.\n",
	    argv[0]);
    exit(1);
  }
  omega_c =( cor_freq !=0 )? 2*PI*cor_freq/(up*samp_rate) :  PI / MAX(up,down);
  (void) gen_filter( dflag, win_type, up, 2*PI*trans/(up*samp_rate),