iir_filt.c.sta

来自「speech signal process tools」· STA 代码 · 共 1,535 行 · 第 1/3 页

/*********************************************************
*
  This material contains proprietary software of Entropic Speech, Inc.   
  Any reproduction, distribution, or publication without the 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 Entropic Speech, Inc.; All rights reserved"
 				

 	
*
*  Module Name: iir_filt.c
*
*  Written By:   David Burton
*  Checked by: 
*
*  DESCRIPTION:  This program designs a recursive filter.
*		 Only Butterworth and Chebyshev1 filters 
*		 are supported at this time.
*		 
*
*
*********************************************************/

static char *sccs_id = "@(#)iir_filt.c	1.2 6/6/88 ESI";

char *Version = "1.2", *Date="6/6/88";

/*
 * System Includes
*/
# include <stdio.h>
# include <math.h>

/* 
 * ESPS Includes
*/
# include <esps/esps.h>
# include <esps/filt.h>
# include <esps/unix.h>
# include <esps/constants.h>

/* 
 * Defines
*/
# define SYNTAX	USAGE ("iir_filt [-P param_file][-x debug_level] filt_file");
# define SQRD(x) (x)*(x)
# define ARCSINH(x) log((x)+sqrt(SQRD(x)+1))

/* 
* ESPS Functions
*/
    char *get_cmd_line();
    double_cplx cdiv(), cmult(), cadd(), csub(), realmult(), csqrt();
    double modul();
    char *calloc();
    int pz_to_co();

/*
 * Globals
*/
    extern char *optarg;
    extern optind;


