mbs_dist.c

speech signal process tools

	if (endrec != LONG_MAX)
	    nan = endrec - startrec1 + 1;

	if (r_arg2 != NULL) {

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

	    if (endrec != LONG_MAX) {

		if (nan == 0)
		    nan = endrec - startrec2 + 1;
		else
		    REQUIRE(nan == endrec - startrec2 + 1,
			    "inconsistent range sizes specified");
	    }

	} else {
	    startrec2 = startrec1;
	}

    } else {

	switch (symtype("start")) {

	case ST_INT:
	    startrec1 = startrec2 = getsym_i("start");
	    break;
	case ST_IARRAY:
	    REQUIRE(getsym_ia("start", startrecs, 2) == 2,
		    "size of array parameter \"start\" is not 2");
	    startrec1 = startrecs[0];
	    startrec2 = startrecs[1];
	    break;
	case ST_UNDEF:
	    break;
	default:
	    ERROR("type of symbol \"start\" is not int or int array");
	}

	switch (symtype("nan")) {

	case ST_INT:
	    nan = getsym_i("nan");
	    break;
	case ST_UNDEF:
	    break;
	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");

    /*
     * 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) = MBSD;

	distfld = "MBSD";
	add_fea_fld(distfld, 1, 0, NULL, FLOAT, NULL, outhd);
	flagfld = "MBSD_included";
	add_fea_fld(flagfld, 1, 0, NULL, CODED, noyes, 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;
    inbufim1 = inrec1->im_spec_val;
    intag1 = inrec1->tag;

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

    if (!a_specified) {

	outrec = allo_fea_rec(outhd);
	outdata = (float *) get_fea_ptr(outrec, distfld, outhd);
	outflag = (short *) get_fea_ptr(outrec, flagfld, outhd);
	outtag = &outrec->tag;
    }

    TRYALLOC(double, num_freqs, PSX, "first power spectrum");
    TRYALLOC(double, num_freqs, PSY, "second power spectrum");

    /*
     * FIRST PASS
     */

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

    if (lseek(fileno(infile1), 0L, SEEK_CUR) < 0) {

	if (debug_level >= 1)
	    fprintf(stderr, "%s: first input is a pipe; making temp file\n",
		    PROG);

	tmpfile1 = tmpfile();
	tmphd1 = copy_header(inhd1);	/* get version without Esignal
					 * header or the like; else we may
					 * get garbage by writing temp
					 * file as FEA and reading back as
					 * another format*/
	write_header(tmphd1, tmpfile1);
    }

    if (lseek(fileno(infile2), 0L, SEEK_CUR) < 0) {

	if (debug_level >= 1)
	    fprintf(stderr, "%s: second input is a pipe; making temp file\n",
		    PROG);

	tmpfile2 = tmpfile();
	tmphd2 = copy_header(inhd2);	/* get version without Esignal
					 * header or the like; else we may
					 * get garbage by writing temp
					 * file as FEA and reading back as
					 * another format*/
	write_header(tmphd2, tmpfile2);
    }

    max_xpower = 0.0;
    max_ypower = 0.0;
    mean_xpower = 0.0;
    mean_ypower = 0.0;

    if (debug_level >= 1)
	 fprintf(stderr, "%s: reading data---first pass\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: read record %ld\n",
		    PROG, n + 1);

	if (tmpfile1 != NULL) {

	    if (debug_level >= 3)
		fprintf(stderr, "%s: writing record to temp file 1\n",
			PROG);

	    put_fea_rec(inrec1->fea_rec, tmphd1, tmpfile1);
	}

	if (tmpfile2 != NULL) {

	    if (debug_level >= 3)
		fprintf(stderr, "%s: writing record to temp file 2\n",
			PROG);

	    put_fea_rec(inrec2->fea_rec, tmphd2, tmpfile2);
	}

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

	spec_to_pwr(num_freqs, inbuf1, inbufim1, spec_type1, PSX);
	spec_to_pwr(num_freqs, inbuf2, inbufim2, spec_type2, PSY);

	if (debug_level >= 4) {
	    print_darray(1, n + 1, "input data (PWR)", num_freqs, PSX);
	    print_darray(2, n + 1, "input data (PWR)", num_freqs, PSY);
	}

	tot_pwr1 = frame_pwr(num_freqs, PSX, delta_f);
	tot_pwr2 = frame_pwr(num_freqs, PSY, delta_f);

	if (debug_level >= 3)
	    fprintf(stderr, "%s: tot_pwr = {%g, %g}\n",
		    PROG, tot_pwr1, tot_pwr2);

	mean_xpower += tot_pwr1;
	mean_ypower += tot_pwr2;

	if (max_xpower < tot_pwr1)
	    max_xpower = tot_pwr1;
	if (max_ypower < tot_pwr2)
	    max_ypower = tot_pwr2;
    }

    if (n != 0) {

	mean_xpower /= n;
	mean_ypower /= n;
    }


    if (debug_level >= 1)
	 fprintf(stderr,
		 "%s: end of first pass; read %ld records\n"
		 "\tmean_xpower = %g, mean_ypower = %g\n"
		 "\tmax_xpower = %g, max_ypower = %g\n",
		 PROG, n,
		 mean_xpower, mean_ypower,
		 max_xpower, max_ypower);

    nan = n;

    if (tmpfile1 != NULL) {

	fclose(infile1);
	infile1 = tmpfile1;
    }

    if (tmpfile2 != NULL) {

	fclose(infile2);
	infile2 = tmpfile2;
    }

    rewind(infile1);
    rewind(infile2);

    /* The factor 375000.0/mean_{x,y}power is the normalizing factor that
     * is applied by the function "normalize".  The factor FRAME converts
     * power to energy for compatibility with program "mbsd".
     */
    XTHRESHOLD =
	max_xpower * FRAME
	    * (375000.0 / mean_xpower) * pow(10.0, -1.5); /* 15 dB below */
    YTHRESHOLD =
	max_ypower * FRAME
	    * (375000.0 / mean_ypower) * pow(10.0, -3.5); /* 35 dB below */
    
    if (debug_level >= 1)
	fprintf(stderr, "%s: XTHRESHOLD = %g, YTHRESHOLD = %g\n",
		PROG, XTHRESHOLD, YTHRESHOLD);

    FREQ = init_freqs(num_freqs, fmax);

    if (debug_level >= 2)
	print_darray(0, 0, "freqs", num_freqs, FREQ);
    
    S = init_spread_func();

    if (debug_level >= 2)
	print_darray(0, 0, "spread function", 2*BSIZE - 1, &S[1 - BSIZE]);

    ABS_TH = init_thresh(num_freqs - 1, FREQ);

    if (debug_level >= 2)
	print_darray(0, 0,
		     "absolute hearing thresh", BSIZE, ABS_TH);

    /*
     * MAIN READ-WRITE LOOP
     */

    inhd1 = read_header(infile1);
    inhd2 = read_header(infile2);

    if (tmpfile1 == NULL) {

	fea_skiprec(infile1, startrec1 - 1, inhd1);
    }

    if (tmpfile2 == NULL) {

	fea_skiprec(infile2, startrec2 - 1, inhd2);
    }

    /* set up to compute overall average distortion */

    sum_distortion = 0.0;
    num_proc = 0;

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

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

	/* READ INPUT FRAMES */

	REQUIRE(get_feaspec_rec(inrec1, inhd1, infile1) != EOF,
		"premature end of file in rereading first input file");
	REQUIRE(get_feaspec_rec(inrec2, inhd2, infile2) != EOF,
		"premature end of file in rereading second input file");

	if (debug_level >= 3)
	    fprintf(stderr, "%s: read 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 PWR */

	spec_to_pwr(num_freqs, inbuf1, inbufim1, spec_type1, PSX);
	spec_to_pwr(num_freqs, inbuf2, inbufim2, spec_type2, PSY);

	if (debug_level >= 4) {
	    print_darray(1, n + 1, "input data (PWR)", num_freqs, PSX);
	    print_darray(2, n + 1, "input data (PWR)", num_freqs, PSY);
	}

	/* NORMALIZE */

	normalize(num_freqs, PSX, mean_xpower, FRAME*delta_f);
	normalize(num_freqs, PSY, mean_ypower, FRAME*delta_f);

	if (debug_level >= 4) {
	    print_darray(1, n + 1,
			 "normalized power spectrum", num_freqs, PSX);
	    print_darray(2, n + 1,
			 "normalized power spectrum", num_freqs, PSY);
	}

	/* CHECK WHETHER FRAME IS TO BE INCLUDED */

	tot_pwr1 = frame_pwr(num_freqs, PSX, 1.0);
	tot_pwr2 = frame_pwr(num_freqs, PSY, 1.0);
	included = check_frame(tot_pwr1, tot_pwr2, XTHRESHOLD, YTHRESHOLD);

	if (debug_level >= 3)
	    fprintf(stderr,
		    "%s: tot_pwr = {%g, %g}\n"
		    "\tframe%sincluded in overall distortion\n",
		    PROG, tot_pwr1, tot_pwr2,
		    (included) ? " " : " NOT ");

	/* BARK SPECTRUM */

	bk_frq(num_freqs, FREQ, PSX, BX);
	bk_frq(num_freqs, FREQ, PSY, BY);

	if (debug_level >= 4) {
	    print_darray(1, n + 1, "critical-band spectrum", BSIZE, BX);
	    print_darray(2, n + 1, "critical-band spectrum", BSIZE, BY);
	}

	/* SPREAD BARK SPECTRUM */

	spread(BX, CX);
	spread(BY, CY);

	if (debug_level >= 4) {
	    print_darray(1, n + 1,
			 "spread critical-band spectrum", BSIZE, CX);
	    print_darray(2, n + 1,
			 "spread critical-band spectrum", BSIZE, CY);
	}

	/* CONVERT TO PHON */

	dbtophon(CX, PX);
	dbtophon(CY, PY);

	if (debug_level >= 4) {
	    print_darray(1, n + 1,
			 "spread spectrum (phon)", BSIZE - 3, PX);
	    print_darray(2, n + 1,
			 "spread spectrum (phon)", BSIZE - 3, PY);
	}

	/* NOISE MASKING THRESHOLD */

	alpha = sfm(num_freqs, PSX);

	if (debug_level >= 3)
	    fprintf(stderr, "%s: spectral flatness %g\n", PROG, alpha);

	thresh2(alpha, CX, CNMT);

	if (debug_level >= 4)
	    print_darray(0, n + 1, "noise masking threshold", BSIZE, CNMT);

	dbtophon(CNMT, PN);

	if (debug_level >= 4)
	    print_darray(0, n + 1,
			 "noise masking threshold (phon)", BSIZE - 3, PN);

	/* CONVERT TO SONE */

	phontoson(PX, SX);