getrmsval.c

speech signal process tools

/*
| This material contains proprietary software of Entropic Processing, Inc.
| Any reproduction, distribution, or publication without the the prior
| written permission of Entropic Processing, Inc. is strictly prohibited.
| Any public distribution of copies of this work authorized in writing by
| Entropic Processing, Inc. must bear the notice
|
|              "Copyright 1986 Entropic Processing, Inc."
|
| Written by:  S. Shankar Narayan
|
| Module:	getrmsval.c
| 
*/

#ifdef SCCS
static char *sccs_id = "@(#)getrmsval.c	1.3 5/12/86 EPI";
#endif
#include <stdio.h>
float   get_rmsval (rc, delta_rc, order, se, state, dempstate, blksze, locn)
float   rc[], se, state[], dempstate[], delta_rc[];
int     order, locn, blksze;
{
    float   gain, c1[100], c2[100], tmp[100], t1, t2, gain2;
    float   s_state[100], e_state[100], ts, te, st_energy, ex_energy, cr_energy;
    float   const, root1, root2, rmsval, lattice_filt();
	extern FILE *hptr;
    double  sqrt ();
    extern float  num[3];

    int     i, j, n;

    for (i = 0; i < order; i++)
	tmp[i+1] = rc[i];
    rctoc (tmp, order, c1, &gain);
/* Compute signal power from the power of pre-emphasized signal */

    n = order + 1;
    /* Compute composite filter coefficients */
    c1[0] = -1.0;
    c1[n] = 0.0;
    for (i = 1; i <= n; i++)
	c2[i] = c1[i] + c1[i - 1] * num[1];
    /* Compute gain of composite filter */
    gain2 = 1.0;
    for (i = n; i > 0; i--)
    {
	t1 = c2[i];
	t2 = 1.0 - t1 * t1;
	gain2 *= t2;
	for (j = 1; j <= i - 1; j++)
	    tmp[j] = (c2[j] + t1 * c2[i - j]) / t2;
	for (j = 1; j <= i - 1; j++)
	    c2[j] = tmp[j];
    }
    /* Signal power to be matched */
    se = se * gain / gain2;

/* Compute energies in the inverse filter response due to synthesizer
		states and impulse */
    for (i = 0; i <= order; i++)
    {
	e_state[i] = 0.0;
	s_state[i] = state[i];
    }
    ts = dempstate[0];
    te = 0.0;

    const = 1.0;
    ts = -ts * num[1] + lattice_filt (rc, delta_rc, order, s_state, 0.0);
    te = -te * num[1] + lattice_filt (rc, delta_rc, order, e_state, const);
    st_energy = ts * ts;
    ex_energy = te * te;
    cr_energy = ts * te;
    const = -1.0 / (blksze - 1);
    for (i = 1; i < blksze; i++)
    {
	ts = -num[1] * ts + lattice_filt (rc, delta_rc, order, s_state, 0.0);
	te = -num[1] * te + lattice_filt (rc, delta_rc, order, e_state, const);
	st_energy += ts * ts;
	ex_energy += te * te;
	cr_energy += ts * te;
    }

/* Solve: A^2 * ex_energy + 2 * A * cr_energy + st_energy = se */

    const = cr_energy * cr_energy
	+ (se - st_energy) * ex_energy;

    if (const < 0)
    {
/* imaginary roots */
	if (hptr) 
		fprintf (hptr,
		 "locn = %d\tInput and output energies not matched\n", locn);
	rmsval = se * gain;
	rmsval = sqrt ((double) rmsval);
    }
    else
    {
	const = sqrt ((double) const);
	root1 = (-cr_energy + const) / ex_energy;
	root2 = (-cr_energy - const) / ex_energy;
	rmsval = root1;
	if (root2 > 0 && root2 < root1)
	    rmsval = root2;

	if (rmsval < 0)
	{
	    const = sqrt (se * gain);
	if (hptr) 
		fprintf (hptr,
		    "OOPs Negative Pulse at locn = %d: Pcomputed =%9.2f\tPused =%9.2f\n",
		    locn, rmsval, const);
	    rmsval = const;
	}
    }
    return (rmsval);
}