mbs_dist.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) 1998 Entropic, Inc.  All rights reserved."
 *
 * The copyright notice above does not evidence any actual or intended
 * publication of this source code.
 *
 * Written by:  Rod Johnson, based on code supplied by Wonho Yang
 * Checked by:
 * Revised by:
 *
 * Brief description:
 *   compute modified Bark spectral distortion
 *
 * Ref:  W Yang, M Benbouchta, & R Yantorno, ICASSP '98 Proc
 * (May 1998) 541-544.
 */

static char *sccs_id = "@(#)mbs_dist.c	1.1	9/10/98	ERL";


/* INCLUDE FILES */

#include <stdlib.h>
#include <math.h>
#include <esps/esps.h>
#include <esps/fea.h>
#include <esps/feaspec.h>

/* LOCAL CONSTANTS */

#define VERSION	("1.1")		/* must be 7 char or less */
#define DATE	("9/10/98")	/* must be 25 char or less */
#define PROG	("mbs_dist")	/* must be 15 char or less */

#define FRAME	(320)		/* frame size (samples) */
#define BSIZE	(18)		/* number of Bark frequencies */
#define PHSIZE	(13)		/* sizes of arrays used for interpolating
				 * values in phons (see dbtophon) */

/* LOCAL MACROS */

#define SYNTAX \
USAGE("mbs_dist [-a] ] [-{pr} range1 [-{pr} range2]] [-x debug_level] [-A]\n" \
      "\t[-P paramFile] inFile1.spec inFile2.spec [outFile.fea]")

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

#define REQUIRE(test, text) {if (!(test)) ERROR(text)}

#define TRYALLOC(type, num, var, txt) \
{if (((var) = (type *) malloc((num)*sizeof(type))) == NULL) \
 {(void) fprintf(stderr, "%s: can't allocate memory--%s - exiting\n", \
		 PROG, (txt)); exit(1);}}

#define GET_GENHD_CODED(name, hd, var, txt) \
{if (genhd_type((name), NULL, (hd)) != CODED) \
 {(void) fprintf(stderr, "%s: header item \"%s\" in %s undefined or" \
		 " wrong type - exiting\n", PROG, (name), (txt)); exit(1);} \
 {(var) = *get_genhd_s((name), (hd));}}

#define STRCMP(a, op, b) (strcmp((a), (b)) op 0)

/* SYSTEM FUNCTIONS AND VARIABLES */

extern int	getopt();		/* parse command options */
extern int	optind;			/* used by getopt */
extern char	*optarg;		/* used by getopt */
long		lseek();		/* used to check for piped input */

/* ESPS FUNCTIONS AND VARIABLES */

char		*get_cmd_line(int argc, char **argv);
int		debug_level = 0;

/* LOCAL TYPEDEFS AND STRUCTURES */

/* This enum type is common to bs_dist and mbs_dist.
 * Keep them consistent.
 */

typedef enum {
    BSD = 1,    MBSD
}		d_type;			/* cf. d_types below */

/* LOCAL FUNCTION DECLARATIONS */

void	spec_to_pwr(int num,
		    float *src, float *src_im, int src_type, double *dest);
double	frame_pwr(int num, double *pwr_spec, double df);
double	*init_freqs(int num, double fmax);
double	*init_spread_func();
double	*init_thresh(int num, double *freqs);
void	normalize(int num, double *pwr_spec, double mean_pwr, double factor);
int	check_frame(double pwr1, double pwr2, double XTHRESH, double YTHRESH);
void	bk_frq(int num, double *freqs, double *pwr_spec, double *bk_spec);
void	spread(double *bk_spec, double *spread_spec);
void	dbtophon(double *spread_spec, double *phon_spec);
void	phontoson(double *phon_spec, double *sone_spec);
double	measure(double *orig, double *dist, double *noise);
double	sfm(int num, double *pwr_spec);
void	thresh2(double alpha, double *spread_spec, double *noise_thr);

void	print_farray(int filenum, long recnum,
		     char *text, long numelem, float *data);
void	print_darray(int filenum, long recnum,
		     char *text, long numelem, double *data);

/* STATIC (LOCAL) GLOBAL VARIABLES */

/* This string array definition is common to bs_dist and mbs_dist.
 * Keep them consistent.
 */

