pkmc_filt.c

speech signal process tools

/*
 * 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-1993 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: Parks-McClellan FIR design
 *
 */

static char *sccs_id = "@(#)pkmc_filt.c	1.7 7/30/97 ERL";

#include <math.h>
#include <string.h>
#include <stdio.h>
#include <errno.h>
#include <stdlib.h>
#include <esps/esps.h>
#include <esps/constants.h>
#include <esps/fea.h>
#include <esps/feafilt.h>
#include <esps/feasd.h>
#include <esps/strlist.h>

#define SYNTAX USAGE ("pkmc_filt [-P param_file][-x debug_level] filt_file");
#define ERROR_EXIT(text) {(void) fprintf(stderr, "%s: %s - exiting\n", \
					 argv[0], text); exit(1);}

/* local definitions*/
#define MAX_ITER 25
#define MAX_NBANDS 10
#define SMALL 0.00000001
#define OPPSIGN(A, B) ((A > 0 && B < 0) || (A < 0 && B > 0))

struct diff{
  double err;    /* value */
  int idx;       /* position on the grid */
};

struct edges{
  int idx1;
  int idx2;
};

struct spec{
  double edge[2];
  double des;
  double wt;
};

int debug_level = 0;

#define MULTIBAND 0
#define DIFFERENTIATOR 1
#define HILBERT 2
static char *jtype_code[] = {"MULTIBAND", "DIFFERENTIATOR", "HILBERT"};

#define CASE1 1
#define CASE2 2
#define CASE3 3
#define CASE4 4


/* system functions */
double pow(), fabs();
void perror();


/*ESPS functions*/
int getopt();
char *get_cmd_line();

/* local function definitions*/
void find_del_interpolate();
void find_extrema();
int check_extrema();
void write_outsd();
void get_coef();
void remez_exchange();
void setup();
void find_impulse();


char *Date="7/30/97";
char *Version="1.7";
char *Program="@(#)pkmc_filt.c	1.7";

