lpc_lib.c

2400bps MELP语音编解码源程序。

/*

2.4 kbps MELP Proposed Federal Standard speech coder

Fixed-point C 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 fixed-point version of the voice codec Mixed Excitation Linear
Prediction (MELP) is based on specifications on the C-language software
simulation contained in GSM 06.06 which is protected by copyright and
is the property of the European Telecommunications Standards Institute
(ETSI). This standard is available from the ETSI publication office
tel. [phone]. ETSI has granted a license to United States
Department of Defense to use the C-language software simulation contained
in GSM 06.06 for the purposes of the development of a fixed-point
version of the voice codec Mixed Excitation Linear Prediction (MELP).
Requests for authorization to make other use of the GSM 06.06 or
otherwise distribute or modify them need to be addressed to the ETSI
Secretariat fax: [phone].

*/
/*  

  lpc_lib.c: LPC subroutines.

*/

/*  compiler include files  */
#include        <stdio.h>
#include        <stdlib.h>
#include        <math.h>
#include "spbstd.h"
#include "mathhalf.h"
#include "mathdp31.h"
#include "math_lib.h"
#include "wmops.h"
#include "mat.h"
#include "lpc.h"
#include "constant.h"

extern int complexity;
extern Shortword frames;
extern FILE *fp_outdat;

#define PRINT 1

#define ALMOST_ONE_Q14   16382    /* ((1<<14)-2) */

/* Lag window coefficients */
Shortword lagw_cof[10] = {
32756,32721,32663,32582,32478,32351,32201,32030,31837,31622};

/*
    LPC_ACOR
        Compute autocorrelations based on windowed speech frame

    Synopsis: lpc_acor(s,window,r,hf_correction,lag_win,p,nx)
        Input:
            input- input vector (npts samples, s[0..npts-1])
            win_cof- window vector (npts samples, s[0..npts-1])
            hf_correction- high frequency correction value
            lag_win- lag window value (for bandwidth expansion
                     factor is BeqT, where T is the sampling rate)
		     NOT used in this function
            order- order of lpc filter
            npts- number of elements in window
        Output:
            r- output autocorrelation vector (order+1 samples, r[0..order])

    Q values:
    input - Q0
    win_cof - Q15
    hf_correction - Q15
*/

void lpc_acor(Shortword input[],Shortword win_cof[],
	      Shortword r[],Shortword hf_correction,
	      Shortword lag_win,Shortword order,Shortword npts)
{
    Shortword i,j;
    Longword L_temp;
    Shortword *inputw;
    Shortword norm_var, scale_fact, temp;

    /* window optimized for speed and readability
       does windowing and autocorrelation sequentially
       and in the usual manner
    */ 

    MEM_ALLOC(MALLOC,inputw,npts,Shortword);

    if (win_cof) {
      for(i=0; i < npts; i++) {
	inputw[i] = mult(win_cof[i],shr(input[i],4));   data_move();
      }
    }
    else {
      v_equ_shr(inputw,input,4,npts);
    }
    
    /* Find scaling factor */
    L_temp = L_v_magsq(inputw,npts,0,1);
    compare_zero();
    if (L_temp) {
      norm_var = norm_l(L_temp);
      norm_var = sub(4,shr(norm_var,1));
      if (norm_var < 0)
	norm_var = 0;
    }
    else
      norm_var = 0;

    /* increment complexity for if statement */
    compare_zero();
    if (win_cof) {
      for(i=0; i < npts; i++) {
	inputw[i] = shr(mult(win_cof[i],input[i]),norm_var);    data_move();
      }
    }
    else
      v_equ_shr(inputw,input,norm_var,npts);

    /* Compute r[0] */
    L_temp = L_v_magsq(inputw,npts,0,1);
    compare_zero();
    if (L_temp > 0) {
      /* normalize with 1 bit of headroom */
      norm_var = norm_l(L_temp);
      norm_var = sub(norm_var,1);
      L_temp = L_shl(L_temp,norm_var);

      /* High frequency correction */
      L_temp = L_add(L_temp,L_mpy_ls(L_temp,hf_correction));

      /* normalize result */
      temp = norm_s(extract_h(L_temp));
      L_temp = L_shl(L_temp,temp);
      norm_var = add(norm_var,temp);
      r[0] = round(L_temp);

      /* Multiply by 1/r[0] for full normalization */
      scale_fact = divide_s(ALMOST_ONE_Q14,r[0]);
      L_temp = L_shl(L_mpy_ls(L_temp,scale_fact),1);
      r[0] = round(L_temp);

    }
    else {
      norm_var = 0;    data_move();
      r[0] = ONE_Q15;   /* 1 in Q15 */    data_move();
      scale_fact = 0;      data_move();
    }	 

    /* Compute remaining autocorrelation terms */
    for(j=1; j <= order; j++)
    {
	L_temp = 0;    data_move();
        for(i=j; i < npts; i++) {
	    L_temp = L_mac(L_temp,inputw[i],inputw[i-j]);
	}
	L_temp = L_shl(L_temp,norm_var);

	/* Scaling */
	L_temp = L_shl(L_mpy_ls(L_temp,scale_fact),1);

	/* Lag windowing */
	L_temp = L_mpy_ls(L_temp,lagw_cof[j-1]);

	r[j] = round(L_temp);
    }
    
    MEM_FREE(FREE,inputw);

} /* lpc_acor */

