ref_cof.c

speech signal process tools

	Fprintf (stderr,"ref_cof: start-1=%ld nframes=%d fsize=%ld order=%d\n", 
start, nframes, fsize, order);
	}
    
/*Filtering initialization.  */

    if (pre == YES)
	{

	nsiz = getsym_i ("nsiz");
	dsiz = getsym_i ("dsiz");

	if (nsiz > dsiz) siz = nsiz;
	else siz = dsiz;

/* CHANGE: xin array is of size da_max. */
	xin = (float*) calloc ((unsigned)da_max, sizeof(float));

	state = (float*) calloc ((unsigned)siz, sizeof(float));
	fzeros = (float*) calloc ((unsigned)nsiz, sizeof(float));
	zeros = (double*) calloc ((unsigned)nsiz, sizeof(double));

/* Initialize the coefficient arrays. */

        (void) sprintf (na_name, "%s_num", filter_name);
        if (getsym_da(na_name,zeros,nsiz) != nsiz)
            {
            Fprintf (stderr, "Wrong number of numerator coefficients in %s\n",
 param_file);
            exit (-1);
            }

	if (dsiz > 0)
	    {
	    fpoles = (float*) calloc ((unsigned)dsiz, sizeof(float));
	    poles = (double*) calloc ((unsigned)dsiz, sizeof(double));

            (void) sprintf(da_name,"%s_den", filter_name);
            if (getsym_da(da_name,poles,dsiz) != dsiz)
                {
	        Fprintf (stderr, "Wrong number of denominator coefficients in %s\n", param_file);
	        exit (-1);
	        }
            for (i=0; i<dsiz; i++)
	        {
	        fpoles[i] = poles[i];
	        }
	    }

        for (i=0; i<nsiz; i++)
	    {
	    fzeros[i] = zeros[i];
	    }

/* Store the coefficients in a zfunc. */

        pzfunc = new_zfunc (nsiz, dsiz, fzeros, fpoles);
	if (pzfunc == NULL) Fprintf(stderr,"pzfunc is NULL\n");

/* Store the order and coefficients in the output header. */

	(void)add_genzfunc("pre_emp", pzfunc, oh);
	}

    if (debug_level)
	{
	Fprintf(stderr, "ref_cof: zfunc stored in output header:\n");
	Fprintf(stderr,"ref_cof: nsiz=%d dsiz=%d\n", nsiz, dsiz);
	    for (i=0; i<pzfunc->nsiz; i++)
	        Fprintf (stderr,"zeros[%d]=%f\n", i, pzfunc->zeros[i]);
	    for (i=0; i<pzfunc->dsiz; i++)
	        Fprintf (stderr,"poles[%d]=%f\n", i, pzfunc->poles[i]);
	}

    if (pre == NO ) {
        if (start > 0) (void) skiprec (fpin, start, size_rec(ih));
    }
    else {

/* Initialize the data arrays from the data if it is an FIR filter. */

        if (dsiz == 0) {
            start_p = start - (long) siz;
            if (start_p > 0) skiprec(fpin, start_p, size_rec (ih));
            for (i=0; i<siz; i++) {
	        if (start_p + i < 0) state[i] = 0;
	        else {
		    more = get_sd(&state[i], 1, ih, fpin);
		    if (!more) 
			ERROR_EXIT("not enough data for initialization");
		}
	    }
	}
/* For IIR filter, initialize the arrays to zero. */

        else {
            skiprec(fpin, start, size_rec(ih)); 
	    for (i=0; i<siz; i++) state[i] = 0;
	}
    }
/* CHANGE: x array is of size da_max. */
    x = (float*) calloc ((unsigned)da_max,sizeof(float));
    r = (double*) calloc ((unsigned)order+1,sizeof(double));

/* CHANGE: dly and dc are double arrays.  */
    dly = (double*) calloc ((unsigned)order+1,sizeof(double));
    dc = (double*) calloc ((unsigned)nframes, sizeof(double));

    rc = (float*) calloc ((unsigned)order+1,sizeof(float));

    oh->common.ndrec = nframes;

    *(float *) get_genhd("src_sf", oh) = ih->hd.sd->sf;
    *(short *) get_genhd("spec_rep", oh) = RC;
/* uncomment after FEA_ANA update
    *(long *) get_genhd("frmlen", oh) = fsize;
    *(long *) get_genhd("start", oh) = start + 1;
    *(long *) get_genhd("nan", oh) = (long) nframes*fsize;
*/

/*
 * Check analysis method and write generic header item
*/
    if((method = lin_search(analysis_methods, analysis)) == -1)
	    {
	    fprintf (stderr, "ref_cof: analisis method is %s\n", analysis);
	    ERROR_EXIT("Invalid spectrum analysis method");
	    }
    switch (method) {
	case BURG:
	case COV:
	case MBURG:
	case VBURG:
	    ERROR_EXIT("Analysis method is not supported yet");
	    break;
	case STRCOV:
	    if (fsize > DA_MAX) { 
		Fprintf(stderr, 
		  "ref_cof: STRCOV method not supported for frames > %d points\n", DA_MAX);
		exit(1);
	    }
	    if (order > 20) Fprintf(stderr,  
		"ref_cof: WARNING -- order > 20 with STRCOV takes long time");
	    break;
	case AUTOC:
	    break;
	default:
	    ERROR_EXIT("Invalid spectrum analysis method");
    }
 
    *(short *)add_genhd_e("method", NULL, 1, analysis_methods, oh) = method;
    *(short *)add_genhd_s("order", NULL, 1, oh) = order;

