getpdat2.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 Proccessing, Inc
	(C) 1985, Entropic Processing, Inc.
         */

#ifdef SCCS
static char *sccsid = "@(#)getpdat2.c	1.4 10/14/86 EPI";
#endif

#define bufsize 	90
#define	bufsize2	bufsize + bufsize
#include <sps/sps.h>
#include <stdio.h>
extern  FILE * hptr, *ifd;
extern int     wflt_order;
extern float lpf_num[], lpf_den[];
getpdata (rdata, pdata, datatype, nsamples)
float   rdata[], pdata[];
char datatype;
int     nsamples;
{
    int     i, j, lpf_order = 2;
    int     beg, lnt;
    float   t, tmpbuf[180], *ptr1, *ptr2;
    double  r[20];
    static int  bufptr = -bufsize, rel_locn = 0;
    static float    rawdata[2 * bufsize];
    static float    ready_pdata[bufsize],
                    ready_rdata[bufsize], rstate[2], pstate[2];
    static float    preemp_cons, preemp_state = 0.0;

    if (bufptr < 0)		/*  initialization */
    {
	for (i = 0; i < lpf_order; i++)
	    rstate[i] = pstate[i] = 0.0;
	get_fdata (ifd, datatype, rawdata, bufsize);
	bufptr = bufsize;
/*	rel_locn = smpnum + frmlen + bufsize/2; */
    }
    for (i = 0; i < nsamples; i++)
    {
	if (bufptr == bufsize)	/*  read new data */
	{
	    get_fdata (ifd, datatype, &rawdata[bufsize], bufsize);
	    rel_locn += bufsize;
	    get_auto (rawdata, bufsize2, r, 1);
	    preemp_cons = -r[1];
	    if (preemp_cons < -0.9375)
		preemp_cons = -0.9375;

	    ptr1 = &rawdata[0];
	    ptr2 = &tmpbuf[0];
	    for (j = 0; j < bufsize2; j++)
	    {
		*ptr2++ = *ptr1 + preemp_cons * preemp_state;
		preemp_state = *ptr1++;
	    }
	    preemp_state = rawdata[bufsize-1];

	    /*  update data memory */
	    ptr1 = &rawdata[0];
	    ptr2 = &rawdata[bufsize];
	    for (j = 0; j < bufsize; j++)
		*ptr1++ = *ptr2++;

	    /* perform lp filtering of raw speech */
	    ptr1 = &tmpbuf[bufsize / 2 + 1];
	    ptr2 = ready_rdata;
	    for (j = 0; j < bufsize; j++)
	    {
		t = *ptr1++ - lpf_den[1] * rstate[0] - lpf_den[2] * rstate[1];
		*ptr2++ = lpf_num[0] * t + lpf_num[1] * rstate[0] + lpf_num[2] * rstate[1];
		rstate[1] = rstate[0];
		rstate[0] = t;
	    }

	    get_frame_boundary (tmpbuf, &beg, &lnt);
	    if (hptr)
		fprintf (hptr, "%d\t%d\n", beg, lnt);

	    get_residue (&tmpbuf[beg], lnt);

	    /* perform lp filtering of residual signal */
	    ptr1 = &tmpbuf[bufsize / 2 + 1];
	    ptr2 = ready_pdata;
	    for (j = 0; j < bufsize; j++)
	    {
		t = *ptr1++ - lpf_den[1] * pstate[0] - lpf_den[2] * pstate[1];
		*ptr2++ = lpf_num[0] * t + lpf_num[1] * pstate[0] + lpf_num[2] * pstate[1];
		pstate[1] = pstate[0];
		pstate[0] = t;
	    }

	    bufptr = 0;
	}
	rdata[i] = ready_rdata[bufptr];
	pdata[i] = ready_pdata[bufptr++];
    }
}


/*  compute optimum frame boundary */
get_frame_boundary (data, ptr_beg, ptr_lnt)
float   data[];
int    *ptr_beg, *ptr_lnt;
{
    double  energy, min_energy;
    float  *ptr1, *ptr2;
    int     beg, end, i;

    energy = 0.0;
    min_energy = energy;
    beg = 0;
    ptr1 = &data[wflt_order];
    ptr2 = &data[0];
    for (i = 1; i < bufsize / 2 - wflt_order; i++, ptr1++, ptr2++)
    {
	energy += *ptr1 * *ptr1 - *ptr2 * *ptr2;
	if (energy < min_energy)
	{
	    beg = i;
	    min_energy = energy;
	}
    }

    ptr1 = &data[3 * bufsize / 2 + wflt_order];
    ptr2 = &data[3 * bufsize / 2];
    energy = 0.0;
    min_energy = energy;
    end = 0;
    for (i = 1; i < bufsize / 2 - wflt_order; i++, ptr1++, ptr2++)
    {
	energy += *ptr1 * *ptr1 - *ptr2 * *ptr2;
	if (energy < min_energy)
	{
	    end = i;
	    min_energy = energy;
	}
    }
    *ptr_beg = beg;
    *ptr_lnt = 3 * bufsize / 2 + wflt_order + end - beg;
}

/*	Compute residual signal using Burg method */
get_residue (f, lnt)
float   f[];
int     lnt;
{
    int     n, stage;
    double register num, den;
    float   b[180], fn, bn, gain, *fptr, *bptr;
#define rcstage num

    for (n = 0, fptr = f, bptr = b; n < lnt; n++)
	*bptr++ = *fptr++;

    gain = 1.0;
    for (stage = 1; stage <= wflt_order; stage++)
    {
	int wt, wt1 = 1, wt2 = wflt_order - stage + 1, wt3 = lnt - stage;
	num = den = 0.0;
	for (n = stage, fptr = &f[n], bptr = &b[n - stage]; n < lnt; n++)
	{
	    wt = wt2;
	    if ( wt > wt1) wt = wt1;
	    if ( wt > wt3) wt = wt3;
	    wt1++;
	    wt3--;

	    bn = *bptr++;
	    fn = *fptr++;
	    num += fn * bn * wt;
	    den += wt * (fn * fn + bn * bn);
	}

	rcstage = 2.0 * num / den;
	for (n = stage, fptr = &f[n], bptr = b; n < lnt; n++)
	{
	    bn = *bptr;
	    fn = *fptr;
	    *bptr++ = bn - rcstage * fn;
	    *fptr++ = fn - rcstage * bn;
	}
	gain *= (1.0 - rcstage * rcstage);
	if (gain <.02)
	{
	    if (hptr)
		fprintf (hptr,
			"WARNING: High LPC Gain(%d\t%g)\n", stage, gain);
	    return;
	}
    }
}