/* 
    Name: lpc_aejw- Compute square of A(z) evaluated at exp(jw)
    Description:
        Compute the magnitude squared of the z-transform of 

        A(z) = 1+a(1)z^-1 + ... +a(p)z^-p

        evaluated at z=exp(jw)
     Inputs:
        a- LPC filter (a[0] is undefined, a[1..p])
        w- radian frequency
        p- predictor order
     Returns:
        |A(exp(jw))|^2
     See_Also: cos(3), sin(3)
     Includes:
        spbstd.h
        lpc.h
     Systems and Info. Science Lab
     Copyright (c) 1995 by Texas Instruments, Inc.  All rights reserved.

*/

/* Q values
   --------
   a - Q12
   w - Q15
   return - Q19
*/

/* lower limit of return value to prevent overflow of weighting function */
#define LOW_LIMIT 54

Longword lpc_aejw(Shortword *a,Shortword w, Shortword p)
{
    Shortword i;
    Shortword c_re,c_im;
    Shortword cs,sn,temp;
    Shortword temp1,temp2;
    Longword L_temp;
    

    if (p==0)
        return((Longword)ONE_Q19);

    /* use horners method
    A(exp(jw)) = 1+ e(-jw)[a(1)+e(-jw)[a(2)+e(-jw)[a(3)+..
            ...[a(p-1)+e(-jw)a(p)]]]]
    */

    cs = cos_fxp(w);     data_move();    /* temp in Q15 */
    sn = negate(sin_fxp(w));     data_move();    /* temp in Q15 */

    c_re = shr(mult(cs,a[p]),3);     data_move();  /* -> Q9 */
    c_im = shr(mult(sn,a[p]),3);     data_move();  /* -> Q9 */

    for(i=p-1; i > 0; i--)
    {
        /* add a[i] */
        temp = shr(a[i],3);
        c_re = add(c_re,temp);     data_move();

        /* multiply by exp(-jw) */
	temp=c_im;
	temp1 = mult(cs,temp);     data_move();   /* temp1 in Q9 */
	temp2 = mult(sn,c_re);     data_move();   /* temp2 in Q9 */
        c_im = add(temp1, temp2);     data_move();
	temp1 = mult(cs,c_re);     data_move();   /* temp1 in Q9 */
	temp2 = mult(sn,temp);     data_move();   /* temp2 in Q9 */
        c_re = sub(temp1,temp2);     data_move();
    }

    /* add one */
    c_re = add(c_re, ONE_Q9);     data_move();

    /* L_temp in Q19 */
    L_temp = L_add(L_mult(c_re,c_re),L_mult(c_im,c_im));
    if (L_temp<LOW_LIMIT)
      L_temp=(Longword)LOW_LIMIT;
    
    return(L_temp);
} /* LPC_AEJW */

