lpcana.c

speech signal process tools

    fea_oh = new_header (FT_FEA);
    fea_oh -> common.tag = YES;
    fea_oh -> variable.refer = input_sd_file;
    if (init_anafea_hd (fea_oh, (long)lpc_order, (long)
			lpc_order, (long)1, (long)1, (long)1, 
			    (short)0, (short)0) != 0)
    {
	Fprintf (stderr, "error filling FEA_ANA header");
	exit (1);
    }
    add_source_file (fea_oh, input_sd_file, sd_ih);
    (void)strcpy (fea_oh -> common.prog, "lpcana");
    (void)strcpy (fea_oh -> common.vers, Version);
    (void)strcpy (fea_oh -> common.progdate, Date);
    (void)add_comment (fea_oh, cmd_line);

    *(float *) get_genhd ("src_sf", fea_oh) = get_genhd_val("record_freq", sd_ih, (double) -1.);

    *(long*) get_genhd ("start", fea_oh) = start;
    *(long*) get_genhd ("nan", fea_oh) = nan;
    *add_genhd_s ("p_offset", (short *) NULL, 1, fea_oh) = P_OFFSET;
    (void)add_genhd_c ("dcrem", "yes", 0, fea_oh);
    (void)add_genhd_c ("psynch", "yes", 0, fea_oh);
    *add_genhd_s ("matsiz", (short *) NULL, 1, fea_oh) = matsiz; 
    *(short *) get_genhd ("spec_rep", fea_oh) = RC;
    *(long *) get_genhd("frmlen", fea_oh) = (long)0; 
    (void) add_genhd_s ("nominal_frame_size", &frame_size, 1, fea_oh);
    *add_genhd_e("method", NULL, 1, ana_methods, fea_oh) = anal_method;
    if (Fflag) 
      (void)add_genhd_c ("old_framing", "yes", 0, fea_oh);
    (void) update_waves_gen(sd_ih, fea_oh, (float)start, (float)0);
    if (sflag) 
	*add_genhd_s("sinc_n", NULL, 1, fea_oh) = sincn;

/* Write Output Header */
    write_header (fea_oh, fptr_fea_ana_file);

/* Allocate records for ESPS files */
    anafea_rec = allo_anafea_rec (fea_oh);

/* tcons = Exponential decay factor with time constant of 2.5 msecs */
    tcons = 400.0 / get_genhd_val("record_freq", sd_ih, (double) -1.);
    if(tcons < 0){
	Fprintf(stderr, "lpcana: record_freq < 0 in input file- exiting.\n");
	exit(1);
    }
    tcons = pow (0.25, tcons);

    if (start > 0)
	fea_skiprec (fptr_sd_file, (long) start - 1, sd_ih);

