fft.c

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

/*
 * This material contains proprietary software of Entropic Speech, Inc.
 * Any reproduction, distribution, or publication without the prior
 * written permission of Entropic Speech, Inc. is strictly prohibited.
 * Any public distribution of copies of this work authorized in writing by
 * Entropic Speech, Inc. must bear the notice
 *
 *   "Copyright (c) 1988, 1990 Entropic Speech, Inc. All rights reserved."
 *
 * Program: fft.c (FEA_SPEC version)
 *
 * This program computes FFTs from SD files to FEA_SPEC files
 *
 *  Rodney Johnson, Entropic Speech, Inc.		
 *  Modified by John Shore to add various options, and remove ndrec 
 *  dependence;
 *  Modified for ESPS 3.0 by John Shore
 *  Revised by Shore for overlapping frames and windowing.
 *  Converted from SPEC to FEA_SPEC output by Rod Johnson.
 *  Converted from SD to FEA_SD input by David Burton
 */
#ifndef lint
    static char *sccs_id = "@(#)fft.c	3.30	7/8/96	 ESI";  
#endif
char *Version = "3.30";
char *Date = "7/8/96";

#include <stdio.h>
#include <esps/esps.h>
#include <esps/unix.h>
#include <esps/limits.h>
#include <esps/feaspec.h>
#include <esps/fea.h>
#include <esps/feasd.h>
#include <esps/window.h>

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

#define Fprintf (void)fprintf

#define SYNTAX \
USAGE("fft [-c] [-l frame_len] [-o order] [-{pr} range] [-w window_type]\n [-x debug_level] [-z] [-P param] [-S step] [-O freq. offset] sdfile specfile")

#define WT_PREFIX "WT_"

/*
 * global variables
 */

char	*ProgName = "fft";
int	debug_level = 0;		/* debug level */
void	get_range();
void	pr_farray();

/*  
 * external functions
 */
double  log10(), pow(), ceil(), sqrt();
char	*get_cmd_line();
void	lrange_switch();
char    *eopen();
void	get_rfft();
void    get_fft();
long	atol();
#if !defined(DEC_ALPHA) && !defined(HP700)
char    *calloc();
char    *malloc();
#endif
#ifndef DEC_ALPHA
char	*strcat();
#endif
/*
 * main program
 */