void
main( argc, argv)
     int argc;
     char **argv;
{
  extern char *optarg;
  extern optind;
  FILE *fpout;
  struct header *oh;
  struct feafilt *frec;
  char *param_file = NULL;
  char *filt_file = NULL;
  filtdesparams *fdp;
  char *str, stri[3], *str_edge1, *str_edge2, *str_des, *str_wt;
  char *stu = "band";
  char *stl1 = "edge1";
  char *stl2 = "edge2";
  char *std = "des";
  char *stw = "wt";
  char *filt_type = NULL;
  int c, i, j, err=0;
  double *des, *cx, *wt;
  int Ngrid, nbands;
  float *hd_edges, *hd_des, *hd_wt;
  struct edges *bands;
  struct spec *firspecs;
  int Ncoef, Ntaps, jtype, fcase;
  double *A, fs, *coef, *h, *f, del, delay = 0;
  float *ripple_db;

/*
 * 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
          }
    }

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

/* 
 *Read parameter file and get design parameters 
 */
  
  if (debug_level)
    Fprintf (stderr,"%s: Reading parameter file %s\n", argv[0],param_file);
  
  if (read_params (param_file, SC_NOCOMMON, (char *)NULL) != 0)
    ERROR_EXIT("read_params could not read the params file.");

  /* filt_type: MULTIBAND, DIFFERENTIATOR, HILBERT */
  if (symtype("filt_type") != ST_UNDEF)
    filt_type = getsym_s ("filt_type");
  else ERROR_EXIT("Filt_type not specified in params file.");
  if((jtype = lin_search(jtype_code, filt_type)) == -1 ){
    Fprintf(stderr, "%s: Invalid filt_type (%s) specified - exiting.\n",
	    argv[0], filt_type);
    exit(1);
  }

  /* filter length parameter: even or odd number */
  if (symtype("filt_length") != ST_UNDEF)
    Ntaps = getsym_i ("filt_length");
  else ERROR_EXIT("Filt_length not specified in params file.");

  /* Ngrid parameter */
  if (symtype("ngrid") != ST_UNDEF)
    Ngrid = getsym_i ("ngrid");
  else Ngrid = 2 * Ntaps;

  /* sampling frequency parameter*/
  if( symtype("samp_freq") != ST_UNDEF)
    fs = getsym_d("samp_freq");
  else ERROR_EXIT("Samp_freq not specified in params file.");
  
  /* nbands parameter */
  if( symtype("nbands") != ST_UNDEF)
    nbands = getsym_i("nbands");
  else ERROR_EXIT("Nbands not specified in params file.");
 
  firspecs = (struct spec *) malloc( sizeof(struct spec) * nbands);
  /* band edges parameters, must be < fs/2 */
  str = (char *) malloc(sizeof(char)*20);
  str_edge1 = (char *) malloc(sizeof(char)*20);
  str_edge2 = (char *) malloc(sizeof(char)*20);
  str_des = (char *) malloc(sizeof(char)*20);
  str_wt = (char *) malloc(sizeof(char)*20);
  for( i = 1 ;i<=nbands; i++){
    str = (char *) strcpy( str, stu);
    sprintf( stri, "%d", i);
    str = (char *) strcat(str, stri);
    str = (char *) strcat(str, "_");
    /* find edges */
    str_edge1 = (char *) strcpy (str_edge1, str);
    str_edge1 = (char *) strcat(str_edge1, stl1);
    str_edge2 = (char *) strcpy (str_edge2, str);
    str_edge2 = (char *) strcat(str_edge2, stl2);
    if( symtype(str_edge1) != ST_UNDEF)
      firspecs[i-1].edge[0] = getsym_d(str_edge1);
    else err = 1;
    if( symtype(str_edge2) != ST_UNDEF)
      firspecs[i-1].edge[1] = getsym_d(str_edge2);
    else err = 1;
    if(err){
      Fprintf(stderr,"%s: missing band edges for band %d\n", argv[0], i);
      exit(1);
    }
    if( firspecs[i-1].edge[1] > fs/2)
      ERROR_EXIT("bandedge can not be greater than the Nyquist rate");
    /* find desired value */
    str_des = (char *) strcpy (str_des, str);
    str_des = (char *) strcat (str_des, std);
    if( symtype(str_des) != ST_UNDEF)
      firspecs[i-1].des = getsym_d(str_des);
    else{
      Fprintf(stderr,"%s: missing desired value for band %d\n", argv[0], i);
      exit(1);
    }
    str_wt = (char *) strcpy (str_wt, str);
    str_wt = (char *) strcat (str_wt, stw);
    if( symtype(str_wt) != ST_UNDEF)
      firspecs[i-1].wt = getsym_d(str_wt);
    else{
      Fprintf(stderr,"%s: missing weight value for band %d\n", argv[0], i);
      exit(1);
    }
  }

  /* Check that upper limit of last band == (sampling rate)/2.
     DB added to avoid numeric-error/core-dump that can occur.
   */
  if( firspecs[nbands-1].edge[1] != fs/2.) {
      Fprintf(stderr, "Upper bandedge of last band (%.2f) != sampling frequency / 2 (%.2f).\n\tPlease correct parameters and try again.\n",
		 firspecs[nbands-1].edge[1], fs/2);
      exit(1);
    }

  /* set up output generic header items */
  hd_edges = (float *) malloc(sizeof(float) * 2* nbands);
  hd_des = (float *) malloc(sizeof(float) *nbands);
  hd_wt = (float *) malloc(sizeof(float) *nbands);
  for(i=0, j=0; i<nbands;i++, j+=2){ 
    hd_edges[j] = firspecs[i].edge[0];
    hd_edges[j+1] = firspecs[i].edge[1];
    hd_des[i] = firspecs[i].des;
    hd_wt[i] = firspecs[i].wt;
  }

  bands = (struct edges *) malloc( sizeof(struct edges) * nbands);

  des = (double *) malloc(sizeof(double) * Ngrid);
  wt = (double *) malloc(sizeof(double) * Ngrid);
  cx = (double *) malloc(sizeof(double) * Ngrid);
  A = (double *) malloc (sizeof(double)* Ngrid);
  coef = (double *) malloc (sizeof(double)* Ntaps);
  h = (double *) malloc( sizeof(double) * Ntaps);
  f = (double *) malloc(sizeof(double) *Ngrid);
  ripple_db = ( float *) malloc(sizeof(float) * nbands);

  setup( f, cx, des, wt, &Ngrid, nbands, bands, Ntaps, &Ncoef, firspecs, jtype, &fcase, fs );
  remez_exchange( Ncoef, Ngrid, cx, des, wt, nbands, bands, A, coef, &del);

  for(i=0;i<nbands;i++) ripple_db[i] = 20*log10( firspecs[i].des + del/wt[i] );
  if(debug_level > 0)
    for(i=0;i<nbands;i++) 
      Fprintf(stderr,"%s: ripple in band %d is %f dB\n", argv[0],  
	      i, ripple_db[i]);
  
  if(debug_level == 15521) {
    write_outsd(A,Ngrid,"pkmc_A");
    write_outsd(des,Ngrid,"pkmc_des");
  }

  find_impulse(coef, Ncoef, h, fcase, Ntaps );

  fdp = (filtdesparams *) calloc( 1, sizeof(filtdesparams));
  fdp->define_pz = NO;
  fdp->method = PARKS_MC;
  fdp->nbands = nbands;
  fdp->bandedges = hd_edges;
  fdp->wts = hd_wt;
  fdp->gains = hd_des;
  fdp->g_size = Ngrid;
  delay = (double) (Ntaps-1) / 2;

  oh = new_header(FT_FEA);
  (void) strcpy (oh->common.prog, argv[0]);
  (void) strcpy (oh->common.vers, Version);
  (void) strcpy (oh->common.progdate, Date);
  oh->common.tag = NO;
  add_comment(oh, get_cmd_line(argc, argv));
  init_feafilt_hd( oh, (long)(Ntaps), 0L, fdp);
  (void)add_genhd_f("band_edges", hd_edges, 2*nbands, oh);
  (void)add_genhd_f("desired_value", hd_des, nbands, oh);
  (void)add_genhd_f("desired_weight", hd_wt, nbands, oh);
  (void)add_genhd_f("ripple_dB", ripple_db, nbands, oh);
  (void)add_genhd_d("samp_freq", &fs, 1, oh);
  (void)add_genhd_d("delay_samples", &delay, 1, oh);

  (void) write_header ( oh, fpout );
  frec = allo_feafilt_rec(oh);
  *frec->num_size = Ntaps;
  for(i=0;i<Ntaps;i++) frec->re_num_coeff[i] = h[i];
  put_feafilt_rec( frec, oh, fpout);


  /* clean up*/
  free(firspecs);
  free(str);
  free(str_edge1);
  free(str_wt);
  free(str_edge2);
  free(str_des);
  free(hd_edges);
  free(hd_des);
  free(hd_wt);
  free(bands);
  free(des);
  free(wt);
  free(cx);
  free(A);
  free(coef);
  free(h);
  free(f);
  free(ripple_db);
  free(fdp);
  
  exit(0);
}

