lpc_lib.c

It is source code for Melp2.4kps vocoder using dsp tms320vc55x of ti

/*

*2.4 kbps MELP Proposed Federal Standard speech coder
*
*TMS320C5x assembly code
*
*version 1.0
*
*Copyright (c) 1998, Texas Instruments, Inc.  
*
*Texas Instruments has intellectual property rights on the MELP
*algorithm.  The Texas Instruments contact for licensing issues for
*commercial and non-government use is William Gordon, Director,
*Government Contracts, Texas Instruments Incorporated, Semiconductor
*Group (phone [phone]).

*/

/*************************************************************************
*
* The following code was hand optimized for the Texas Instuments
* TMS320C5x DSP by DSPCon, Inc.  For information, please contact DSPCon
* at:
* 
*                       DSPCon, Inc.
*                       380 Foothill Road
*                       Bridgewater, New Jersey 08807
*                       [phone]
*                       [e-mail]
*                       www.dspcon.com
*
*************************************************************************/



/*

  lpc_lib.c: LPC subroutines.

*/

/*  compiler include files  */

#include        <stdlib.h>
#include        <math.h>
#include "spbstd.h"
#include "mathhalf.h"
#include "mathdp31.h"
#include "math_lib.h"
#include "mat.h"
#include "lpc.h"



extern Shortword frames;

#define LPC_ORD         10
#define DFTLENGTH       512
#define DFTLENGTH_D2    (DFTLENGTH/2)
#define DFTLENGTH_D4    (DFTLENGTH/4)
#define ONE_Q11 (1<<11)
#define ONE_Q15   (Shortword)(((Longword)1<<15)-1)
#define ALMOST_ONE_Q14   ((1<<14)-2)

static Shortword default_w, default_w0;
       Shortword lsp_cos[DFTLENGTH * 3];    /* cosine table */

static Shortword cos_tab[129] =
    {
	32767, 32766, 32758, 32746, 32729, 32706, 32679, 32647, 32610,
	32568, 32522, 32470, 32413, 32352, 32286, 32214, 32138, 32058,
	31972, 31881, 31786, 31686, 31581, 31471, 31357, 31238, 31114,
	30986, 30853, 30715, 30572, 30425, 30274, 30118, 29957, 29792,
	29622, 29448, 29269, 29086, 28899, 28707, 28511, 28311, 28106,
	27897, 27684, 27467, 27246, 27020, 26791, 26557, 26320, 26078,
	25833, 25583, 25330, 25073, 24812, 24548, 24279, 24008, 23732,
	23453, 23170, 22884, 22595, 22302, 22006, 21706, 21403, 21097,
	20788, 20475, 20160, 19841, 19520, 19195, 18868, 18538, 18205,
	17869, 17531, 17190, 16846, 16500, 16151, 15800, 15447, 15091,
	14733, 14373, 14010, 13646, 13279, 12910, 12540, 12167, 11793,
	11417, 11039, 10660, 10279,  9896,  9512,  9127,  8740,  8351,
	7962,   7571,  7180,  6787,  6393,  5998,  5602,  5205,  4808,
	4410,   4011,  3612,  3212,  2811,  2411,  2009,  1608,  1206,
	804,     402,     0
    };


    static  Shortword   inputw[200];



#undef LOW_LIMIT




/*
    Name: lpc_clmp- Sort and ensure minimum separation in LSPs.
    Aliases: lpc_clamp
    Description:
	Ensure that all LSPs are ordered and separated
	by at least delta.  The algorithm isn't guarenteed
	to work, so it prints an error message when it fails
	to sort the LSPs properly.
    Inputs:
	w- lsp vector (order p, w[1..p])
	delta- the clamping factor
	p- order of lpc filter
    Outputs:
	w- the sorted and clamped lsps
    Returns: NULL
    See_Also:
    Includes:
	spbstd.h
	lpc.h
    Bugs:
	Currently only supports 10 loops, which is too
	complex and perhaps unneccesary.

    Systems and Info. Science Lab
    Copyright (c) 1995 by Texas Instruments, Inc.  All rights reserved.

*/
/*  Q values:
    w - Q15
    delta - Q15 */

#define MAX_LOOPS 10