/*
 * MAIN LOOP
*/
    
    /* save for putsym in common */

    start_orig = start; 

    for (frame_num = 1; frame_num <= frame_count; frame_num++)
    {
	if (debug_level)
	    Fprintf (stderr, "Locn = %d\n", start);

/* Read Speech Data */
	if(get_sd_recf(&raw_data[frame_beg], 
	    frame_size, sd_ih, fptr_sd_file) == 0)
	{
	    Fprintf(stderr, "EOF hit while reading data\n");
	    exit(1);
	}

/* Remove DC in the speech */
	remove_dc (&raw_data[frame_beg], frame_size);

/* Compute Residual Signal */

	whiten_data (&raw_data[frame_beg - lpc_order], frame_size, lpc_order, matsiz, &res_data[frame_beg]);

/* Backward Scanning of residual data */
	j = frame_beg + frame_size - 1;
	bpeak = tmpbuf[frame_size - 1] = res_data[j--];
	for (dataptr = frame_size - 2; dataptr > -1; dataptr--)
	{
	    bpeak = tcons * bpeak;
	    t = res_data[j--];
	    if (t > bpeak)
		bpeak = tmpbuf[dataptr] = t;
	    else
		tmpbuf[dataptr] = 0;
	}

/* Forward Scanning of residual data */
	j = frame_beg - P_OFFSET;
	for (dataptr = 0; dataptr < frame_size; dataptr++)
	{
	    fpeak = tcons * fpeak;
	    if (tmpbuf[dataptr] > fpeak)
	    {
		fpeak = tmpbuf[dataptr];
		if (dataptr - plocn < min_pp || zcross_count < 2)
		{
		    if (peakval < fpeak)
		    {
			plocn = dataptr;
			peakval = fpeak;
			zcross_count = 0;
		    }
		}
		else
		{
		    if (plocn < 0)
		    {
			double  ratio = tcons;
			int     lnt = dataptr - plocn;

			ratio = pow (ratio, (double) lnt);
			if (fpeak * ratio < peakval)
			    build_frame (plocn);
		    }
		    else
			build_frame (plocn);

		    peakval = fpeak;
		    plocn = dataptr;
		    zcross_count = 0;
		}
	    }
	    if (raw_data[j] * raw_data[j - 1] < 0)
		zcross_count++;
	    j++;
	}

/* Shift Data in */
	j = frame_size;
	for (i = 0; i < 3 * frame_size; i++)
	{
	    raw_data[i] = raw_data[j];
	    res_data[i] = res_data[j++];
	}

	for (i = 0; i <= pulse_count; i++)
	    fplocn[i] -= frame_size;

	plocn -= frame_size;
	start += frame_size;
    }

/*
 * Write Common info, if appropriate
 */
    if(strcmp(fea_ana_file, "<stdout>") !=0){
	if(putsym_s("filename", input_sd_file) != 0)
	    Fprintf(stderr, "Could not write into ESPS Common");
	(void)putsym_s("prog", argv[0]);
	(void)putsym_i("start", (int)start_orig);
	(void)putsym_i("nan", (int)nan);
    }
    (void)fclose (fptr_sd_file);
    exit(0);
    return(0);
}