#undef LOW_LIMIT


/*
    Name: lpc_bwex- Move the zeros of A(z) toward the origin.
    Aliases: lpc_bw_expand
    Description:
        Expand the zeros of the LPC filter by gamma, which
        moves each zero radially into the origin.
<nf>
        for j = 1 to p
            aw[j] = a[j]*gamma^j
</nf>
        (Can also be used to perform an exponential windowing procedure).
    Inputs:
        a- lpc vector (order p, a[1..p])
        gamma- the bandwidth expansion factor
        p- order of lpc filter
    Outputs:
        aw- the bandwidth expanded LPC filter
    Returns: NULL
    See_Also: lpc_lagw(3l)
    Includes:
        spbstd.h
        lpc.h

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

    Q values:
    *a,*aw - Q12
    gamma - Q15
    gk - Q15
*/

Shortword lpc_bwex(Shortword *a, Shortword *aw, Shortword gamma, Shortword p)
{
    Shortword i;
    Shortword gk;

    gk = gamma;      data_move();
    for(i=1; i <= p; i++) {
        aw[i] = mult(a[i],gk);
	gk = mult(gk,gamma);
    }
    return((Shortword)0);
}


/*
    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 */
	        compare_nonzero();
		d = sub(w[i+1],w[i]);     data_move();
                if (d < delta)
                {
		    data_move();    data_move();
                    step1 = step2 = shr(sub(delta,d),1);

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

		    /* increment complexity for if statement */
		    compare_nonzero();    compare_nonzero();
                    if (i==(p-1) && (w[i+1] > sub(ONE_Q15,delta))) {
                        step2 = shr(sub(ONE_Q15,w[i+1]),1);    data_move();
		    }
                    else {
		      /* increment complexity for if statement */
		      compare_nonzero();
		      if (i < (p-1))
		      {
		        /* increment complexity for if statement */
		        compare_nonzero();
			temp = sub(w[i+2],w[i+1]);
                        if (temp < delta) {
			  step2 = 0;    data_move();
			}
                        else {
			  /* increment complexity for if statement */
			  compare_nonzero();
			  if (temp < shl(delta,1)) {
                            step2 = shr(sub(temp,delta),1);
			  }
			}
		      }
                    }
                    w[i] = sub(w[i],step1);    data_move();
                    w[i+1] = add(w[i+1],step2);    data_move();
                    
                }
            }
        }
    }  
    
    /* Debug: check if the minimum separation rule was met */
    temp = mult(32440,delta);   /* temp = 0.99*delta */
    for(j=1; j < p; j++)
      if ((w[j+1]-w[j]) < temp)
          (void)fprintf(stderr,"%s: LSPs not separated by enough (line %d)\n",
              __FILE__,__LINE__);
    
    if (unsorted)
      (void)fprintf(stderr,"%s: Fxp LSPs still unsorted (line %d)\n",
                    __FILE__,__LINE__);

    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 i,j;
    Shortword shift,alphap,*k;
    Longword L_temp, *y1, *y2;
    Shortword n_temp,n_y2,nr1,nr0;
    
 
    MEM_ALLOC(MALLOC,y1,p+2,Longword);
    MEM_ALLOC(MALLOC,y2,p+2,Longword);

    if (k_tmp == NULL)
 	MEM_ALLOC(MALLOC,k,p+1,Shortword);
    else
	k = k_tmp;

    nr1 = abs_s(r[1]);    data_move();
    nr0 = abs_s(r[0]);    data_move();

    k[1] = divide_s(nr1,nr0);    data_move();

    /* increment complexity for if statement */
    logic(); compare_zero();
    /* if (((r[1] < 0) && (r[0] < 0)) || ((r[1] > 0) && (r[0] > 0))) */
    if ((r[1] ^ r[0]) >= 0)
    {
      k[1] = negate(k[1]);    data_move();
    }
    alphap = mult(r[0],sub(ONE_Q15,mult(k[1],k[1])));    data_move();

    y2[1] = L_deposit_h(r[1]);
    y2[2] = L_add(L_deposit_h(r[0]),L_mult(k[1],r[1]));