impulse_rs.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."		|
|								|
+---------------------------------------------------------------+
|								|
|  impulse_resp -- impulse response of filters			|
|  								|
|  by Rodney Johnson, EPI					|
|  								|
|  Module:  impulse_resp.c					|
|  								|
|  This, the main module of impulse_resp, interprets the        |
|  command line, sets up and checks input and output files,	|
|  reads the input data, computes the output values with the	|
|  help of block_filter, and writes the output.			|
|  								|
|  Input is an SPS filter file.  Output is and SPS filter	|
|  file or sampled-data file.					|
|								|
+--------------------------------------------------------------*/

#ifndef lint
static char *sccs_id = "@(#)impulse_rs.c	1.1 12/17/86";
#endif

#include <stdio.h>
#include <sps/sps.h>
#include <sps/filt.h>
#include <sps/sd.h>

#define Version "1.1"
#define Date "12/17/86"
#define ProgName "impulse_resp"
#define SYNTAX USAGE("impulse_resp -n size [-x debug_level]\
[-f range][-i][-s] input.filt output")
#define FILETYPE {Fprintf(stderr, \
        "%s: internal error--confused about output file type.", \
        ProgName); exit(1);}

char    *strcpy(), *calloc();
void    perror(), exit();

char    *get_cmd_line();
void    lrange_switch();
int     block_filter();
void    set_sd_type();

void    pr_zfunc(), tab();

char    head_errtxt[100];