main(argc, argv)
    int	    argc;
    char    **argv;
{
    char    *param_file = NULL;	/* parameter file name */
    FILE    *ifile = stdin,	/* input and output file streams */
            *ofile = stdout;
    struct header  		/* input and output file headers */
            *ih, *oh;
    char    *iname, 		/* input and output file names */
	    *oname;
	    
    struct feaspec *spec_rec;   /* record for spectral data */

    struct feasd *sd_rec;       /* record for input data */

    extern int
            optind;
    extern char
            *optarg;
    int     ch;
    int     order;		/* order of fft */
    char    *prange = NULL;	/* string for range specification (-p) */
    int	    zflag = 0;		/* flag for -z (silent) option */
    int     oflag = 0;		/* flag for order option */
    int	    pflag = 0;		/* flag for range option */
    int	    cflag = 0;		/* flag for complex option */
    int	    wflag = 0;		/* flag for window option */
    int     lflag = 0;		/* flag for -l frame_len option */
    int	    Sflag = 0;		/* flag for -S option (step size) */
    char    *window_type;	/* name of type of window to apply to data */
    char    *pref_w_type;	/* window type name with added prefix */
    long    i, j;
    short   freq_format;
    short   spec_type;
    long    num_freqs;
    long    frame_len = -1;	/* number of sampled data points */
    long    tr_length;		/* number of points in transform */
    long    hf;			/* half number of frequencies */
    float   sf;
    int	    actsize;		/* actual number points read */
    int	    data_length;	/* actual number points transformed */
    long    nframes = 1;	/* number fixed length frames to process */
    float   fpow;		/* power at one frequency */
    float   tot_power;		/* total power */
    float   *x, *wx, *y, *wy;	/* arrays for data, windowed data, and fft */
    float   scale;		/* scaling constant */
    int	    more = 1;		/* flag to indicate more sampled data */
    long    position;		/* current position in SD file */
    long    step;		/* step size for shifting frames */
    long    start = 1;
    long    nan = 0; 
    int	    win;		/* window type code */
    extern char
	    *window_types[];	/* standard window type names */
    int	    init_err;		/* return code from init_sd_rec */
    int	    complex;		/* NO or YES for data type */
    int	    in_type;		/* input data type */
    int	    num_channels;	/* number of input data channels */
    float   *fdata;		/* pointer to data in fea_sd record */
    float_cplx
	    *cdata;		/* pointer to data in fea_sd record */
     double	offset;		/* frequency offset (-O option) */
     int	Oflag = 0;	/* offset was specified via -O */
     float	*freqs = NULL;	/* vector of freqs for freq. offset */
     double	freq;

    /*
     * process command line options
     */
    while ((ch = getopt(argc, argv, "cl:o:p:r:w:x:O:P:S:z")) != EOF)
        switch (ch)  {
	case 'c':
	    cflag++;
	    break;
	case 'l':
	    frame_len = atoi(optarg);
	    lflag++;
	    break;
        case 'o':
            order = atoi(optarg);
	    oflag++;
	    break;
        case 'p': 
        case 'r': 
            prange = optarg;
	    pflag++;
            break;
	case 'w':
	    window_type = optarg;
	    wflag++;
	    break;
	case 'x':
	    debug_level = atoi(optarg);
	    break;
 	case 'O':
 	    offset = atof(optarg);
 	    Oflag++;
 	    break;
	case 'z':
	    zflag++;
	    break;
	case 'P':
	    param_file = optarg;
	    break;
	case 'S':
	    step = atol(optarg);
	    Sflag++;
	    break;
        default: 
            SYNTAX;
            break;
        }

/*
 * check number of range specifications
 */
    if(pflag > 1){
	Fprintf(stderr,
		"fft: range should only be specified once - exiting.\n");
	exit(1);
    }

/*
 * process file arguments
 */
    if (optind < argc)
	iname = eopen(ProgName,
		      argv[optind++], "r", FT_FEA, FEA_SD, &ih, &ifile);
    else
    {
	Fprintf(stderr, "fft: no input sampled data file specified.\n");
	SYNTAX;
    }
    if (debug_level)
	Fprintf(stderr, "Input file is %s\n", iname);
/* 
 * check if input file is multichannel
 */
    if ((num_channels = 
	 get_fea_siz("samples", ih,(short *) NULL, (long **) NULL)) != 1)
    {
        Fprintf(stderr, "fft: Multichannel data not supported yet\n");
	exit(1);
    }

/*
 * Check if input data is complex
 */

    in_type = get_fea_type("samples", ih);
    switch (in_type)
    {
    case DOUBLE:
    case FLOAT:
    case LONG:
    case SHORT:
    case BYTE:
        complex = NO;
        break;
    case DOUBLE_CPLX:
    case FLOAT_CPLX:
    case LONG_CPLX:
    case SHORT_CPLX:
    case BYTE_CPLX:
        complex = YES;
        break;
    default:
        Fprintf(stderr, "%s: bad type code (%d) in input file header.\n",
                ProgName, in_type);
        exit(1);
        break;
    }
/*
 * setup output fea file
 */
    if (optind < argc)
	oname = eopen(ProgName, argv[optind++], "w", NONE, NONE, &oh, &ofile);
    else {
	Fprintf(stderr, "%s: no output FEA_SPEC file specified.\n", ProgName);
	SYNTAX;
    }
    if (debug_level)
	Fprintf(stderr, "Output file is %s\n", oname);
/*
 * process parameters
 */
    if (strcmp(iname, "<stdin>") != 0)
	(void) read_params(param_file, SC_CHECK_FILE, iname);
    else
        (void) read_params(param_file, SC_NOCOMMON, iname);

    if (!oflag)
	order = 
	    (symtype("order") != ST_UNDEF)
	    ? getsym_i("order")
	    : 10;
    
    if ((order >= 20) && !zflag)
	Fprintf(stderr,
		"%s: WARNING - HEY!!!, order = %ld makes for a huge transform\n",
		ProgName, order);

    tr_length = LROUND(pow(2.0, (double) order));

    get_range(&start, &nan, &frame_len, &step, &nframes, prange, 
	      pflag, lflag, Sflag, tr_length, ih, zflag);

    if (!wflag)
	window_type =
	    (symtype("window_type") != ST_UNDEF)
	    ? getsym_s("window_type")
	    : "RECT";

    pref_w_type = 
	malloc((unsigned)(strlen(WT_PREFIX) + strlen(window_type) + 1));
    spsassert(pref_w_type, "can't allocate space for window type name");
    (void) strcpy(pref_w_type, WT_PREFIX);
    (void) strcat(pref_w_type, window_type);

    win = lin_search(window_types, pref_w_type);

    spsassert(win > -1, "window type not recognized");

    if (debug_level)
	Fprintf(stderr, "%s: window_type = %s, win = %d\n",
		ProgName, window_type, win);

    if (debug_level)
	Fprintf(stderr,
		"%s: start = %ld, nan = %ld, order = %d, tr_length = %d\n", 
		ProgName, start, nan, order, tr_length);
    if (debug_level) 
	Fprintf(stderr, "%s: frame size = %ld, step = %ld, nframes = %ld\n",
		ProgName, frame_len, step, nframes);

/*
 * Allocate data arrays
 */
    x = (float *) calloc((unsigned) MAX(frame_len, tr_length), sizeof(float));
    y = (float *) calloc((unsigned) MAX(frame_len, tr_length), sizeof(float));
    spsassert(x != NULL && y != NULL, 
	      "fft: couldn't allocate enough memory");
    wx=(float *) calloc((unsigned) MAX(frame_len, tr_length), sizeof(float)); 
    wy=(float *) calloc((unsigned) MAX(frame_len, tr_length), sizeof(float)); 
    spsassert(wx != NULL && wy != NULL, 
	      "fft: couldn't allocate enough memory");

    if (start < 1)
	ERROR_EXIT("can't have negative starting point");

    /*
     * only read as many points as length of transform
     */
    if ((tr_length < frame_len) && !zflag) {
	Fprintf(stderr, 
	    "fft: WARNING - transform length %ld less than frame length %ld\n",
	    tr_length, frame_len);
	Fprintf(stderr,"\t\t(framing determined by step size)\n");
    }	
    if ((tr_length > frame_len) && !zflag) {
	Fprintf(stderr, 
		"fft: WARNING - transform length (%ld) greater than frame length (%ld)\n",
		tr_length, frame_len);
	Fprintf(stderr,
		"\t %ld zeros padded for each frame\n", tr_length-frame_len);
    }
    /*
     * create, fill in, and write output file header
     */
    oh = new_header(FT_FEA);
    add_source_file(oh, iname, ih);
    oh->common.tag = YES;
    (void) strcpy(oh->common.prog, ProgName);
    (void) strcpy(oh->common.vers, Version);
    (void) strcpy(oh->common.progdate, Date);
    oh->variable.refer = iname;
    add_comment(oh, get_cmd_line(argc, argv));

    update_waves_gen(ih, oh, (float) (start + frame_len/2), (float)step);

    sf = *get_genhd_d("record_freq", ih);

      if (cflag)
	{
	  freq_format = SPFMT_ASYM_EDGE;
	  spec_type = SPTYP_CPLX;
	  num_freqs = tr_length + 1;
	  hf = tr_length/2;
	}
      else if(complex) {
	if (Oflag) {
	  freq_format = SPFMT_ARB_FIXED;
	  spec_type = SPTYP_DB;
	  num_freqs = tr_length + 1;
	  hf = tr_length/2;
	  freqs = (float *) malloc(num_freqs * sizeof(float));
	  /* DEBUG */
	  if (debug_level) fprintf(stderr, "offset %g, sf %g\n", offset, sf);
	  for (freq = offset - sf/2.0, i = 0; i < num_freqs; i++, freq += sf/tr_length)
	    freqs[i] = freq;
	} else {
	  freq_format = SPFMT_ASYM_EDGE;
	  spec_type = SPTYP_DB;
	  num_freqs = tr_length + 1;
	  hf = tr_length/2;
	}