mbs_dist.c

speech signal process tools

	phontoson(PY, SY);
	phontoson(PN, SN);

	if (debug_level >= 4) {
	    print_darray(1, n + 1,
			 "spread spectrum (sone)", BSIZE - 3, SX);
	    print_darray(2, n + 1,
			 "spread spectrum (sone)", BSIZE - 3, SY);
	    print_darray(0, 0,
			 "noise spectrum (sone)", BSIZE - 3, SN);
	}

	/* COMPUTE FRAME DISTORTION */

	frame_distortion = measure(SX, SY, SN);

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

	/* OVERALL AVERAGE DISTORTION */

	if (included) {

	    num_proc++;
	    sum_distortion += frame_distortion;
	}

	/* WRITE FRAME OUTPUT */

	if (!a_specified) {

	    *outdata = frame_distortion;
	    *outflag = included;

	    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: sum of %ld distortions out of %ld processed is %g.\n",
		PROG, num_proc, n, sum_distortion);

    if (a_specified || A_specified) {

	printf("%g\n", (num_proc == 0) ? 0 : sum_distortion/num_proc);
    }

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


/*
 * Convert FEA_SPEC data to spec_type PWR
 */

void
spec_to_pwr(int     num,
	    float   *src,
	    float   *src_im,
	    int     src_type,
	    double  *dest)
{
    int     j;

    switch (src_type)
    {
    case SPTYP_PWR:
	for (j = 0; j < num; j++)
	    dest[j] = src[j];
	break;
    case SPTYP_DB:
	for (j = 0; j < num; j++)
	    dest[j] = pow(10.0, 0.1*src[j]);
	break;
    case SPTYP_REAL:
	for (j = 0; j < num; j++)
	    dest[j] = src[j]*src[j];
	break;
    case SPTYP_CPLX:
	for (j = 0; j < num; j++)
	    dest[j] = src[j]*src[j] + src_im[j]*src_im[j];
	break;
    default:
	ERROR("unrecognized spec_type");
    }
}

/*
 * Integrate power spectral density over frame
 */

double
frame_pwr(int       num,
	  double    *pwr_spec,
	  double    df)
{
    int     j;
    double  sum;

    sum = 0.5*(pwr_spec[0] + pwr_spec[num-1]);

    for (j = 1; j < num - 1; j++)
	sum += pwr_spec[j];

    return sum*df;
}

/* Initialize an array of equally spaced frequencies */

double *
init_freqs(int      num,
	   double   fmax)
{
    int     j;
    double  *freqs;
    double  den;

    TRYALLOC(double, num, freqs, "frequencies");

    den = (double) (num - 1);
    for (j = 0; j < num; j++)
	freqs[j] = (j / den) * fmax;

    return freqs;
}

/*
 * Initialize spread function.
 */

double *
init_spread_func()
{
    static double spread_func[2*BSIZE-1];
    double  *S = &spread_func[BSIZE-1];
    int     k;
    double  x;

    for (k = 1 - BSIZE; k < BSIZE; k++) {

	x = (double) k + 0.474;
	x = 15.81 + 7.5*x - 17.5*sqrt(1.0 + x*x);
	S[k] = pow(10.0, 0.1*x);
    }

    return S;
}


/*
 * Initialize absolute hearing threshold in dB by the formula of Terhardt.
 *
 *  thrshld(f) = { 3.64(f/1000)^(-0.8) - 6.5exp[-0.6(f/1000 - 3.3)^2]
 *                 + 0.001(f/1000)^4 }
 *
 * Reference :  Terhardt, E., Stoll, G. and Seewann, M, "Algorithm for
 *        extraction of pitch and pitch salience from complex tonal
 *        signals", J. Acoust. Soc. Am., vol. 71(3), Mar., 1982.
 */

double *
init_thresh(int num, double *freqs)
{
    static double ABS_TH[BSIZE];
    int     i, j, k;
    double  f;
    double  *L;
    double  fsq, t, SUM;

    TRYALLOC(double, num, L, "working storage for \"init_thresh\"");

    L[0] = 0.0;

    for (i = 1; i < num; i++) {

        f = freqs[i]/1000.0;
        fsq = f * f;
        t = f - 3.3;
        L[i] =
	    3.64 / pow(f, 0.8)
		- 6.5 * exp(-0.6 * t * t)
		    + 0.001 * fsq * fsq;
    }

    for (i = 1; i <= BSIZE; i++) {

        SUM = 0.0;
        k = 0;

        for (j = 1; j < num; j++) {

            if (BARK[i-1] <= freqs[j] && freqs[j] < BARK[i]) {

                SUM += L[j];
                k++;
            }
        }

        if (k == 0)
            ABS_TH[i-1] = 0.0;
        else
            ABS_TH[i-1] = SUM / k;
    }

    return ABS_TH;
}