static char	*d_types[] = {
    "NONE", "BSD",  "MBSD",
    NULL
};				/* cf. d_type above */

static char	*noyes[] = {"NO", "YES", NULL};

static double	*FREQ;		/* frequency scale */

static int	BARK[BSIZE+1] = {
    0, 100, 200, 300, 400, 510, 630, 770, 920, 1080,
    1270, 1480, 1720, 2000, 2320, 2700, 3150, 3700, 4400
};				/* Bark frequency */

static double	*S;		/* spreading function */

static double	*ABS_TH;	/* absolute hearing threshold */

static double	eqlcon[PHSIZE][BSIZE-3] = {
{12.0,7.0,4.0,1.0,0.0,0.0,0.0,-0.5,-2.0,-3.0,-7.0,-8.0,-8.5,-8.5,-8.5},
{20.0,17.0,14.0,12.0,10.0,9.5,9.0,8.5,7.5,6.5,4.0,3.0,2.5,2.0,2.5},
{29.0,26.0,23.0,21.0,20.0,19.5,19.5,19.0,18.0,17.0,15.0,14.0,13.5,13.0,13.5},
{36.0,34.0,32.0,30.0,29.0,28.5,28.5,28.5,28.0,27.5,26.0,25.0,24.5,24.0,24.5},
{45.0,43.0,41.0,40.0,40.0,40.0,40.0,40.0,40.0,39.5,38.0,37.0,36.5,36.0,36.5},
{53.0,51.0,50.0,49.0,48.5,48.5,49.0,49.0,49.0,49.0,48.0,47.0,46.5,45.5,46.0},
{62.0,60.0,59.0,58.0,58.0,58.5,59.0,59.0,59.0,59.0,58.0,57.5,57.0,56.0,56.0},
{70.0,69.0,68.0,67.5,67.5,68.0,68.0,68.0,68.0,68.0,67.0,66.0,65.5,64.5,64.5},
{79.0,79.0,79.0,79.0,79.0,79.0,79.0,79.0,78.0,77.5,76.0,75.0,74.5,73.0,73.0},
{89.0,89.0,89.0,89.5,90.0,90.0,90.0,89.5,89.0,88.5,87.0,86.0,85.5,84.0,83.5},
{100.0,100.0,100.0,100.0,100.0,99.5,99.0,99.0,98.5,98.0,96.0,95.0,94.5,93.5,93.0},
{112.0,112.0,112.0,112.0,111.0,110.5,109.5,109.0,108.5,108.0,106.0,105.0,104.5,103.0,102.5},
{122.0,122.0,121.0,121.0,120.5,120.0,119.0,118.0,117.0,116.5,114.5,113.5,113.0,111.0,110.5}
};

static double	phons[PHSIZE] = {
0.0,10.0,20.0,30.0,40.0,50.0,60.0,70.0,80.0,90.0,100.0,110.0,120.0
};


/* MAIN PROGRAM */