Shortword lpc_clmp(Shortword *w, Shortword delta, Shortword p)
{
    Shortword i,j,unsorted;
    Shortword temp,d,step1,step2;

    /* sort the LSPs for 10 loops */
    for (j=0,unsorted=TRUE; unsorted && (j < MAX_LOOPS); j++)
    {
	for(i=1,unsorted=FALSE; i < p; i++)
	    if (w[i] > w[i+1])
	    {
		temp = w[i+1];
		w[i+1] = w[i];
		w[i] = temp;
		unsorted = TRUE;
	    }
    }

    /* ensure minimum separation */
    if (!unsorted)
    {
	for(j=0; j < MAX_LOOPS; j++)
	{
	    for(i=1; i < p; i++)
	    {
		/* increment complexity for if statement */
                
                d = sub(w[i+1],w[i]);     
		if (d < delta)
		{
                        
		    step1 = step2 = shr(sub(delta,d),1);

		    /* increment complexity for if statement */
                        
		    if (i==1 && (w[i] < delta)) {
                      step1 = shr(w[i],1);    
		    }
		    else {
		      /* increment complexity for if statement */
                      
		      if (i > 1)
		      {
			/* increment complexity for if statement */
                        
			temp = sub(w[i],w[i-1]);
			if (temp < delta) {
                          step1 = 0;    
			}
			else {
			  /* increment complexity for if statement */
                              
			  if (temp < shl(delta,1)) {
			    step1 = shr(sub(temp,delta),1);
			  }
			}
		      }
		    }

		    /* increment complexity for if statement */
                        
		    if (i==(p-1) && (w[i+1] > sub(ONE_Q15,delta))) {
                        step2 = shr(sub(ONE_Q15,w[i+1]),1);    
		    }
		    else {
		      /* increment complexity for if statement */
                      
		      if (i < (p-1))
		      {
			/* increment complexity for if statement */
                        
			temp = sub(w[i+2],w[i+1]);
			if (temp < delta) {
                          step2 = 0;    
			}
			else {
			  /* increment complexity for if statement */
                          
			  if (temp < shl(delta,1)) {
			    step2 = shr(sub(temp,delta),1);
			  }
			}
		      }
		    }
                    w[i] = sub(w[i],step1);    
                    w[i+1] = add(w[i+1],step2);    

		}
	    }
	}
    }

    return((Shortword)0);
} /* LPC_CLMP */


/*
    Name: lpc_schr- Schur recursion (autocorrelations to refl coef)
    Aliases: lpc_schur
    Description:
	Compute reflection coefficients from autocorrelations
	based on schur recursion.  Will also compute predictor
	parameters by calling lpc_refl2pred(3l) if necessary.
    Inputs:
	r- autocorrelation vector (r[0..p]).
	p- order of lpc filter.
    Outputs:
	a-     predictor parameters    (can be NULL)
	k_tmp- reflection coefficients (can be NULL)
    Returns:
	alphap- the minimum residual energy
    Includes:
	spbstd.h
	lpc.h
    See_Also:
	lpc_refl2pred(3l) in lpc.h or lpc(3l)

    Q values:
    r - Q0
    a - Q12
    k_tmp - Q15
*/

Shortword lpc_schr(Shortword *r, Shortword *a, Shortword *k_tmp, Shortword p)
{
    Shortword   alphap;


    alphap = lpc_schur_asm(r, k_tmp, p);

    if (a != NULL)
    {
	lpc_refl2pred(k_tmp,a,p);
    }

    return(alphap);   /* alphap in Q15 */
} /* LPC_SCHR */




/* minimum LSP separation-- a global variable */
Shortword lsp_delta = 0;

/* LPC_PRED2LSP
      get LSP coeffs from the predictor coeffs
      Input:
	 a- the predictor coefficients
	 p- the predictor order
      Output:
	 w- the lsp coefficients

      This function uses a DFT to evaluate the P and Q polynomials,
      and is hard-coded to work only for 10th order LPC.

   Q values:
   a - Q12
   w - Q15
*/


#define ONE_Q26 ((Longword)1<<26)
#define X05_Q14 (0.5*(1<<14))

