lpc_lib.c

MELPe 1200 bps, fixed point

    Shortword	p_cos;
    Shortword	temp1, temp2;
    Longword	L_temp1, L_temp2;


	if (firstTime){
		/* for (i = 0; i < DFTLENGTH; i++)
			 lsp_cos[i] = cos(i*(TWOPI / DFTLENGTH)); */

		/* cos_fxp() takes Q15 input.  (TWO/DFTLENGTH) above is DFTLENGTH_D4  */
		/* in Q15.  The first for loop fills lsp_cos[] in the first quadrant, */
		/* and the next loop fills lsp_cos[] for the other three quadrants.   */

		temp1 = 0;
		for (i = 0; i <= DFTLENGTH_D4; i++){
			lsp_cos[i] = cos_fxp(temp1);
			lsp_cos[i + DFTLENGTH_D2] = negate(lsp_cos[i]);
			temp1 = add(temp1, DFTLENGTH_D4);
		}
		temp1 = DFTLENGTH_D4;
		temp2 = DFTLENGTH_D4;
		for (i = 0; i < DFTLENGTH_D4; i++){
			lsp_cos[temp1] = negate(lsp_cos[temp2]);
			lsp_cos[temp1 + DFTLENGTH_D2] = lsp_cos[temp2];
			temp1 = add(temp1, 1);
			temp2 = sub(temp2, 1);
		}

		/* compute default values for freq[] */
		/* (1./p2) in Q15 */
		default_w = divide_s(ONE_Q11, shl(order, 10));
		/* freq[0] = (0.5/p2) in Q15 */
		default_w0 = shr(default_w, 1);

		firstTime = FALSE;
	}

	prev_less = TRUE;
	L_fill(mag, 0x7fffffff, 2);
	fill(s_mag, 0x7fff, 2);

	/*	p2 = p/2; */
	p2 = shr(order, 1);
	count = 0;

	/*	Search all frequencies for minima of Pc(w) */
	for (i = 0; i <= DFTLENGTH_D2; i++){
		p_cos = i;
		/* mag2 = 0.5 * lsp[p2]; */
		L_temp2 = L_mult(lsp[p2], X05_Q14);                  /* mag[2] in Q25 */
		for (j = sub(p2, 1); j >= 0; j--){
			L_temp1 = L_shr(L_mult(lsp[j], lsp_cos[p_cos]), 1);
			L_temp2 = L_add(L_temp2, L_temp1);
			p_cos = add(p_cos, i);
			if (p_cos > DFTLENGTH - 1)
				p_cos -= DFTLENGTH;
		}
		s_mag[2] = extract_h(L_temp2);
		mag[2] = L_abs(L_temp2);

		if (mag[2] < mag[1]){
			prev_less = TRUE;
		} else {
			if (prev_less){
				if ((s_mag[0] ^ s_mag[2]) < 0){
					/* Minimum frequency found */
					/*	freq[count] = i - 1 + (0.5 *
							 (mag[0] - mag[2]) / (mag[0] + mag[2] - 2*mag[1]));
							 freq[count] *= (2. / DFTLENGTH) ;*/
					L_temp1 = L_shr(L_sub(mag[0], mag[2]), 1);
					L_temp2 = L_sub(mag[0], L_shl(mag[1], 1));
					L_temp2 = L_add(L_temp2, mag[2]);
					temp1 = L_divider2(L_temp1, L_temp2, 0, 0);
                                                              /* temp1 in Q15 */
					temp1 = shr(temp1, 9);                              /* Q6 */
					temp2 = sub(i, 1);
					temp1 = add(shl(temp2, 6), temp1);
					freq[count] = divide_s(temp1, shl(DFTLENGTH, 5));
					count = add(count, 1);
				}
			}
			prev_less = FALSE;
		}
		L_v_equ(mag, &(mag[1]), 2);
		v_equ(s_mag, &(s_mag[1]), 2);
	}

	/* Verify that all roots were found.  Under normal conditions the condi-  */
	/* tion for the following IF statement should never be true.              */

	if (count != p2){
		/* use default values */
		freq[0] = default_w0;
		for (i = 1; i < p2; i++)
			freq[i] = add(freq[i - 1], default_w);
	}
}


