bs_dist.c

来自「speech signal process tools」· C语言 代码 · 共 1,029 行 · 第 1/2 页

	default:
	    ERROR("type of symbol \"nan\" is not int");
	}
    }

    if (debug_level >= 1)
	fprintf(stderr,
		"%s: start = {%ld, %ld}, nan = %ld\n",
		PROG, startrec1, startrec2, nan);

    REQUIRE(startrec1 >= 1,
	    "range starts before beginning of first input file");
    REQUIRE(startrec2 >= 1,
	    "range starts before beginning of second input file");

    /* DISTORTION TYPE */

    distortion_type = BSD;
    if (m_arg != NULL) {

	distortion_type = lin_search(d_types, m_arg);

    } else if (symtype("distortion_type") != ST_UNDEF) {

	distypesym = getsym_s("distortion_type");

	REQUIRE(distypesym != NULL,
		"type of symbol \"distortion_type\" is not STRING");
	distortion_type = lin_search(d_types, distypesym);

    } else {
	distortion_type = BSD;
    }

    REQUIRE(distortion_type == BSD,
	    "distortion_type is not BSD");

    if (debug_level >= 1)
	fprintf(stderr,
		"%s: distortion_type = \"%s\"\n",
		PROG, d_types[distortion_type]);

    /* THRESHOLD */

    check_thresh = NO;	/* this becomes YES only if all conditions for
			 * checking input frame powers against the
			 * threshold are met:  we are computing BSD, the
			 * threshold is greater than 0, and the tot_power
			 * field is defined in both input files
			 */

    if (distortion_type == BSD) {

	if (e_arg != NULL) {

	    threshold = atof(e_arg);

	} else if (symtype("threshold") != ST_UNDEF) {

	    REQUIRE(symtype("threshold") == ST_FLOAT,
		    "type of symbol \"threshold\" is not FLOAT");
	    threshold = getsym_d("threshold");

	} else {
	    threshold = 0.0;
	}

	REQUIRE(threshold >= 0.0,
		"threshold is negative");

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

	if (threshold > 0.0) {

	    REQUIRE(tot_power_defined1,
		    "non-zero threshold specified, "
		    "but no tot_power field defined in first input file");
	    REQUIRE(tot_power_defined2,
		    "non-zero threshold specified, "
		    "but no tot_power field defined in second input file");

	    check_thresh = YES;
	}

	if (debug_level >= 1)
	    fprintf(stderr,
		    "%s: total power of input frames "
		    "will%sbe compared with threshold\n",
		    PROG, (check_thresh ? " " : " not "));
    }

    /* PERCEPTUAL WEIGHTING FACTORS */

    if (symtype("perceptual_weights") != ST_UNDEF) {

	REQUIRE(symtype("perceptual_weights") == ST_FARRAY
		&& symsize("perceptual_weights") >= num_freqs,
		"parameter \"perceptual_weights\" is too short "
		"or not a FLOAT array");

	TRYALLOC(float, symsize("perceptual_weights"), p_wts, "p_wts");
	TRYALLOC(float, symsize("perceptual_weights"), log_wts, "log_wts");
	getsym_fa("perceptual_weights", p_wts, num_freqs);

	for (i = 0; i < num_freqs; i++) {
	    REQUIRE(p_wts[i] > 0.0,
		    "non-positive value in perceptual_weights array");
	    log_wts[i] = 10 * log10(p_wts[i]);
	}

    } else {

	TRYALLOC(float, num_freqs, p_wts, "p_wts");
	TRYALLOC(float, num_freqs, log_wts, "log_wts");

	for (i = 0; i < num_freqs; i++)
	{
	    double  a = 2.6;	/* Numerator coefficient ... */
	    double  b = 1.6;	/* ... and denominator coefficient of bilinear
				 * preemphasis filter.  These magic numbers are
				 * from section VI.B of Wang, Sekey, & Gersho,
				 * (1992).
				 */
	    double  w0 = (a + 1.0)*(a + 1.0)/((b + 1.0)*(b + 1.0));
	    double  two_x = 2.0 * cos(PI*freqs1[i]/4000.0);
	    double  num = (a + two_x) * a + 1.0;
	    double  den = (b + two_x) * b + 1.0;
	    /* 
	     * numerator = |a+z|^2
	     * denominator = |b+z|^2
	     * where z = exp (2i pi f) and f = freqs1[i]/8000
	     */
	    
	    p_wts[i] = (num/den)/w0;
	    log_wts[i] = 10.0 * log10(p_wts[i]);
	}
    }

    if (debug_level >= 2) {

	fprintf(stderr,
		"%s: log perceptual weighting factors (dB) ...",
		PROG);

	for (i = 0; i < num_freqs; i++) {
	    if (i%5 == 0)
		fprintf(stderr, "\n [%2ld]\t", i);
	    fprintf(stderr, " %13.6g", log_wts[i]);
	}

	fprintf(stderr, "\n");
    }

    /*
     * CREATE OUTPUT FILE
     */

    if (!a_specified) {

	/* CREATE HEADER */

	if (debug_level >= 1)
	    fprintf(stderr, "%s: creating output header\n", PROG);

	outhd = new_header(FT_FEA);
	add_source_file(outhd, inname1, inhd1);
	add_source_file(outhd, inname2, inhd2);
	add_comment(outhd, get_cmd_line(argc, argv));
	(void) strcpy(outhd->common.prog, PROG);
	(void) strcpy(outhd->common.vers, VERSION);
	(void) strcpy(outhd->common.progdate, DATE);
	outhd->variable.refer = inhd1->variable.refer;

	if (tag_match) {

	    outhd->common.tag = YES;
	    *add_genhd_d("src_sf", NULL, 1, outhd) = src_sf1;

	    if (debug_level >= 1)
		fprintf(stderr,
			"%s: output file is tagged; src_sf = %g\n",
			PROG, src_sf1);

	} else {

	    outhd->common.tag = NO;

	    if (debug_level >= 1)
		fprintf(stderr, "%s: output file is not tagged\n",
			PROG);
	}

	if (rec_freq_match) {

	    update_waves_gen(inhd1, outhd, (float) startrec1, 1.0f);

	    if (debug_level >= 1)
		fprintf(stderr, "%s: output record_freq defined; "
			"value is %g\n",
			PROG, record_freq1);

	} else {

	    if (debug_level >= 1)
		fprintf(stderr, "%s: output record_freq undefined\n",
			PROG);
	}

	start = add_genhd_l("start", NULL, 2, outhd);
	start[0] = startrec1;
	start[1] = startrec2;
	*add_genhd_l("nan", NULL, 1, outhd) = nan;
	*add_genhd_e("distortion_type",
		     NULL, 1, d_types, outhd) = distortion_type;

	if (distortion_type == BSD) {
	    *add_genhd_d("threshold", NULL, 1, outhd) = threshold;
	}

	(void) add_genhd_f("perceptual_weights", p_wts, num_freqs, outhd);

	outfld = (distortion_type == BSD) ? "BSD" : "MBSD";

	add_fea_fld(outfld, 1, 0, NULL, FLOAT, NULL, outhd);

	/* OPEN FILE */

	if (debug_level >= 1)
	    fprintf(stderr, "%s: opening output file\n", PROG);

	outname = eopen(PROG, outname, "w", NONE, NONE, NULL, &outfile);

	if (debug_level >= 1)
	    fprintf(stderr, "%s: output file = \"%s\"\n", PROG, outname);

	/* WRITE HEADER */

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

	write_header(outhd, outfile);

    } /* end if (!a_specified) */

   /*
    * SET UP I/O RECORD STRUCTS AND DATA BUFFERS
    */

    if (debug_level >= 1)
	fprintf(stderr,
		"%s: allocating input and output record structrures\n",
		PROG);

    inrec1 = allo_feaspec_rec(inhd1, FLOAT);
    inbuf1 = inrec1->re_spec_val;
    intag1 = inrec1->tag;

    inrec2 = allo_feaspec_rec(inhd2, FLOAT);
    inbuf2 = inrec2->re_spec_val;
    intag2 = inrec2->tag;

    if (!a_specified) {

	outrec = allo_fea_rec(outhd);
	outdata = (float *) get_fea_ptr(outrec, outfld, outhd);
	outtag = &outrec->tag;
    }

    if (check_thresh) {

	tot_power1 = inrec1->tot_power;
	tot_power2 = inrec2->tot_power;
    }

    /*
     * MAIN READ-WRITE LOOP
     */

    fea_skiprec(infile1, startrec1 - 1, inhd1);
    fea_skiprec(infile2, startrec2 - 1, inhd2);

    if (a_specified || A_specified) {

	/* set up to compute overall average distortion */

	sum_distortion = 0.0;
	norm = 0.0;
    }

    if (debug_level >= 1)
	 fprintf(stderr, "%s: reading and processing data\n", PROG);

    /* If nan == 0, the loop terminates, via a break, when the end of
     * either input file is reached.
     * If nan != 0, the loop terminates after nan records have been
     * read and processed, unless the end of an input file is reached
     * sooner.  In the latter case, a "premature end of file" warning
     * is issued. 
     */

    for (n = 0; nan == 0 || n < nan; n++) {

	/* READ INPUT FRAMES */

	if (get_feaspec_rec(inrec1, inhd1, infile1) == EOF) {

	    if (nan != 0)
		fprintf(stderr,
			"%s: WARNING: premature end of file in \"%s\"; "
			"only %d records read\n",
			PROG, inname1, n);
	    break;
	}

	if (get_feaspec_rec(inrec2, inhd2, infile2) == EOF) {

	    if (nan != 0)
		fprintf(stderr,
			"%s: WARNING: premature end of file in \"%s\"; "
			"only %d records read\n",
			PROG, inname2, n);

	    break;
	}

	if (tag_match && inhd1->common.tag && *intag1 != *intag2) {

	    fprintf(stderr, "%s: WARNING: "
		    "tags mismatch in record %ld. (%ld != %ld)\n",
		    PROG, n + 1, *intag1, *intag2);

	    tag_match = NO;
	}

	if (debug_level >= 3)
	    fprintf(stderr, "%s: processing record %ld\n",
		    PROG, n + 1);

	if (debug_level >= 4) {
	    print_farray(1, n + 1, "input data", num_freqs, inbuf1);
	    print_farray(2, n + 1, "input data", num_freqs, inbuf2);
	}

	/* CONVERT TO DB */

	switch (spec_type1) {

	case SPTYP_PWR:
	    for (j = 0; j < num_freqs; j++)
		inbuf1[j] = 10.0 * log10(inbuf1[j]);
	    break;
	case SPTYP_DB:
	    /* nothing to do */
	    break;
	default:
	    /* Should be impossible; we've already checked spec_type1 */
	    ERROR("ERROR determining spec_type of first input file");
	}

	switch (spec_type2) {

	case SPTYP_PWR:
	    for (j = 0; j < num_freqs; j++)
		inbuf2[j] = 10.0 * log10(inbuf2[j]);
	    break;
	case SPTYP_DB:
	    /* nothing to do */
	    break;
	default:
	    /* Should be impossible; we've already checked spec_type2 */
	    ERROR("ERROR determining spec_type of first input file");
	}

	if (debug_level >= 4) {
	    print_farray(1, n + 1, "input data (DB)", num_freqs, inbuf1);
	    print_farray(2, n + 1, "input data (DB)", num_freqs, inbuf2);
	}

	/* APPLY PERCEPTUAL WEIGHTS */

	for (j = 0; j < num_freqs; j++) {

	    inbuf1[j] += log_wts[j];
	    inbuf2[j] += log_wts[j];
	}

	if (debug_level >= 4) {
	    print_farray(1, n + 1,
			 "input data (phones)", num_freqs, inbuf1);
	    print_farray(2, n + 1,
			 "input data (phones)", num_freqs, inbuf2);
	}

	/* PHONES TO SONES */

	for (j = 0; j < num_freqs; j++) {

	    inbuf1[j] = exp((inbuf1[j] - 40) * LOG2_OVER_10);
	    inbuf2[j] = exp((inbuf2[j] - 40) * LOG2_OVER_10);
	}

	if (debug_level >= 4) {
	    print_farray(1, n + 1,
			 "input data (sones)", num_freqs, inbuf1);
	    print_farray(2, n + 1,
			 "input data (sones)", num_freqs, inbuf2);
	}

	/* SQUARED EUCLIDEAN DISTANCE */

	frame_distortion = 0.0;

	for (j = 0; j < num_freqs; j++) {

	    double  diff = inbuf1[j] - inbuf2[j];

	    frame_distortion += diff * diff;
	}

	if (debug_level >= 3)
	    fprintf(stderr, "%s: frame %ld distortion = %g\n",
		    PROG, n + 1, frame_distortion);

	/* OVERALL AVERAGE DISTORTION */

	if (a_specified || A_specified)
	{
	    if (!check_thresh
		|| (*tot_power1 >= threshold && *tot_power2 >= threshold))
	    {
		sum_distortion += frame_distortion;
	    }
	    else /* (check_thresh
		  *  && (*tot_power1 < thresh || *tot_power2 < thresh) */
	    {
		if (debug_level >= 3)
		    fprintf(stderr,
			    "%s: excluding frame %ld distortion---"
			    "power below threshold\n",
			    PROG, n + 1);
	    }

	    for (j = 0; j < num_freqs; j++) {
		norm += inbuf1[j] * inbuf1[j];
	    }

	    if (debug_level >= 3)
		fprintf(stderr,
			"%s: frame %ld cumulative distortion = %g, "
			"norm factor %g\n",
			PROG, n + 1, sum_distortion, norm);
	}

	/* WRITE FRAME OUTPUT */

	if (!a_specified) {

	    *outdata = frame_distortion;

	    if (outhd->common.tag) {

		*outtag = *intag1;
	    }

	    if (debug_level >= 3)
		fprintf(stderr, "%s: writing output record %ld\n",
			PROG, n + 1);

	    put_fea_rec(outrec, outhd, outfile);
	}

    } /* end for (n = 0; ... ) */
    
    /*
     * WRAP UP
     */

    if (debug_level >= 1)
	fprintf(stderr, "%s: processed %ld frames\n", PROG, n);

    if (a_specified || A_specified) {

	printf("%g\n", sum_distortion/norm);
    }

    return 0;	/* equivalent to exit(0); */
}


/*
 * formatted printout of float array for debugging
 */

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

    fprintf(stderr, "%s: file %d record %ld %s ...",
	    PROG, filenum, recnum, text);

    for (j = 0; j < numelem; j++) {

	if (j%5 == 0)
	    fprintf(stderr, "\n [%ld]\t", j);
	fprintf(stderr, " %13.6g", data[j]);
    }

    fprintf(stderr, "\n");
}