build_frame (next_pulse_locn)
short   next_pulse_locn;
{
    static short    flnt = 0;
    long    i, j, k, m, pulse_length;
    float   lpc_filter[TWENTY], rc[TWENTY], rcf[TWENTY], gain, fgain;
    static float *frame_data = NULL;
    short   vcount = 0;
    short   noz;
    float   se, t, peak, threshold;
    static float    prev_t = 0;
    double  r[TWENTY], rf[TWENTY];

    if (frame_data == NULL) 
	frame_data = (float *) malloc(buf_size * sizeof *frame_data);

    pulse_length = next_pulse_locn - fplocn[pulse_count];

    if (flnt + pulse_length / 2 > dynamic_frame_size)
    {
	dynamic_frame_size -= flnt - frame_size;
	if (debug_level)
	{
	    Fprintf (stderr, "Locn = %d\tSize = %d\n", fplocn[0] + start, flnt);
	    for (i = 0; i < pulse_count; i++)
		Fprintf (stderr, "%d\tLocn = %d\tSize = %d\n", i, fplocn[i], fpsize[i]);
	}
/* Spectrum Analysis of Data */

	j = fplocn[0] + frame_beg - lpc_order - P_OFFSET;

	(void)compute_rc(&raw_data[j], flnt + lpc_order, anal_method, 0, 
			 WT_RECT, lpc_order, sincn, 20 , 1e-5, rc, r, &gain); 

/* Make voicing decision */

	for (k = 0; k < pulse_count; k++)
	{
	    se = 0;
	    j = frame_beg + fplocn[k];
	    peak = res_data[j];
	    j -= P_OFFSET;
	    for (i = 0; i < fpsize[k]; i++)
		se += res_data[j++];
	    fpse[k] = se;

	    noz = 0;
	    j = frame_beg + fplocn[k] - P_OFFSET;
	    for (i = 0; i < fpsize[k]; i++)
	    {
		t = raw_data[j++];
		if (prev_t * t < 0)
		    noz++;
		prev_t = t;
	    }

	    fptype[k] = UNVOICED_FRM;
	    t = fpsize[k];
	    threshold = PEAK_THRESHOLD * (1.0 - pow (tcons, t));
	    if (noz <= ZCROSS_THRESHOLD * fpsize[k] || se < peak * threshold)
		fptype[k] = VOICED_FRM;

	    if (noz < 2 || noz < 0.03 * fpsize[k])
		fptype[k] = UNVOICED_FRM;

	    if (rc[1] < K1_THRESHOLD || gain > LPC_GAIN_THRESHOLD)
		fptype[k] = UNVOICED_FRM;

	    if (se < 0.4 * peak * threshold)
		fptype[k] = VOICED_FRM;

	    if (debug_level && peak && gain)
		Fprintf (stderr, "%6d\t%3d\t%6.3f\t%3d\t%8.2f\t%6.3f\t%6.3f\t%d\n",
			start + fplocn[k], fpsize[k], se / peak, noz, 1.0 / gain, rc[1], rc[2], fptype[k]);
	}

	vcount = 0;
	for (i = 0; i < pulse_count; i++)
	    if (fptype[i] == VOICED_FRM)
		vcount++;

/* Voiced frames should have all blocks voiced */
	if (2 * vcount > pulse_count)
	    for (i = 0; i < pulse_count; i++)
		fptype[i] = VOICED_FRM;

/* Output Analysis Data */

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

	    if (fpsize[i])
	    {
		*anafea_rec -> tag = start + fplocn[i];
		*anafea_rec -> frame_len = fpsize[i];
		if (fptype[i] == VOICED_FRM)
		{
		    *anafea_rec -> frame_type = VOICED;
		    *anafea_rec -> voiced_fraction = 1.0;
		}
		else
		{
		    *anafea_rec -> frame_type = UNVOICED;
		    *anafea_rec -> voiced_fraction = 0.0;
		}

		anafea_rec -> p_pulse_len[0] = fpsize[i];
		*anafea_rec -> num_pulses = 1;

		fpse[i] = fpse[i] / fpsize[i];

		if (fpse[i] == 0)
			fpse[i] = 1.0;

		if (!Fflag && (fpsize[i] > 12)) {
		    for (m = 0; m < fpsize[i]; m++) 
			frame_data[m] = raw_data[m + fplocn[i] + frame_beg];

		    if (debug_level > 1) 
			pr_farray(frame_data, fpsize[i], "DC-removed frame");

		    (void)compute_rc(frame_data, fpsize[i], anal_method, 0, 
				     WT_RECT, lpc_order, sincn, 20 , 1e-5, 
				     rcf, rf, &fgain); 

		    for (j = 0; j < lpc_order; j++)
			anafea_rec -> spec_param[j] = rcf[j + 1];
		}
		else {

		    if (debug_level && (!Fflag) && (fpsize[i] <= 12))
			Fprintf(stderr, 
		    "lpcana: very short pulse, didn't recalculate spectrum\n");


		    for (j = 0; j < lpc_order; j++)
			anafea_rec -> spec_param[j] = rc[j + 1];
		}
		anafea_rec -> raw_power[0] = fpse[i] / gain;
		anafea_rec -> lpc_power[0] = fpse[i];

		put_anafea_rec (anafea_rec, fea_oh, fptr_fea_ana_file);
	    }
	}
	fplocn[0] = fplocn[pulse_count];
	pulse_count = 0;
	flnt = 0;
    }
    fpsize[pulse_count] = pulse_length;
    flnt += pulse_length;
    pulse_count++;
    fplocn[pulse_count] = next_pulse_locn;
}

/* Remove dc value in a data sequence */
remove_dc (x, size)
    float   x[];
    short   size;
{
#define c0	63.0 / 64.0
#define c1	1.0 / 64.0
    short   i;
    static float    dc = 0.0;
    for (i = 0; i < size; i++)
    {
	dc = c0 * dc + c1 * x[i];
	x[i] -= dc;
    }
}