void
main(argc, argv)
    int     argc;
    char    **argv;
{
    struct header
            *ih, *oh;
    FILE    *ifile, *ofile;
    char    *iname, *oname;
    struct filt_data
            *i_filt_rec, *o_filt_rec;

    char    *cmd_line = get_cmd_line(argc, argv);
    extern  optind;
    extern char
            *optarg;
    int     c;
    int     size = -1;      /* invalid value; no default */
    int     debug_level = 0;
    char    *range = NULL;
    long    start_rec, end_rec, num_rec;
    short   invert = NO;
    int     out_type = FT_FILT;

    struct zfunc
            zf;
    static float
            zero[1] = {0.0},
            one[1] = {1.0};
    float   *impulse, *state, *imp_resp;
    int     st_size;
    int     n, i;

/* read command-line options */

    while ((c = getopt(argc, argv, "n:x:f:is")) != EOF)
    switch (c)
    {
    case 'n':
        size = atoi(optarg);
        break;
    case 'x':
        debug_level = atoi(optarg);
        break;
    case 'f':
        range = optarg;
        break;
    case 'i':
        invert = YES;
        break;
    case 's':
        out_type = FT_SD;
        break;
    default:
        SYNTAX;
        break;
    }

/* open files and read headers */

    if (optind >= argc) SYNTAX
    else
    {
        iname = argv[optind++];
        if (strcmp(iname, "-") == 0)
        {
            iname = "<stdin>";
            ifile = stdin;
        }
        else
            TRYOPEN(ProgName, iname, "r", ifile)
    }
    if ((ih = read_header(ifile)) == NULL)
    {
        Fprintf(stderr, head_errtxt);
        NOTSPS(ProgName, iname);
    }
    if (ih->common.type != FT_FILT)
    {
        Fprintf(stderr,
                "%s: input file %s is not an SPS filter file\n",
                ProgName, iname);
        exit(1);
    }

    if (optind >= argc) SYNTAX
    else
    {
        oname = argv[optind++];
        if (strcmp(oname, "-") == 0)
        {
            oname = "<stdout>";
            ofile = stdout;
        }
        else
            TRYOPEN(ProgName, oname, "w", ofile)
    }
    
    if (optind < argc) SYNTAX

    if (size < 0)
    {
        Fprintf(stderr,
            "%s: output length unspecified or negative\n",
            ProgName);
        exit(1);
    }

    start_rec = 1;
    end_rec = num_rec = ih->common.ndrec;
    lrange_switch(range, &start_rec, &end_rec, 0L);
    if (end_rec < start_rec
        || end_rec > num_rec || start_rec < 1)
    {
        Fprintf(stderr,
                "%s: bad range: %ld records in %s\n",
                ProgName, num_rec, iname);
        exit(1);
    }
    num_rec = end_rec - start_rec + 1;

/* Make output header */

    oh = new_header(out_type);
    Strcpy(oh->common.prog, ProgName);
    Strcpy(oh->common.vers, Version);
    Strcpy(oh->common.progdate, Date);
    add_source_file(oh, iname, ih);
    if (add_comment(oh, cmd_line) == 0)
        Fprintf(stderr,	"%s: WARNING - %s\n",
            ProgName,
            "not enough space for command line in comment");
    oh->variable.refer = ih->variable.refer;
    (void) add_genhd_l("start_rec", &start_rec, 1, oh);
        /* record number of first input record in output header */
    (void) add_genhd_s("invert", &invert, 1, oh);
        /* record -i option in output header */

    switch (out_type)
    {
    case FT_FILT:
        oh->common.ndrec = num_rec;

        oh->hd.filt->max_num = size;
        oh->hd.filt->max_den = 0;
        oh->hd.filt->nbands = 0;
        oh->hd.filt->npoints = 0;
        oh->hd.filt->g_size = 0;
        oh->hd.filt->nbits = 0;
        oh->hd.filt->type = NONE;
        oh->hd.filt->method = NONE;
        oh->hd.filt->bandedges = NULL;
        oh->hd.filt->points = NULL;
        oh->hd.filt->gains = NULL;
        oh->hd.filt->wts = NULL;
        break;
    case FT_SD:
        set_sd_type(oh, FLOAT);
        oh->common.ndrec = num_rec * size;
        oh->hd.sd->sf = 1.0;
        break;
    default:
        FILETYPE
        break;
    }

    write_header(oh, ofile);

/* set up input and output filter records */

    i_filt_rec = allo_filt_rec(ih);
    imp_resp = (float *) calloc((unsigned) size, sizeof (float));
    switch (out_type)
    {
    case FT_FILT:
        o_filt_rec = allo_filt_rec(oh);
        o_filt_rec->filt_func.nsiz = size;
        o_filt_rec->filt_func.dsiz = 0;
        o_filt_rec->filt_func.zeros = imp_resp;
        o_filt_rec->filt_func.poles = NULL;
        break;
    case FT_SD:
        break;
    default:
        FILETYPE
        break;
    }

/* other initializations */

    impulse = (float *) calloc((unsigned) size, sizeof (float));
    if (size > 0) impulse[0] = 1.0;
    for (i = 1; i < size; i++) impulse[i] = 0.0;

    st_size = ih->hd.filt->max_num;
    if (st_size < ih->hd.filt->max_den)
        st_size = ih->hd.filt->max_den;
    state = (float *) calloc((unsigned) st_size, sizeof (float));

    if (debug_level >= 1)
    {
        Fprintf(stderr,
            "%ld records. first: %ld, last: %ld.\n",
            num_rec, start_rec, end_rec);
        Fprintf(stderr,
            "size: %d, max_num: %d, max_den: %d, st_size: %d\n",
            size, ih->hd.filt->max_num,
            ih->hd.filt->max_den, st_size);
    }

/* skip to first input record to process */

    for (n = 1; n < start_rec; n++)
        if (get_filt_rec(i_filt_rec, ih, ifile) == EOF)
        {
            Fprintf(stderr,
                "%s: premature EOF encountered\n", ProgName);
            exit(1);
        }

/* main loop */

    for ( ; n <= end_rec; n++)
    {
        if (get_filt_rec(i_filt_rec, ih, ifile) == EOF)
        {
            Fprintf(stderr,
                "%s: premature EOF encountered\n", ProgName);
            exit(1);
        }
        if (invert)
        {
            zf.nsiz = i_filt_rec->filt_func.dsiz;
            zf.dsiz = i_filt_rec->filt_func.nsiz;
            if (zf.nsiz == 0) /* Null denominator means 1. */
            {
                zf.nsiz = 1;
                zf.zeros = one;
            }
            else
                zf.zeros = i_filt_rec->filt_func.poles;
            if (zf.dsiz == 0) /* Null numerator means 0. */
            {        /* Let block_filter handle this error. */
                zf.dsiz = 1;
                zf.poles = zero;
            }
            else
                zf.poles = i_filt_rec->filt_func.zeros;
        }
        else
            zf = i_filt_rec->filt_func;

        if (debug_level >= 1)
        {
            Fprintf(stderr, "record number %d.  tag = %ld.\n",
                n, i_filt_rec->tag);
            if (debug_level >= 2)
            pr_zfunc(0, "zfunc", &zf);
        }

        for (i = 0; i < st_size; i++)
            state[i] = 0.0;
        (void) block_filter(size,
            impulse, imp_resp, &zf, state);
        switch (out_type)
        {
        case FT_FILT:
            o_filt_rec->tag = i_filt_rec->tag;
            put_filt_rec(o_filt_rec, oh, ofile);
            break;
        case FT_SD:
            put_sd_recf(imp_resp, size, oh, ofile);
            break;
        default:
            FILETYPE
            break;
        }
    }

    exit(0);
    /* NOTREACHED */
}

/*
 * for debug printout of zfuncs
 * pr_zfunc and tab lifted from psps source directory
 * printf( ) changed to Fprintf(stderr, )
 */
void
pr_zfunc(lev,label,zf)
    int lev;
    char *label;
    struct zfunc *zf;
{
    int i;

    if (zf == NULL)
    {
        tab(lev);
	Fprintf(stderr, "No %s zfunc\n",label);
	return;
    }
    tab(lev);
    Fprintf(stderr,
        "%s: nsiz: %d, dsiz: %d, { ", label, zf->nsiz, zf->dsiz);
    for (i = 0; i < zf->nsiz; i++)
        Fprintf(stderr, "%.4g ", zf->zeros[i]);
    Fprintf(stderr, "} / { ");
    for (i = 0; i < zf->dsiz; i++)
        Fprintf(stderr, "%.4g ", zf->poles[i]);
    Fprintf(stderr, "}\n");
}

void
tab(indent)
{
    while (indent--) Fprintf(stderr, "  ");
}