/* LPC_PRED2REFL                                                              */
/*	  get refl coeffs from the predictor coeffs                               */
/*	  Input:                                                                  */
/*		 lpc- the predictor coefficients                                      */
/*		 order- the predictor order                                           */
/*	  Output:                                                                 */
/*		 refc- the reflection coefficients                                    */
/*	  Returns:                                                                */
/*		 energy - energy of residual signal                                   */
/* Reference:  Markel and Gray, Linear Prediction of Speech                   */
/*                                                                            */
/* Q values:                                                                  */
/* lpc[] - Q12, *refc - Q15,                                                  */

Shortword lpc_pred2refl(Shortword lpc[], Shortword *refc, Shortword order)
{
	register Shortword	i, j;
	Longword	acc;
	Shortword	*b, *b1, e;
	Shortword	energy = ONE_Q15;
	Shortword	shift, shift1, sign;
	Shortword	temp;


	b = v_get(order);
	b1 = v_get((Shortword) (order - 1));

	/* equate temporary variables (b = lpc) */
	v_equ(b, lpc, order);

	/* compute reflection coefficients */
	for (i = sub(order, 1); i >= 0; i--){

		if( b[i] >= 4096 )
			b[i] = 4095;
		if( b[i] <= -4096 )
			b[i] = -4095;

		acc = L_mult(b[i], b[i]);
		acc = L_sub(ONE_Q25, acc);
		acc = L_shl(acc, 6);					/* Q31 */
		energy = mult(energy, extract_h(acc));

		shift = norm_l(acc);
		e = extract_h(L_shl(acc, shift));

		v_equ(b1, b, i);

		for (j = 0; j < i; j++){
			/*	b[j] = (b1[j] - local_refc*b1[i - j])/e; */
			acc = L_mult(b[i], b1[i-j-1]);						/* Q25 */
			acc = L_sub(L_shl(L_deposit_l(b1[j]), 13), acc);	/* Q25 */

			/* check signs of temp and e before division */
			sign = extract_h(acc);
			acc = L_abs(acc);
			shift1 = norm_l(acc);	
			temp = extract_h(L_shl(acc, shift1));
			if( temp > e ){
				temp = shr(temp, 1);
				shift1 = sub(shift1, 1);
			}
			b[j] = divide_s(temp, e);
			shift1 = sub(shift1, 3);
			shift1 = sub(shift1, shift);
			b[j] = shr(b[j], shift1);					/* b[j] in Q12 */

			if ( sign < 0)
				b[j] = negate(b[j]);
		}
	}

	*refc = shl(b[0], 3);
	v_free(b1);
	v_free(b);
	return(energy);
}

/* LPC_LSP2PRED                                                               */
/*	  get predictor coefficients from the LSPs                                */
/* Synopsis: lpc_lsp2pred(w,a,p)                                              */
/*	  Input:                                                                  */
/*		 lsf- the LSPs                                                        */
/*		 order- the predictor order                                           */
/*	  Output:                                                                 */
/*		 lpc- the predictor coefficients                                      */
/* Reference:  Kabal and Ramachandran                                         */
/*                                                                            */
/* Q values:                                                                  */
/* lsf - Q15, lpc - Q12, c - Q14                                              */

