fft.c

speech signal process tools

/*--------------------------------------------------------------+
|								|
|   This material contains proprietary software of Entropic	|
|   Processing, Inc.  Any reproduction, distribution, or	|
|   publication without the the prior written permission of	|
|   Entropic Processing, Inc. is strictly prohibited.  Any	|
|   public distribution of copies of this work authorized in	|
|   writing by Entropic Processing, Inc. must bear the notice	|
|								|
|          "Copyright 1986 Entropic Processing, Inc."		|
|								|
+---------------------------------------------------------------+
|								|
|  fft -- transform an SPS sampled-data file to an SPS		|
|         spectral file by discrete Fourier transform		|
|								|
|  Rodney Johnson, Entropic Processing, Inc.			|
|  Modified by John Shore to add -r option, and remove ndrec    |
|  dependence
|								|
+--------------------------------------------------------------*/
#ifdef SCCS
static char *sccs_id = "@(#)fft.c	1.1 6/8/87 ESI";  
#endif
char *Version = "1.1";
char *Date = "6/8/87";

#include <sps/sps.h>
#include <stdio.h>
#include <sps/sd.h>
#include <sps/spec.h>
/*
 * these defines for max quantities will be removed when 
 * sps.h gets them
 */
#define SPSHIBITL	(1L << (8 * (int)sizeof(long)) - 1)
#define SPSMAXLONG	(~SPSHIBITL)

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

#define Fprintf (void)fprintf
#undef ROUND
#define ROUND(x) (((x) > 0)? (long)(0.5+(x)) : -(long)(0.5-(x)))
#define SYNTAX USAGE("fft [-p range] [-r nframes] [-o order] sdfile specfile")
/*
 * global variables
 */
char *ProgName = "fft";
void get_range();
/*
 * external functions
 */
void    perror(), exit();
char    *strcpy();
double  log10(), pow();
char	*get_cmd_line();
void	lrange_switch();
/*
 * main program
 */