/*
 * Normalize a power spectrum to an overall average power of
 * 375000.0, assuming that "mean_pwr" gives the value before
 * normalization, and "pwr_spec" indicates the current frame.
 * The original "mbsd" normalizes the _unwindowed_ sampled
 * data to an overall RMS of 1000.0.  This is as close as we
 * can come in terms of the mean square of the hanning-windowed
 * signal.  (The mean square value of the hanning window function
 * is 0.375).
 */

void
normalize(int       num,
	  double    *pwr_spec,
	  double    mean_pwr,
	  double    factor)
{
    int     j;
    double  ratio;

    if (mean_pwr == 0)
	return;

    ratio = factor*375000.0/mean_pwr;

    for (j = 0; j < num; j++)
	pwr_spec[j] *= ratio;

    pwr_spec[0] *= 0.5;
}


/*
 * Decide whether to include current frame in computation of overall
 * average distortion.
 */

int
check_frame(double  pwr1,
	    double  pwr2,
	    double  XTHRESH,
	    double  YTHRESH)
{
    return (pwr1 > XTHRESH && pwr2 > YTHRESH);
}


/*
 * Compute critical-band spectrum from power spectrum.
 */

void
bk_frq(int	num,
       double	*freqs,
       double	*pwr_spec,
       double	*bk_spec)
{
    int     i, j;

    for (i = 0; i < BSIZE; i++) {

	bk_spec[i] = 0.0;

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

	    if (BARK[i] <= freqs[j] && freqs[j] < BARK[i+1])
		bk_spec[i] += pwr_spec[j];
	}
    }
}


/*
 * Apply spreading function to critical-band spectrum
 * to compute spread critical-band spectrum.
 */

void
spread(double *bk_spec, double *spread_spec)
{
    int     i, j;

    for (i = 0; i < BSIZE; i++) {

	spread_spec[i] = 0;

	for (j = 0; j < BSIZE; j++)
	    spread_spec[i] += S[i-j] * bk_spec[j];
    }
}


/*
 * Convert spread Bark spectrum to phons.
 */

void
dbtophon(double *spread_spec, double *phon_spec)
{
    int     i, j;
    double  t1;
    double  Ti;

    for (i = 0; i < BSIZE-3; i++) {

	Ti = 10.0 * log10(spread_spec[i+3]);

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

	    if (Ti < eqlcon[j][i])
		break;
	}

	if (j == 0) {

	    phon_spec[i] = phons[0];

	} else if (j == PHSIZE) {

	    phon_spec[i] = phons[PHSIZE-1];

	} else {

	    t1 = (Ti - eqlcon[j-1][i]) / (eqlcon[j][i] - eqlcon[j-1][i]);
	    phon_spec[i] = phons[j-1] + t1 * (phons[j] - phons[j-1]);
	}
    }
}


/*
 * Convert spectrum in phons to sones.
 */

void
phontoson(double *phon_spec, double *sone_spec)
{
    int     i;

    for (i = 0; i < BSIZE-3; i++) {

        if (phon_spec[i] >= 40.0) {

	    sone_spec[i] = pow(2.0, 0.1*(phon_spec[i] - 40.0));
	}
	else {

	    sone_spec[i] = pow(phon_spec[i] / 40.0, 2.642);
	}
    }
}


/*
 * Compute frame distortion.
 */

double
measure(double *orig, double *dist, double *noise)
{
    int     i;
    double  d;
    double  t;

    d = 0.0;

    for (i = 0; i < BSIZE - 3; i++) {

	t = fabs(orig[i] - dist[i]) - noise[i];

	if (t > 0.0)
	    d += t;
    }

    return d;
}


/*
 * Spectral flatness measure (based on the ratio of the arithmetic
 * and geometric means).
 */

double
sfm(int num, double *pwr_spec)
{
    double  alpha;
    double  a_mean;		/* arithmetic mean */
    double  g_mean;		/* geometric mean */
    int     i;

    a_mean = 0;
    g_mean = 0;

    for (i = 0; i < num; i++) {

	a_mean += pwr_spec[i];
	g_mean += log(pwr_spec[i]);
    }

    a_mean = a_mean / (double) num;
    g_mean = exp(g_mean / (double) num);

    alpha = 10.0 * log10(g_mean / a_mean) / (-60.0);

    if (alpha > 1.0)
	alpha = 1.0;

    return alpha;
}


/*
 * Noise masking threshold.
 */

void
thresh2(double alpha, double *spread_spec, double *noise_thr)
{
    int     i;
    double  t;

    for (i = 0; i < BSIZE; i++) {

        t = alpha * (14.5 + (double) i + 1.0) + (1.0 - alpha) * 5.5;
	t = 10.0 * log10(spread_spec[i]) - t;

	if (t < ABS_TH[i])
	    t = ABS_TH[i];

	noise_thr[i] = pow(10.0, t/10.0);
    }
}


/*
 * Formatted printout of float array for debugging.
 */

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

    if (filenum == 0)
	fprintf(stderr, "%s: %s ...", PROG, text);
    else
	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");
}


/*
 * Formatted printout of double array for debugging
 */

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

    if (filenum == 0)
	fprintf(stderr, "%s: %s ...", PROG, text);
    else
	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");
}