void
find_impulse(coef, Ncoef, h, fcase, Ntaps)
  double *coef, *h;
  int Ncoef, fcase, Ntaps;
{
  int i;
  switch(fcase){
  case CASE1:
    if( Ncoef != (Ntaps - 1)/2 + 1) 
      Fprintf(stderr,"ERROR: number of coefficients mismatch.\n");
    for(i=0; i< Ncoef-1; i++) h[i] = 0.5 * coef[Ncoef- 1 - i];
    h[Ncoef-1] = coef[0];
    for(i=Ncoef; i< Ntaps; i++) h[i] = h[Ntaps -i -1];
    break;
  case CASE2:
    if( Ncoef != Ntaps / 2)
      Fprintf(stderr,"ERROR: number of coefficients mismatch.\n");
    if( Ncoef > 1){
      h[0] = 0.25 * coef[Ncoef - 1];
      for(i=1; i< Ncoef-1; i++) h[i] = 0.25* (coef[Ncoef-1-i] + coef[Ncoef-i]);
      h[Ncoef - 1] = 0.5*(coef[0] + 0.5 * coef[1]);
    }
    else
      h[0] = 0.5*coef[0];
    for(i=Ncoef; i< Ntaps; i++) h[i] = h[Ntaps -i -1];
    break;
  case CASE3:
    if( Ncoef != (Ntaps-1) / 2)
      Fprintf(stderr,"ERROR: number of coefficients mismatch.\n");
    if( Ncoef >= 3){
      h[0] = .25 * coef[Ncoef -1];
      h[1] =  0.25 * coef[Ncoef -1 ];
      for(i=2; i< Ncoef-1; i++) h[i]= 0.25*(coef[Ncoef-i-1] - coef[Ncoef-i+1]);
      h[Ncoef-1] = 0.5 *(coef[0] -0.5*coef[2]);
    }
    else if(Ncoef==2){
      h[0] = 0.25 * coef[1];
      h[1] = 0.5 * coef[0];
    }
    else if(Ncoef==1)
      h[0] = 0.5 * coef[0];
    
    h[Ncoef] = 0;
    for(i=Ncoef+1; i< Ntaps; i++) h[i] = -h[Ntaps -i-1];
    break;
  case CASE4:
    if( Ncoef != Ntaps/ 2)
      Fprintf(stderr,"ERROR: number of coefficients mismatch.\n");
    if( Ncoef > 1){
      h[0] = 0.25 * coef[Ncoef-1];
      for(i=1; i< Ncoef-1; i++) h[i]= 0.25*(coef[Ncoef-i-1]-coef[Ncoef-i]);
      h[Ncoef-1] = 0.5 * (coef[0] - 0.5 * coef[1]);
    }
    else
      h[0] = 0.5 * coef[0];

    for(i=Ncoef; i<Ntaps; i++) h[i] = - h[Ntaps-i-1];
    break;
  }
}

void
setup( w, cx, des, wt, Ngrid, nbands, bands, Ntaps, Ncoef, 
	   firspecs, jtype, fcase, nfs )
     double *w, *cx, *des, *wt;
     int *Ngrid, nbands, *Ncoef, Ntaps;
     struct edges *bands;
     struct spec *firspecs;
     int jtype, *fcase;
     double nfs;        /* Nyquist rate, not sampling rate */
{
  double pi2, dfreq, minf, maxf, f, desired, weight;
  int k,m,i,j,d;
  

  pi2 = 2 * PI;
  dfreq = 0.5 / (*Ngrid-1);

  if( jtype == MULTIBAND ){
    if( 2*(Ntaps/2) != Ntaps) { 
      *fcase = CASE1;
      minf = 0;
      maxf = 0.5;