main (argc, argv)
int argc;
char *argv[];
    {
    
    FILE *fpout=NULL;		    /*output stream pointer*/
    struct header *oh=NULL;		    /*output file header pointer*/
    struct filt_data *frec=NULL;	    /*output data rescord structure*/
    char *param_file = "params";    /*parameter file name*/
    char *filt_file = NULL;	    /*output file name*/

    int debug_level = 0;
    short order;		    /*holds final filter order*/
    int nroots_num, nroots_den;	    /*# roots in final filter*/

    float *num_coeff=NULL, *den_coeff=NULL;   /*holds designed coefficients*/
    double *r_pole=NULL, *r_zero=NULL, *i_pole=NULL, *i_zero=NULL;
    float gain = 1.0;		    /*scaling gain*/
    float sf = 8000;		    /*sampling frequency*/
    char *filt_mthd = NULL;	    /*filter polynomial type*/
    int filter_poly = 0;	    /* integer representation of type*/
    char *filt_resp_type =  NULL;   /*filter response type*/
    int filter_response = 0;	     /* integer representation of type*/
    char *tmp = NULL;		    /* temporary character ptr*/

    int i, c;
    int filter_order = 0;	    /* filter polynomial order */
    float three_dBs[2];		    /* holds 3 dB values */
    float warped_3dBs[3];	    /* holds prewarped 3 dB values*/
    int num_of_3dBs;		    /* number of 3 dB values*/
    float band_edges[4];	    /* holds band edges for ELLIPTICAL*/
    float warped_edges[4];	    /* holds prewarped edges for ELLIPTICAL*/
    int num_of_band_edges;	    /*number of band edge freqs*/
    float pass_loss=0;		    /*max. pass band loss*/
    float stop_loss=0;		    /*min. stop band loss*/
    short num_num_co, num_den_co;   /* number of filter coeffs*/
    float gain_factor;		    /*values used to scale numerator coeffs*/

    float prewarp();
    void low_pass_proto();
    void poles_and_zeros();
    void elliptical_p_and_z();
    void freq_xfrm();
    void blt();
    void pz_to_filt_co();
    void set_gain();
    double sqrt(), pow(), tan(), cos(), atan2();
    double sin(), atof(), sinh(), cosh(), fabs();


/* Check the command line options. */

    while ((c = getopt (argc, argv, "P:x:")) != EOF)
	{
	switch (c)
	    {
	    case 'P':
		param_file = optarg;
		break;
	    case 'x':
		debug_level = atoi(optarg);
		break;
	    default:
		SYNTAX
	    }
	}


/* allocate space for filt_resp_type */

    filt_resp_type = (char *) calloc((unsigned)8, sizeof(char));


/* Get the output filter filename */

    if (optind < argc)
	{
	filt_file = argv [optind];
	if (strcmp (filt_file, "-") == 0) 
	    {
	    filt_file = "<stdout>";
	    fpout = stdout;
	    }
	else TRYOPEN("iir_filt", filt_file, "w", fpout);
	if (debug_level) 
	    Fprintf (stderr,"iir_filt: Output file is %s\n", filt_file);
	}
    else
	{
	Fprintf (stderr,"iir_filt: No output file specified.\n");
	SYNTAX
	exit (1);
	}

/* Read parameter file and get design parameters */

    if (debug_level)
	Fprintf (stderr,"iir_filt: Reading parameter file %s\n", param_file);

    if (read_params (param_file, SC_NOCOMMON, (char *)NULL) != 0)
	{
	Fprintf (stderr, 
    "iir_filt: read_params could not read the params file (%s) - exiting.\n",
		param_file);
        exit(1);
	}


	if(symtype("filt_order") != ST_UNDEF){
	    filter_order = getsym_i("filt_order");
	}
	else{
	    Fprintf(stderr,	
           "irr_filt: Filt_order not specified in params file - exiting.\n");
	    exit(1);
	}

	
	if (symtype("samp_freq") != ST_UNDEF)
	    sf = (float)getsym_d ("samp_freq");
	else{
	    Fprintf(stderr, 
    "iir_filt: Sampling frequency not specified in params file - exiting\n");
	    exit(1);
	}

	if (symtype("gain") != ST_UNDEF)
	    gain = (float)getsym_d ("gain");
	else{
	    Fprintf(stderr, 
    "iir_filt: Pass band gain not specified in params file - exiting\n");
	    exit(1);
	}

	if (symtype("filt_method") != ST_UNDEF)
	    filt_mthd = getsym_s ("filt_method");
	else{
	    Fprintf(stderr, 
    "iir_filt: Filt_method not specified in params file - exiting\n");
	    exit(1);
	}


	if (symtype("filt_type") != ST_UNDEF){
	    tmp = getsym_s ("filt_type");
	    /* add FILT_ prefix*/
	    (void)strcat(filt_resp_type, "FILT_");
	    (void)strcat(filt_resp_type, tmp);
	}
	else{
	    Fprintf(stderr, 
    "iir_filt: Filt_type not specified in params file - exiting\n");
	    exit(1);
	}

/* convert character strings to symbolic constants*/

    if((filter_response = lin_search(filt_type, filt_resp_type)) == -1){
       Fprintf(stderr, 
       "iir_filt: Invalid filt_type (%s) specified - exiting.\n", 
	tmp);
       exit(1);
    }

    if((filter_poly = lin_search(filt_method, filt_mthd)) == -1){
	Fprintf(stderr, 
	    "iir_filt: Invalid filt_method (%s) specified - exiting.\n",
	    filt_mthd);
	exit(1);
    }

	
/* Get filter design parameters - 
    params values depend on filter polynomial type */

    if(filter_poly == ELLIPTICAL || filter_poly == CHEBYSHEV2){
	Fprintf(stderr, 
  "iir_filt: %s filters not supported at this time - exiting.\n", filt_mthd);
	exit(1);
    }
    switch(filter_poly){
	case BUTTERWORTH:
	    if(symtype("three_dB_freqs") != ST_UNDEF){
		num_of_3dBs = getsym_fa("three_dB_freqs", three_dBs, (int)2);
	    }
	    else{
		Fprintf(stderr,
		    "iir_filt: three_db_freqs not specified - exiting.\n");
		exit(1);
	    }
	    break;

	case CHEBYSHEV1:
	    if(symtype("three_dB_freqs") != ST_UNDEF){
		num_of_3dBs = getsym_fa("three_dB_freqs", three_dBs, (int)2);
	    }
	    else{
		Fprintf(stderr,
		    "iir_filt: three_db_freqs not specified - exiting.\n");
		exit(1);
	    }

	    if(symtype("pass_band_loss") != ST_UNDEF)
		pass_loss = (float)getsym_d("pass_band_loss");
	    else{
		Fprintf(stderr,
		    "iir_filt: Pass_band_loss not specified - exiting.\n");
		exit(1);
	    }

	    if( pass_loss <= 0.){
		Fprintf(stderr, 
		    "iir_filt: Pass band loss must be > 0 - exiting.\n");
		exit(1);
	    }

	    break;

	case CHEBYSHEV2:
	    if(symtype("three_dB_freqs") != ST_UNDEF){
		num_of_3dBs = getsym_fa("three_dB_freqs", three_dBs, (int)2);
	    }
	    else{
		Fprintf(stderr,
		    "iir_filt: three_db_freqs not specified - exiting.\n");
		exit(1);
	    }


	    if(symtype("stop_band_loss") != ST_UNDEF)
		stop_loss = (float)getsym_d("stop_band_loss");
	    else{
		Fprintf(stderr,
		    "iir_filt: Stop_band_loss not specified - exiting.\n");
		exit(1);
	    }

	    if( stop_loss <= 0.){
		Fprintf(stderr, 
		    "iir_filt: stop band loss must be > 0 - exiting.\n");
		exit(1);
	    }

	    break;

	case ELLIPTICAL:
	    if(symtype("band_edges") != ST_UNDEF){
	       num_of_band_edges = getsym_fa("band_edges",band_edges,(int)4);
	    }
	    else{
		Fprintf(stderr,
	        "iir_filt: band_edges not specified - exiting.\n");
		exit(1);
	    }

	    if(symtype("pass_band_loss") != ST_UNDEF)
		pass_loss = (float)getsym_d("pass_band_loss");
	    else{
		Fprintf(stderr,
		    "iir_filt: Pass_band_loss not specified - exiting.\n");
		exit(1);
	    }

	    if( pass_loss <= 0.){
		Fprintf(stderr, 
		    "iir_filt: Pass band loss must be > 0 - exiting.\n");
		exit(1);
	    }

	    if(symtype("stop_band_loss") != ST_UNDEF)
		stop_loss = (float)getsym_d("stop_band_loss");
	    else{
		Fprintf(stderr,
		    "iir_filt: Stop_band_loss not specified - exiting.\n");
		exit(1);
	    }

	    if( stop_loss <= 0.){
		Fprintf(stderr, 
		    "iir_filt: stop band loss must be > 0 - exiting.\n");
		exit(1);
	    }

	    break;
	default:
	    Fprintf(stderr, 
"iir_filt: Invalid filter design method (%s) specified - exiting\n", 
	    filt_mthd);
	    exit(1);
    }

/* Sanity check on specified parameters */


	if(sf <= 0.){
	    Fprintf(stderr, "iir_filt: Sampling Frequency must be > 0.\n");
	    exit(1);
	}

	switch(filter_poly){
	    int num_of_freqs;

	    case BUTTERWORTH:
	    case CHEBYSHEV1:
	    case CHEBYSHEV2:

		num_of_freqs = 2;
		if(filter_response == FILT_HP || filter_response == FILT_LP)
		    num_of_freqs = 1;

		for(i = 0;  i < num_of_freqs; i++){
		    if(three_dBs[i] >= sf/2){
			Fprintf(stderr, 
"iir_filt: Critical frequency (%g Hz) >= (sampling frequency)/2 - exiting\n",
			three_dBs[i]);
			exit(1);
		    }
		}

		for(i = 0; i < num_of_freqs; i++){
		    if(three_dBs[i] <= 0.){
			Fprintf(stderr,
     "iir_filt: Critical frequency (%g Hz) <= 0 - exiting.\n", three_dBs[i]);
			exit(1);
		    }
		}
		break;

	    case ELLIPTICAL:

		num_of_freqs = 4;
		if(filter_response == FILT_HP || filter_response == FILT_LP)
		    num_of_freqs = 2;


		for(i = 0;  i < num_of_freqs; i++){
		    if(band_edges[i] >= sf/2){
			Fprintf(stderr, 
"iir_filt: Critical frequency (%g Hz) >= (sampling frequency)/2 - exiting\n",
			band_edges[i]);
			exit(1);
		    }
		}

		for(i = 0;  i < num_of_freqs; i++){
		    if(band_edges[i] <= 0.){
			Fprintf(stderr, 
     "iir_filt: Critical frequency (%g Hz) <= 0 - exiting\n", band_edges[i]);
			exit(1);
		    }
		}
		break;
	    default:
		Fprintf(stderr, 
	           "iir_filt: Invalid filter design type (%s)\n", filt_mthd);
		exit(1);
	}

	if(gain <= 0){
	    Fprintf(stderr, "iir_filt: Gain must be > 0 - exiting.\n");
	    exit(1);
	}

	if(filter_order <= 0){
          Fprintf(stderr, "iir_filt: Filter order must be > 0 - exiting.\n");
	    exit(1);
	}

	switch(filter_response){
	    case FILT_LP:
	    case FILT_HP:
		switch(filter_poly){
		    case ELLIPTICAL:
			if(band_edges[1] <= band_edges[0]){
			    Fprintf(stderr,
"iir_filt: Band edges must be ordered from smallest to largest - exiting.\n");
			    exit(1);
			}

			break;
		    case BUTTERWORTH:
		    case CHEBYSHEV1:
		    case CHEBYSHEV2:
			break;
		    default:
			break;
		}
		break;
	    case FILT_BP:
	    case FILT_BS:
		switch(filter_poly){
		    case ELLIPTICAL:
			if(band_edges[1] <= band_edges[0]  ||
			    band_edges[2] <= band_edges[1] ||
			    band_edges[2] <= band_edges[3]){
				Fprintf(stderr,
"iir_filt: Band edges must be ordered from smallest to largest - exiting.\n");
			    exit(1);
			 }
			 break;
		    case BUTTERWORTH:
		    case CHEBYSHEV1:
		    case CHEBYSHEV2:
			if(three_dBs[1] <= three_dBs[0]){
			    Fprintf(stderr,
"iir_filt: Band edges must be ordered from smallest to largest - exiting.\n");
			    exit(1);
			}
			break;
		    default:
			break;
		}
		break;
	    default:
		break;
	}


/*DEBUG*/
    if(debug_level > 0){
	Fprintf(stderr, "\n");
	Fprintf(stderr, "\tiir_filt: LISTING OF SPECIFIED PARAMETERS\n\n");
	Fprintf(stderr, 
	    "\n\niir_filt: Samp. freq. = %f, pass band gain = %f\n",
	    sf, gain);
	Fprintf(stderr, 
	    "iir_filt: filt_method = %s, filt_type = %s\n",
	    filt_mthd, filt_resp_type);
        Fprintf(stderr,
		"\niir_filt: order = %d, 3 dB freqs = %f and %f\n",
		filter_order, three_dBs[0], three_dBs[1]);
        Fprintf(stderr,
"iir_filt: pass band loss = %f dB,\nband edges = %f, %f, %f and %f Hz\n",
    pass_loss, band_edges[0], band_edges[1], band_edges[2], band_edges[3]);
	    Fprintf(stderr,"iir_filt: stop band loss = %f dB\n", stop_loss);
    }
/*END DEBUG*/



/* Create output header */

    oh = new_header (FT_FILT);

/* Allocate space for filter data record */