main(argc, argv)
    int  argc;
    char **argv;
{
    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 spec_data		/* record for spectral data */
            *spec_rec;
    extern int
            optind;
    extern char
            *optarg;
    int     ch;
    int     order = 10;		/* order of fft */
    char    *prange = NULL;	/* string for range specification (-p) */
    int	    p_flag = 0;		/* flag for range option */
    long    first, last, nan, length, i, j;
    long    nframes = 1;   /*number fixed length frames to process*/
    float   x[DFTMAX], y[DFTMAX];   /* arrays for data and fft of data */
    float   logscale;		/* scaling constant */
    int	    more = 1;		/* flag to indicate more sampled data*/
    long    position;		/* current position in SD file */
    /*
     * process command line options
     */
    while ((ch = getopt(argc, argv, "o:p:x:r:")) != EOF)
        switch (ch)    
        {
        case 'o':
            order = atoi(optarg);
	    break;
        case 'p': 
            prange = optarg;
	    p_flag++;
            break;
	case 'r':
	    nframes = atoi(optarg);
	    if (nframes == 0) nframes = SPSMAXLONG;
	    break;
        default: 
            SYNTAX;
            break;
        }
    /*
     * process file arguments
     */
    if (optind < argc) {
	iname = argv[optind++];
	if (strcmp (iname, "-") == 0)
	    iname = "<stdin>";
	else
	    TRYOPEN (argv[0], iname, "r", ifile);
    }
    else {
	Fprintf(stderr, "fft: no input SD file specified.\n");
	SYNTAX;
    }
    if (optind < argc) {
	oname = argv[optind++];
    	if (strcmp (oname, "-") == 0)
	    oname = "<stdout>";
	else
	    TRYOPEN (argv[0], oname, "w", ofile);
    }
    else {
	Fprintf(stderr, "fft: no output file specified.\n");
	SYNTAX;
    }

/*
 * read range from SPS common, if range option not used; 
 * the filename in common match that of the input SD file
 */
    if (!p_flag) read_params((char *)NULL, SC_CHECK_FILE, iname);
    get_range(&first, &last, prange, p_flag);
    symerr_exit();  /*exit if any of the parameters were missing*/

    if (first < 1)
	ERROR_EXIT("can't have negative starting point");
    if (last == first)
	ERROR_EXIT("range must specify more than one point");
    nan = last - first + 1;
    length = ROUND(pow(2.0, (double) order));
    if (length > DFTMAX)
    {
        Fprintf(stderr,
                "%s: order %d too high; %d points is max.\n",
                ProgName, order, DFTMAX);
        exit(1);
    }
    /*
     * only read as many points as length of transform
     */
    if (length < nan) nan = length;
    /*
     * read input SD file header
     */
    if ((ih = read_header(ifile)) == NULL)
	ERROR_EXIT("couldn't read input file header");
    if (ih->common.type != FT_SD)
	ERROR_EXIT("input file not a sampled-data file");
    /*
     * create, fill in, and write SPEC file header
     */
    oh = new_header(FT_SPEC);
    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;
    (void) add_comment(oh, get_cmd_line(argc, argv));
    oh->hd.spec->start = first;
    oh->hd.spec->nan = nan;
    oh->hd.spec->frmlen = nan;
    oh->hd.spec->win_type = NONE;
    oh->hd.spec->sf = ih->hd.sd->sf;
    oh->hd.spec->dcrem = 0.0;
    oh->hd.spec->voicing = NO;
    oh->hd.spec->freq_format = SYM_EDGE;
    oh->hd.spec->spec_type = ST_DB;
    oh->hd.spec->contin = YES;
    oh->hd.spec->num_freqs = length/2 + 1;
    oh->hd.spec->pre_emp = NULL;
    oh->hd.spec->frame_meth = FM_FIXED;
    *add_genhd_l("fft_length", NULL, 1, oh) = length;
    write_header(oh, ofile);
    /*
     * allocate spectral record and move to starting position
     */
    spec_rec = allo_spec_rec(oh); 
    if (first > 1) skiprec(ifile, first - 1, size_rec(ih));
    /*
     * main loop
     */
    logscale = 10 * log10(1.0);
    position = first;
    for (i = 0; i < nframes && more; i++) {
	/*
	 * initialize data arrays
	 */
	for (j = 0; j < length; j++) x[j] =  y[j] = 0.0;
	/*
	 * get sampled data and perform fft
	 */
    	if (get_sd_recf(x, (int) nan, ih, ifile) == EOF) {
	    more = 0;
	    Fprintf(stderr, "%s: WARNING, EOF reached in input file\n", ProgName);
	}
	get_fft(x, y, order);
	/*
	* fill in and write out spectral record
	*/
	spec_rec->tag = first;
	for (j = 0; j < oh->hd.spec->num_freqs; j++)
	    spec_rec->re_spec_val[j] =
                logscale + 10 * log10(x[j]*x[j] + y[j]*y[j]);
        spec_rec->tag = position;
	position += nan;
	put_spec_rec(spec_rec, oh, ofile);
    }
/*
 * put info in ESPS common
 */
    if (strcmp(oname, "<stdout>") != 0) {
	(void) putsym_s("filename", iname);
	(void) putsym_s("prog", ProgName);
	(void) putsym_i("start", (int) first);
	(void) putsym_i("nan", (int) nan);
    }    
    exit(0);
}



void
get_range(srec, erec, rng, rflag)
long *srec;			/* starting record */
long *erec;			/* end record */
char *rng;			/* range string from range option */
int rflag;			/* flag for whether range option used */
{
    *srec = 1;
    *erec = SPSMAXLONG;
    if (rflag) 
	lrange_switch (rng, srec, erec, 1);	
    else {
	if(symtype("start") != ST_UNDEF) *srec = getsym_i("start");
	if(symtype("nan") != ST_UNDEF) *erec = *srec + getsym_i("nan") - 1; 
    }
}