/* Check sinc_flag stuff */
    
    if (sinc_flag == YES) {
	*add_genhd_e("sinc_flag", NULL, 1, local_no_yes, oh) = YES;
	(void) sprintf (comment,"The autocorrelation function has been multiplied\n by sinc (t/%d).\n", sinc);	
	if (add_comment(oh, comment) == 0) Fprintf(stderr, 
	  "ref_cof: WARNING -- not enough space for comment\n");
    }
    else
	*add_genhd_e("sinc_flag", NULL, 1, local_no_yes, oh) = NO;

/* add DC_REMOVED info */

    switch(method) {
	case AUTOC:
	    add_genhd_c("DC_removed", "YES", (int)NULL, oh);
	    break;
	case STRCOV:
	    add_genhd_c("DC_removed", "NO", (int)NULL, oh);
	    break;
    }

/* Write out header */
    write_header (oh, fpout);

/* CHANGE: 8 lines deleted here.  */

/* Allocate space to store the analysis records. */

    recp = allo_anafea_rec (oh);
    *recp->frame_len = fsize;

/* CHANGE: Open up the temporary file. */
    TRYOPEN ("ref_cof", mktemp(temp_file), "w+", fptemp);
    if (debug_level)
	fprintf (stderr,"Temporary file is %s\n", temp_file);

/* CHANGE: This whole section to compute the d.c. offset
	   is new.    */
/* Compute the D.C. offset.   */

    nget = da_max;
    nleft = nframes*fsize;
    knt = 0;
    for (k=0; k < nframes; k++)
	{
	dc [k] = 0;
	for (j=0; j<fsize; j++)
	    {
	    if (knt == 0 && nleft > 0)
		{
		if (nleft < da_max) nget = nleft;
		else nget = da_max;
		if (pre == YES)
		    {
		    get_sd_recf (xin, (short)nget, ih, fpin);
		    block_filter (nget, xin, x, pzfunc, state);
		    }
		else
		    {
		    get_sd_recf (x, (short)nget, ih, fpin);
		    }
		nleft -= nget;
		}
	    dc [k] += x[knt];
	    knt = (knt + 1) % nget;
	    if (knt == 0)
		{
		put_sd_recf (x, nget, ih, fptemp);
		}
	    }
	dc [k] = dc [k]/fsize;
	}
    rewind (fptemp);

/*reset file pointer to beginning of data*/
    rewind (fpin);
    if((ih = read_header(fpin)) == NULL)
	NOTSPS("ref_cof", in_file);
    (void)skiprec(fpin, start, size_rec(ih));

    if (debug_level)
	{
	fprintf (stderr, "ref_cof: D.C. offsets:\n");
	for (k=0; k<nframes; k++)
	    fprintf (stderr, "ref_cof: dc[%d] = %g\n", k, dc[k]);
	}

/****************************************************************/

/* Begin the spectral analysis.  */

    nget = da_max;
    nleft = nframes*fsize;
    knt = 0;
    for (k=0; k<nframes && more; k++)
	{

/* Zero the arrays.  */

        for (i=0; i<= order; i++)
    	    {
	    r[i] = 0.0;
	    dly[i] = 0.0;
	    }
/* CHANGE: two lines deleted here. */
/* compute spectrum */ 

/* CHANGE: this section changed extensively because of temp file. */
	 switch (method) {
	    case AUTOC: 
		for (j=0; j<fsize; j++) 
		    {
		    if (knt == 0 && nleft > 0)
			{
			if (nleft < da_max) nget = nleft;
			else nget = da_max;
			get_sd_recf (x, (short)nget, ih, fptemp);
			nleft -= nget;
			}
		    xn = x[knt] - dc[k];
		    knt = (knt + 1) % nget;
		    for (i=order; i>0; i--) dly[i] = dly[i-1];
		    dly[0] = xn;
		    for (i=0; i<= order; i++) r[i] += xn*dly[i];
		    }
		for (i=1; i<= order; i++) r[i] = r[i]/r[0];
		r[0] = r[0]/fsize;
		break;
	    case STRCOV:
		/* This only works because fsize is limited to DA_MAX
		   for STRCOV method.   */
		more = get_sd (x, (int)fsize, ih, fpin);
		bestauto(x, (int)fsize, r, order, order+1, 0);
		break;

	    default: /*should never get here*/
		ERROR_EXIT("unsupported analysis method");
	 }
	if (sinc_flag == YES) 
	    for (i=1; i<=order; i++) {
		arg = (double) PI*i/sinc;
		r[i] = r[i]*(sin (arg)/(arg));
	    }
	get_atal (r, order, lpcfilter, rc, &gain);
	for (i=0; i<order; i++) recp->spec_param[i] = rc[i+1];
	*recp->tag = (long) (start + k*fsize);
	recp->raw_power[0] = r[0];
	res_energy = r[0]*gain;
	recp->lpc_power[0] = res_energy;
	recp->p_pulse_len[0] = 0;
	*recp->frame_len = fsize;
	if (debug_level == 1) {
	    for (i=0; i<= order; i++) 
		Fprintf (stderr, 
		    "r[%d]=%f rc[%d]=%f lpcfilter[%d]=%f\n", 
		    i,r[i],i,recp->spec_param[i],i,lpcfilter[i]);
		Fprintf (stderr,"\n\n");
		}
	put_anafea_rec (recp, oh, fpout);
	}
    fclose (fpin);
    fclose (fpout);

/* CHANGE: must unlink temp file.  */
    (void) unlink (temp_file);
}


int
get_sd(data, samp, hd, file)
float *data;
int samp;
struct header *hd;
FILE *file;
{
    if (get_sd_recf(data, samp, hd, file) == EOF)
	return 0;
    else
	return 1;
}