Shortword lpc_pred2lsp(Shortword *a,Shortword *w,Shortword p)
{
    Shortword i,j;
    Shortword p_cof[LPC_ORD+1], q_cof[LPC_ORD+1],
	      p_freq[LPC_ORD+1], q_freq[LPC_ORD+1];
    Longword L_p_cof[LPC_ORD+1], L_q_cof[LPC_ORD+1];
    Longword L_ai,L_api,L_temp;

    /* Generate P' and Q' polynomials   */
    L_p_cof[0] = ONE_Q26;
    L_q_cof[0] = ONE_Q26;

    for ( i = 1; i <= p; i++ ) {
	/* temp = sub(a[i],a[p+1-i]); */     /* temp in Q12 */
        L_ai = L_shr(L_deposit_h(a[i]),2);     
        L_api = L_shr(L_deposit_h(a[p+1-i]),2);     
	L_temp = L_sub(L_ai,L_api);       /* L_temp in Q26 */
	L_p_cof[i] = L_add(L_temp,L_p_cof[i-1]);      /* L_p_cof in Q26 */
	/* temp = add(a[i],a[p+1-i]); */
	L_temp = L_add(L_ai,L_api);
	L_q_cof[i] = L_sub(L_temp,L_q_cof[i-1]);      /* L_q_cof in Q26 */
    }

    /* Convert p_cof and q_cof to short */
    for ( i = 0; i <= p; i++ ) {
      p_cof[i] = round(L_p_cof[i]);     /* p_cof in Q10 */
      q_cof[i] = round(L_q_cof[i]);     /* q_cof in Q10 */
    }

    /* Find root frequencies of LSP polynomials */
    lsp_to_freq(p_cof,p_freq,p);
    lsp_to_freq(q_cof,q_freq,p);

    /* Combine frequencies into single array    */
    w[0] = 0;
    j = 1;
    for (i = 0; i < (shr(p,1)); i++) {
        w[j++] = q_freq[i];      
        w[j++] = p_freq[i];      
    }

    return((Shortword)0);

} /* LPC_PRED2LSP */


/*      Subroutine LSP_init:  Initializes the cos table for lsp_to_freq  */
void    LSP_init()
{
    Shortword i;
    Shortword p2, count, prev_less, j;
    Shortword *p_cos, *p_endlsp, *p_lsp, *p_begcos, *p_endcos;

    for (i = 0; i < DFTLENGTH_D4; i++) {
      lsp_cos[i] = cos_tab[i];
      lsp_cos[i+DFTLENGTH_D4] = -cos_tab[DFTLENGTH_D4-i];
      lsp_cos[i+2*DFTLENGTH_D4] = -cos_tab[i];
      lsp_cos[i+3*DFTLENGTH_D4] = cos_tab[DFTLENGTH_D4-i];
    }

    v_equ( &lsp_cos[DFTLENGTH], lsp_cos, DFTLENGTH);
    v_equ( &lsp_cos[DFTLENGTH * 2], lsp_cos, DFTLENGTH);


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

} /* LPC_INIT */




#undef  DFTLENGTH
#undef  DFTLENGTH_D2
#undef  DFTLENGTH_D4


/* LPC_PRED2REFL
      get refl coeffs from the predictor coeffs
      Input:
	 a- the predictor coefficients
	 p- the predictor order
      Output:
	 k- the reflection coefficients
      Returns:
	 energy - energy of residual signal
   Reference:  Markel and Gray, Linear Prediction of Speech

   Q values:
   *a - Q12
   *k - Q15
   energy - Q15
*/

Shortword lpc_pred2refl(Shortword *a,Shortword *k,Shortword p)
{
    Shortword e;
    Shortword energy = ONE_Q15;
    Shortword i,j;
    Shortword shift;
    Shortword temp, n_temp, n_e;

    Shortword   b[15], b1[15];

    /* equate temporary variables (b = a) */
    for(i=1; i <= p; i++) {
      b[i] = a[i];      
    }

    /* compute reflection coefficients */
    for(i=p; i > 0; i--)
    {
	k[i] = shl(b[i],3);     /* Q12 -> Q15 */
	e = sub(ONE_Q15,mult(k[i],k[i]));
	energy = mult(energy,e);
	for(j=1; j < i; j++) {
            b1[j] = b[j];      
	}

	for(j=1; j < i; j++) {
	    /* b[j] = (b1[j] - k[i]*b1[i-j])/e; */
	    temp = mult(k[i],b1[i-j]);    /* temp in Q12 */
	    temp = sub(b1[j],temp);

	    /* check signs of temp and e before division */
	    n_temp = abs_s(temp);
	    n_e = abs_s(e);
	    shift = norm_s(n_e);
	    n_e = shl(n_e,shift);
	    if (n_temp > n_e) {
	      n_temp = shr(n_temp,1);
	      shift = add(shift,1);
	    }
	    b[j] = shl(divide_s(n_temp,n_e),shift);    /* b[j] in Q12 */

	    if ((temp ^ e) < 0)
	      b[j] = negate(b[j]);
	}
    }


    return(energy);
} /* LPC_PRED2REFL */