Shortword lpc_lsp2pred(Shortword lsf[], Shortword lpc[], Shortword order)
{
	register Shortword	i, j, k;
	Shortword	p2;
	Shortword	c0, c1;
	Longword	*f0, *f1;
	Longword	L_temp;


	/* ensure minimum separation and sort */
	lpc_clmp(lsf, 0, order);

	/* p2 = p/2 */
	p2 = shr(order, 1);

	f0 = L_v_get((Shortword) (p2 + 1));
	f1 = L_v_get((Shortword) (p2 + 1));

	/* f is Q25 */
	f0[0] = f1[0] = (Longword) ONE_Q25;

	/* -2.0*cos((double) lsf[0]*M_PI) */
	f0[1] = L_shr(L_deposit_h(negate(cos_fxp(lsf[0]))), 5);
	f1[1] = L_shr(L_deposit_h(negate(cos_fxp(lsf[1]))), 5);

	k = 2;

	for (i = 2; i <= p2; i++){
		/* c is Q14 */
		/* multiply by 2 is considered as Q15->Q14 */
		c0 = negate(cos_fxp(lsf[k]));
		k ++;
		c1 = negate(cos_fxp(lsf[k]));
		k ++;

		f0[i] = f0[i - 2];
		f1[i] = f1[i - 2];

		for (j = i; j >= 2; j--){
			/* f0[j] += c0*f0[j - 1] + f0[j - 2] */
			L_temp = L_mpy_ls(f0[j - 1], c0);
			L_temp = L_add(L_shl(L_temp, 1), f0[j - 2]);
			f0[j] = L_add(f0[j], L_temp);

			/* f1[j] += c1*f1[j - 1] + f1[j - 2] */
			L_temp = L_mpy_ls(f1[j - 1], c1);
			L_temp = L_add(L_shl(L_temp, 1), f1[j - 2]);
			f1[j] = L_add(f1[j], L_temp);
		}

		f0[1] = L_add(f0[1], L_shl(L_mpy_ls(f0[0], c0), 1));
		f1[1] = L_add(f1[1], L_shl(L_mpy_ls(f1[0], c1), 1));
	}

	for (i = sub(p2, 1); i >= 0; i--){
		/* short f (f is a Q14) */
		f0[i + 1] = L_add(f0[i + 1], f0[i]);
		f1[i + 1] = L_sub(f1[i + 1], f1[i]);

		/* lpc[] is Q12 */
		/* lpc[i] = 0.50*(f0[i] + f1[i]) */
		/* lpc[p + 1 - i] = 0.50*(f0[i] - f1[i]) */
		/* Q25 -> Q27 -> Q12 */
		lpc[i] = extract_h(L_shl(L_add(f0[i + 1], f1[i + 1]), 2));
		lpc[order - 1 - i] = extract_h(L_shl(L_sub(f0[i + 1], f1[i + 1]), 2));
	}

	v_free(f0);
	v_free(f1);
	return(0);
}


/* Name: lpc_syn- LPC synthesis filter.                                       */
/*	Aliases: lpc_synthesis                                                    */
/*	Description:                                                              */
/*		LPC all-pole synthesis filter                                         */
/*                                                                            */
/* 		for j = 0 to n-1                                                      */
/* 			y[j] = x[j] - sum(k=1 to p) y[j-k] a[k]                           */
/*                                                                            */
/*	Inputs:                                                                   */
/*		x- input vector (n samples, x[0..n-1])                                */
/*		a- lpc vector (order p, a[1..p])                                      */
/*		order- order of lpc filter                                            */
/*		length- number of elements in vector which is to be filtered          */
/*		y[-p..-1]- filter memory (past outputs)                               */
/*	Outputs:                                                                  */
/*		y- output vector (n samples, y[0..n-1]) (Q0)                          */
/*	Returns: NULL                                                             */
/*	Includes:                                                                 */
/*		lpc.h                                                                 */
/*                                                                            */
/*	Systems and Info. Science Lab                                             */
/*	Copyright (c) 1995 by Texas Instruments, Inc.  All rights reserved.       */

Shortword lpc_syn(Shortword x[], Shortword y[], Shortword a[], Shortword order,
				  Shortword length)
{
	register Shortword	i, j;
	Longword	accum;


/* Tung-chiang believes a[] is Q12, x[] and y[] are Q0. */

	for (j = 0; j < length; j++){
		accum = L_shr(L_deposit_h(x[j]), 3);
		for (i = order; i > 0; i--)
			accum = L_msu(accum, y[j - i], a[i - 1]);
		/* Round off output */
		accum = L_shl(accum, 3);
		y[j] = round(accum);
	}
	return(0);
}