/* Whiten a given segment of data using linear predicition */
whiten_data (x, lnt, order, matsiz, output)
    float   x[], output[];
    int     order, lnt, matsiz;
{
/*
Whiten a given segment of data

input:
	x - input data array of size "lnt"

parameters:
	order - No. of reflection coefficients
	matsiz - sample covariance matrix size

output:
	output - square of the residual signal
	
*/
    register double num, den, dtemp;
    register float  stemp;
    double  f0Tf0, b0Tb0, f0Tb0;
    float  *ptr1, *f, *b;
    int     i, n, stage, vctsiz;
#define rcstage	num

    lnt += order;
/* Allocate memory for f and b arrays */
    f = (float *) malloc ((unsigned)((lnt+1) * sizeof *f));
    if (f == NULL)
    {
	Fprintf (stderr, "whitendata: couldn't allocate dynamic memory for array - f\n");
	exit (1);
    }

    b = (float *) malloc ((unsigned)((lnt+1) * sizeof *b));
    if (b == NULL)
    {
	Fprintf (stderr, "whitendata: couldn't allocate dynamic memory for array - b\n");
	exit (1);
    }

    for (n = 0; n < lnt; n++)
	f[n] = b[n] = x[n];

    f[lnt] = b[lnt] = 0.0;

    vctsiz = lnt - matsiz + 1;

    /* Computesignal energy */

    for (dtemp = 0.0, n = 0; n < vctsiz; n++)
	dtemp += f[n] * f[n];

    b0Tb0 = den = dtemp;
    f0Tf0 = den + f[vctsiz] * f[vctsiz] - f[0] * f[0];

    for (ptr1 = &f[vctsiz], n = 0; n < matsiz - 1; n++, ptr1++)
    {
	dtemp += *ptr1 * *ptr1 - f[n] * f[n];
	den += dtemp;
    }
    den = 2 * den - b0Tb0 - dtemp;

    for (stage = 1; stage <= order; stage++)
    {

	/* Compute first inner product */

	for (ptr1 = &f[stage], dtemp = 0.0, n = 0; n < vctsiz; n++)
	    dtemp += *ptr1++ * b[n];

	/*  Compute recursively the inner product of other vectors */

	f0Tb0 = num = dtemp;
	for (n = stage; n < matsiz - 1; n++)
	{
	    i = vctsiz + n;
	    dtemp += f[i] * b[i - stage] - f[n] * b[n - stage];
	    num += dtemp;
	}

	/*  Compute reflection coefficient and the residuals */

	rcstage = 2.0 * num / den;
	for (ptr1 = b, n = stage; n < lnt; n++)
	{
	    stemp = *ptr1;
	    *ptr1++ = stemp - rcstage * f[n];
	    f[n] -= rcstage * stemp;
	}

	/* Update b0Tb0 and compute of bnTbn */

	dtemp = b0Tb0 + rcstage * rcstage * f0Tf0 - 2.0 * rcstage * f0Tb0;
	f0Tf0 = f0Tf0 + rcstage * rcstage * b0Tb0 - 2.0 * rcstage * f0Tb0;
	i = matsiz - stage - 1;
	for (b0Tb0 = dtemp, n = 0; n < i; n++)
	{
	    stemp = b[n + vctsiz];
	    dtemp += stemp * stemp - b[n] * b[n];
	}

	/* Update den */

	den = den * (1.0 - rcstage * rcstage) - f0Tf0 - dtemp;

	/* Update f0Tf0 for next stage */

	stemp = f[vctsiz + stage];
	f0Tf0 += stemp * stemp - f[stage] * f[stage];
    }

/* Compute residual square and store in output array */
    n = 0;
    for (i = order; i < lnt; i++)
	output[n++] = f[i] * f[i];

    free ((char *) f);
    free ((char *) b);
}