int
main(int    argc,
     char   **argv)
{
    int		    ch;			/* command-line option letter */

    long	    i, j, k, n;		/* loop counters and array indices */
    long	    num_proc;		/* number of "processed" frames */

    int		    a_specified = NO;	/* a-option specified? (output only
					 * average final distortion with no
					 * output file?) */

    char	    *r_arg1 = NULL;	/* first r-option argument */
    char	    *r_arg2 = NULL;	/* second r-option argument */
    long	    startrec1;		/* first record to process from
					 * first input file */
    long	    startrec2;		/* first record to process from
					 * second input file */
    int		    startrecs[2]; 	/* startrec1 and startrec2 as
					 * an array */
    long	    nan;		/* number of records to process,
					 * but 0 means range extends to
					 * end of file */

    int		    A_specified = NO;	/* A-option specified? (output
					 * average final distortion?) */

    char	    *param_name = NULL;	/* parameter file name */

    char	    *inname1;		/* first input file name */
    FILE	    *infile1;		/* first input stream */
    struct header   *inhd1;		/* first input file header */
    struct feaspec  *inrec1;		/* input FEA_SPEC record from
					 * infile1 */
    float	    *inbuf1;		/* input buffer--data from inrec1 */
    float	    *inbufim1;		/* data from inrec1 (imag part) */
    long	    *intag1;		/* pointer to tag in input record */
    int		    freq_format1;	/* freq_format of infile1 */
    int		    spec_type1;		/* spec_type of infile1 */
    int		    contin1;		/* is first input spectrum a
					 * continuous density? */
    double	    src_sf1;		/* src_sf in infile1 */
    double	    sf1;		/* sf in infile1 */
    int		    num_freqs1;		/* num_freqs in infile1 */
    double	    record_freq1;	/* record_freq in infile1 */
    FILE	    *tmpfile1 = NULL;	/* temp file to allow two passes
					 * over data when infile1 is pipe */
    struct header   *tmphd1 = NULL;	/* header for use with tmpfile1 */
    double	    tot_pwr1;		/* power in frame of infile1 */

    char	    *inname2;		/* second input file name */
    FILE	    *infile2;		/* second input stream */
    struct header   *inhd2;		/* second input file header */
    struct feaspec  *inrec2;		/* input FEA_SPEC record from
					 * infile2 */
    float	    *inbuf2;		/* input buffer--data from inrec2 */
    float	    *inbufim2;		/* data from inrec2 (imag part) */
    long	    *intag2;		/* pointer to tag in input record */
    int		    freq_format2;	/* freq_format of infile2 */
    int		    spec_type2;		/* spec_type of infile2 */
    int		    contin2;		/* is second input spectrum a
					 * continuous density? */
    double	    src_sf2;		/* src_sf in infile2 */
    double	    sf2;		/* sf in infile2 */
    int		    num_freqs2;		/* num_freqs in infile2 */
    double	    record_freq2;	/* record_freq in infile2 */
    FILE	    *tmpfile2 = NULL;	/* temp file to allow two passes
					 * over data when infile2 is pipe */
    struct header   *tmphd2 = NULL;	/* header for use with tmpfile2 */
    double	    tot_pwr2;		/* power in frame of infile2 */

    double	    delta_f;		/* interval between frequencies */
    int		    num_freqs;		/* number of frequencies used in
					 * computation */
    double	    fmax;		/* largest frequency used in
					 * computation */
    int		    tag_match;	 	/* are input files tagged files with
					 * the same src_sf? */
    int		    rec_freq_match;	/* do input files have the same
					 * non-zero record_freq? */

    double	    max_xpower; 	/* max frame power in infile1 */
    double	    max_ypower; 	/* max frame power in infile2 */
    double	    mean_xpower;	/* total power in infile1 */
    double	    mean_ypower;	/* total power in infile2 */
    double	    XTHRESHOLD;		/* original-speech threshold for
					 * including frame in overall
					 * average distortion */
    double	    YTHRESHOLD;		/* processed-speech threshold for
					 * including frame in overall
					 * average distortion */
    double	    *PSX;		/* power spectrum from infile1 */
    double	    *PSY;		/* power spectrum from infile2 */
    double	    BX[BSIZE];		/* Bark spectrum from infile1 */
    double	    BY[BSIZE];		/* Bark spectrum from infile2 */
    double	    CX[BSIZE];		/* spread Bark spect from infile1 */
    double	    CY[BSIZE];		/* spread Bark spect from infile2 */
    double	    CNMT[BSIZE];	/* spread Bark spect noise masking
					 * threshold */
    double	    PX[BSIZE-3]; 	/* spread Bark spect from infile1
					 * in phons */
    double	    PY[BSIZE-3]; 	/* spread Bark spect from infile2
					 * in phons */
    double	    PN[BSIZE-3]; 	/* spread Bark spect of noise
					 * in phons */
    double	    SX[BSIZE-3]; 	/* spread Bark spect from infile1
					 * in sones */
    double	    SY[BSIZE-3]; 	/* spread Bark spect from infile2
					 * in sones */
    double	    SN[BSIZE-3]; 	/* spread Bark spect of noise
					 * in sones */

    double	    alpha;		/* noise masking threshold */

    int 	    included;		/* include current frame in overall
					 * average distortion? */

    double	    frame_distortion;	/* distortion for one frame */
    double	    sum_distortion;	/* distortion summed over frames */

    char	    *outname;		/* output file name */
    FILE	    *outfile;		/* output stream */
    struct header   *outhd;		/* output file header */
    struct fea_data *outrec;		/* output file record */
    char	    *distfld;		/* name of frame distortion field */
    char	    *flagfld;		/* name of "included" flag field */
    float	    *outdata;		/* pointer to frame distortion field
					 * in output record */
    short	    *outflag;		/* pointer to "included" flag field
					 * in output record */
    long	    *outtag;		/* pointer to tag in output record */
    long	    *start;		/* pointer to "start" in output hdr */

    /*
     * PARSE COMMAND-LINE OPTIONS.
     */

    while ((ch = getopt(argc, argv, "ap:r:x:AP:")) != EOF)
	switch (ch) {

	case 'a':
	    a_specified = YES;
	    break;
	case 'p':	/* -p is a synonym for -r */
	    /* FALL THROUGH */
	case 'r':
	    if (r_arg1 == NULL) {
		r_arg1 = optarg;
	    } else if (r_arg2 == NULL) {
		r_arg2 = optarg;
	    } else {
		fprintf(stderr, "%s: too may -r options.\n", PROG);
		SYNTAX;
	    }
	    break;
	case 'x':
	    debug_level = atoi(optarg);
	    break;
	case 'A':
	    A_specified = YES;
	    break;
	case 'P':
	    param_name = optarg;
	    break;
	default:
	    SYNTAX;
	    break;
	}

    /*
     * PROCESS FILE NAMES.
     */

    if (argc - optind > (a_specified ? 2 : 3)) {

	fprintf(stderr, "%s: too many file names.\n", PROG);
	SYNTAX;
    }

    if (argc - optind < (a_specified ? 2 : 3)) {

	fprintf(stderr, "%s: not enough file names.\n", PROG);
	SYNTAX;
    }

    inname1 = argv[optind++];
    inname2 = argv[optind++];

    REQUIRE(STRCMP(inname1, !=, "-") || STRCMP(inname2, !=, "-"),
	    "input files both standard input");

    if (!a_specified) {

	outname = argv[optind++];

	if (STRCMP(outname, !=, "-")) {

	    REQUIRE(STRCMP(inname1, !=, outname),
		    "output file same as first input file");
	    REQUIRE(STRCMP(inname2, !=, outname),
		    "output file same as second input file");
	}

	if (A_specified) {

	    REQUIRE(STRCMP(outname, !=, "-"),
		    "-A specified, but output file is standard output");
	}
    }

    /*
     * OPEN AND CHECK INPUT FILES.
     */

    inname1 = eopen(PROG, inname1, "r", FT_FEA, FEA_SPEC, &inhd1, &infile1);
    inname2 = eopen(PROG, inname2, "r", FT_FEA, FEA_SPEC, &inhd2, &infile2);

    if (debug_level >= 1) {

	fprintf(stderr, "%s: first input file = \"%s\"\n", PROG, inname1);
	fprintf(stderr, "%s: second input file = \"%s\"\n", PROG, inname2);
    }

    GET_GENHD_CODED("freq_format", inhd1, freq_format1, inname1);
    REQUIRE(freq_format1 == SPFMT_SYM_EDGE,
	    "freq_format in first input file is not SYM_EDGE");
    GET_GENHD_CODED("freq_format", inhd2, freq_format2, inname2);
    REQUIRE(freq_format1 == SPFMT_SYM_EDGE,
	    "freq_format in second input file is not SYM_EDGE");

    GET_GENHD_CODED("spec_type", inhd1, spec_type1, inname1);
    GET_GENHD_CODED("spec_type", inhd2, spec_type2, inname2);

    GET_GENHD_CODED("contin", inhd1, contin1, inname1);
    GET_GENHD_CODED("contin", inhd2, contin2, inname2);
    REQUIRE((contin1 && contin2) || (!contin1 && !contin2),
	    "one input contains a spectral density, "
	    "and one contains discrete powers");

    src_sf1 = get_genhd_val("src_sf", inhd1, 0.0);
    if (src_sf1 == 0.0)
	src_sf1 = get_genhd_val("sf", inhd1, 0.0);
    src_sf2 = get_genhd_val("src_sf", inhd2, 0.0);
    if (src_sf2 == 0.0)
	src_sf2 = get_genhd_val("sf", inhd2, 0.0);

    sf1 = get_genhd_val("sf", inhd1, 0.0);
    REQUIRE(sf1 > 0.0,
	    "sf missing or not positive in first input file");
    sf2 = get_genhd_val("sf", inhd2, 0.0);
    REQUIRE(sf2 > 0.0,
	    "sf missing or not positive in second input file");

    num_freqs1 = (long) get_genhd_val("num_freqs", inhd1, -1.0);
    REQUIRE(num_freqs1 > 0,
	    "num_freqs missing or not positive in first input file");
    num_freqs2 = (long) get_genhd_val("num_freqs", inhd2, -1.0);
    REQUIRE(num_freqs2 > 0,
	    "num_freqs missing or not positive in second input file");

    delta_f = sf1 / (2.0 * (num_freqs1 - 1));
    REQUIRE(delta_f < 100.0,
	    "frequency spacing not small enough");
    /* 100.0 is the minimum width in Hz of the 1-bark frequency bins
     * (cf. array BARK).  This just assures that no bins turn up empty.
     * Considerably smaller values are preferable in practice---e.g.
     * 8000.0/1024.
     */
    REQUIRE(sf2 / (2.0 * (num_freqs2 - 1)) < 1.0001*delta_f
	    && sf2 /(2.0 * (num_freqs2 - 1)) > 0.9999*delta_f,
	    "inconsistent frequency spacing in the two input files");

    num_freqs = 1 + ROUND(4000.0/delta_f);
    fmax = (num_freqs - 1) * delta_f;

    if (num_freqs > num_freqs1) {

	num_freqs = num_freqs1;
	fmax = 0.5 * sf1;
    }

    if (num_freqs > num_freqs2) {

	num_freqs = num_freqs2;
	fmax = 0.5 * sf2;
    }

    record_freq1 = get_genhd_val("record_freq", inhd1, 0.0);
    record_freq2 = get_genhd_val("record_freq", inhd2, 0.0);

    if (debug_level >= 1)
	fprintf(stderr, "%s: in %s---\n"
		"\tspec_type = \"%s\", contin = %d, src_sf = %g,\n"
		"\tsf = %g, num_freqs = %ld, record_freq = %g,\n"
		"\tfile is%stagged\n",
		PROG, inname1,
		sptyp_names[spec_type1],
		contin1,
		src_sf1,
		sf1,
		num_freqs1,
		record_freq1,
		(inhd1->common.tag) ? " " : " not ");

    if (debug_level >= 1)
	fprintf(stderr, "%s: in %s---\n"
		"\tspec_type = \"%s\", contin = %d, src_sf = %g,\n"
		"\tsf = %g, num_freqs = %ld, record_freq = %g,\n"
		"\tfile is%stagged\n",
		PROG, inname2,
		sptyp_names[spec_type2],
		contin2,
		src_sf2,
		sf2,
		num_freqs2,
		record_freq2,
		(inhd2->common.tag) ? " " : " not ");

    if (debug_level >= 1)
	fprintf(stderr, "%s: num_freqs = %d\n", PROG, num_freqs);

    tag_match = (inhd1->common.tag && inhd2->common.tag
		 && src_sf1 > 0.0
		 && src_sf2 > 0.9999*src_sf1
		 && src_sf2 < 1.0001*src_sf1);
    rec_freq_match = (record_freq1 > 0.0
		      && record_freq2 > 0.9999*record_freq1
		      && record_freq2 < 1.0001*record_freq1);

    REQUIRE(tag_match || rec_freq_match,
	    "input files are mismatched in record_freq "
	    "and (if tagged) in src_sf");

    /*
     * PROCESS OPTIONS AND PARAMETERS.
     */

    (void) read_params(param_name, SC_NOCOMMON, NULL);

    /* OUTPUT SELECTION */

    REQUIRE(!a_specified || !A_specified,
	    "both -a and -A specified");

    /* RECORD RANGE */

    startrec1 = 1;
    startrec2 = 1;
    nan = 0;

    if (r_arg1 != NULL) {

	long	endrec;

	endrec = LONG_MAX;
	lrange_switch(r_arg1, &startrec1, &endrec, NO);
	REQUIRE(endrec >= startrec1,
		"empty range of records specified");