语音压缩

开发平台：
Visual C++

vq_lib.c：源码内容
							/*
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 972 480 7442).
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. +33 (0)4 92 94 42 58. 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: +33 493 65 47 16.
*/
#include <stdio.h>
#include "spbstd.h"
#include "mathhalf.h"
#include "wmops.h"
#include "mat.h"
#include "math_lib.h"
#include "lpc.h"
#include "vq.h"
#include "constant.h"
/* VQ_LSPW- compute LSP weighting vector-
    Atal's method:
        From Atal and Paliwal (ICASSP 1991)
        (Note: Paliwal and Atal used w(k)^2(u(k)-u_hat(k))^2,
         and we use w(k)(u(k)-u_hat(k))^2 so our weights are different
         but the method (i.e. the weighted MSE) is the same.
                */
/* Q values:
   w is Q11
   lsp is Q15
   a is Q12 */
Shortword *vq_lspw(Shortword *w, Shortword *lsp, Shortword *a, 
		   Shortword p)
{
    Shortword j;
    Longword L_temp;
    
    for(j = 0; j < p; j++)
      {
	L_temp = lpc_aejw(a,lsp[j],p);    /* L_temp in Q19 */
	w[j] = L_pow_fxp(L_temp,(Shortword)XN03_Q15,19,11);
	if (j == 8)
	  w[j] = mult(w[j], (Shortword)X064_Q15);
	else if (j == 9)
	  w[j] = mult(w[j], (Shortword)X016_Q15);
      }
    return(w);
} /* VQ_LSPW */
/*
    VQ_MS4-
        Tree search multi-stage VQ encoder (optimized for speed
        should be compatible with VQ_MS)
    Synopsis: vq_ms4(cb,u,u_est,levels,ma,stages,p,w,u_hat,a_indices)
        Input:
            cb- one dimensional linear codebook array (codebook is structured 
                as [stages][levels for each stage][p])
            u- dimension p, the parameters to be encoded (u[0..p-1])
            u_est- dimension p, the estimated parameters or mean (if NULL, assume
               estimate is the all zero vector) (u_est[0..p-1])
            levels- the number of levels at each stage (levels[0..stages-1])
            ma- the tree search size (keep ma candidates from one stage to the
               next)
            stages- the number of stages of msvq
            p- the predictor order
            w- the weighting vector (w[0..p-1])
            max_inner- the maximum number of times the swapping procedure
                       in the inner loop can be executed
        Output:
            u_hat-   the reconstruction vector (if non null)
            a_indices- the codebook indices (for each stage) a_indices[0..stages-1]
        Parameters:
*/
#define P_SWAP(x,y,type) do{type u__p;u__p = x;x = y;y = u__p;}while(0)
/* cb is Q17
   u is Q15
   u_est is Q15
   w is Q11
   u_hat is Q15 */
Shortword vq_ms4(Shortword *cb, Shortword *u, Shortword *u_est, 
		 Shortword *levels, Shortword ma, Shortword stages, 
		 Shortword p, Shortword *w, Shortword *u_hat, 
		 Shortword *a_indices, Shortword max_inner)
{
    Shortword tmp,*u_tmp,*uhatw,uhatw_sq;
    Shortword d_cj,d_opt;
    Shortword *d,*p_d,*n_d,*p_distortion,*cb_currentstage,*cbp, *cb_table;
    Shortword *errors,*p_errors,*n_errors,*p_e;
    Shortword i,j,m,s,c,p_max,inner_counter;
    Shortword *indices,*p_indices,*n_indices;
    Shortword *parents,*p_parents,*n_parents;
    Shortword *tmp_p_e;
    Shortword temp;
    Longword L_temp,L_temp1;
    
    /* make sure weights don't get too big */
#define MAXWT  4096           /* w[i] < 2.0 to avoid saturation */
#define MAXWT2 (MAXWT*2)
#define MAXWT4 (MAXWT*4)
    j = 0;    data_move();
    for (i = 0; i < p; i++) {
      if (w[i] > MAXWT4) {
        /* increment complexity for if statement */
        compare_nonzero();
	j = 3;    data_move();
	break;
      }
      else if (w[i] > MAXWT2) {
        /* increment complexity for if statement */
        compare_nonzero();
	j = 2;    data_move();
      }
      else if (w[i] > MAXWT) {
        /* increment complexity for if statement */
        compare_nonzero();
	/* increment complexity for if statement */
        compare_zero();
	if (j == 0) {
	  j = 1;    data_move();
	}
      }
    }
    for (i = 0; i < p; i++) {
      w[i] = shr(w[i],j);    data_move();
    }
    /* allocate memory for the current node and
       parent node (thus, the factors of two everywhere)
       The parents and current nodes are allocated contiguously */
    MEM_ALLOC(MALLOC,indices,2*ma*stages,Shortword);
    MEM_ALLOC(MALLOC,errors,2*ma*p,Shortword);
    MEM_ALLOC(MALLOC,uhatw,p,Shortword);
    MEM_ALLOC(MALLOC,d,2*ma,Shortword);
    MEM_ALLOC(MALLOC,parents,2*ma,Shortword);
    MEM_ALLOC(MALLOC,tmp_p_e,ma*p,Shortword);
    
    /* initialize memory */
    v_zap(indices,(Shortword)(2*stages*ma));
    v_zap(parents,(Shortword)(2*ma));
    
    /* initialize inner loop counter */
    inner_counter = 0;  data_move();
    
    /* set up memory */
    p_indices = &indices[0];  data_move();
    n_indices = &indices[ma*stages];  data_move();
    p_errors = &errors[0];  data_move();
    n_errors = &errors[ma*p];  data_move();
    p_d = &d[0];  data_move();
    n_d = &d[ma];  data_move();
    p_parents = &parents[0];  data_move();
    n_parents = &parents[ma];  data_move();
    
    /* u_tmp is the input vector (i.e. if u_est is non-null, it
       is subtracted off) */
    MEM_ALLOC(MALLOC,u_tmp,p+1,Shortword);
    /* u_tmp is Q15 */
    (void)v_equ(u_tmp,u,p);
    if (u_est)
      (void)v_sub(u_tmp,u_est,p);
    /* change u_tmp from Q15 to Q17 */
    for (j=0; j<p; j++)
      u_tmp[j] = shl(u_tmp[j], 2);  data_move();
    /* L_temp is Q31 */
    L_temp = 0;  data_move();
    for(j = 0 ; j < p; j++) {
      temp = mult(u_tmp[j],w[j]);    /* temp = Q13 */
      L_temp = L_mac(L_temp,temp,u_tmp[j]);
    }
    /* tmp in Q15 */
    tmp = extract_h(L_temp);
    /* set up inital error vectors (i.e. error vectors = u_tmp) */
    for(c=0; c < ma; c++)
    {
	/* n_errors is Q17
	   n_d is Q15 */
	(void)v_equ(&n_errors[c*p],u_tmp,p);
	n_d[c] = tmp;	  data_move();
    }
    /* no longer need memory so free it here */
    MEM_FREE(FREE,u_tmp);
    /* codebook pointer is set to point to first stage */
    cbp = cb;
    /* set m to 1 for the first stage
       and loop over all stages */
    for(m=1,s=0; s < stages; s++)
    {
        /* Save the pointer to the beginning of the
           current stage.  Note: cbp is only incremented in
           one spot, and it is incremented all the way through
           all the stages. */
        cb_currentstage = cbp;
        /* set up pointers to the parent and current nodes */
        P_SWAP(p_indices,n_indices,Shortword*);
        P_SWAP(p_parents,n_parents,Shortword*);
	P_SWAP(p_errors,n_errors,Shortword*);
	P_SWAP(p_d,n_d,Shortword*);
	
        /* p_max is the pointer to the maximum distortion
           node over all candidates.  The so-called worst
           of the best. */
        p_max = 0;  data_move();
	/* store errors in Q15 in tmp_p_e */
	for (i = 0; i < m*p; i++) {
	  tmp_p_e[i] = shr(p_errors[i],2);  data_move();
	}
        /* set the distortions to a large value */
        for(c=0; c < ma; c++) {
	    n_d[c] = SW_MAX;  data_move();
	}
        for(j=0; j < levels[s]; j++)
        {
            /* compute weighted codebook element, increment codebook pointer */
	    /* L_temp is Q31 */
	    L_temp = 0;    data_move();
	    for(i=0 ; i < p; i++,cbp++)
            {
	        /* Q17*Q11<<1 = Q29 */
	        L_temp1 = L_mult(*cbp,w[i]);
		/* uhatw[i] = -2*tmp */
		/* uhatw is Q15 (shift 3 to take care of *2)*/
		uhatw[i] = negate(extract_h(L_shl(L_temp1,3)));
		/* tmp is now Q13 */
		tmp = extract_h(L_temp1);
		L_temp = L_mac(L_temp,*cbp,tmp);
            }
	    /* uhatw_sq is Q15 */
	    uhatw_sq = extract_h(L_temp);
            /* p_e points to the error vectors and p_distortion
               points to the node distortions.  Note: the error
               vectors are contiguous in memory, as are the distortions.
               Thus, the error vector for the 2nd node comes immediately
               after the error for the first node.  (This saves on
               repeated initializations) */
	    /* p_e is Q15 */
	    p_e = tmp_p_e;
            p_distortion = p_d;
	    
            /* iterate over all parent nodes */
            for(c=0; c < m; c++)
            {
		/* L_temp is Q31
		   p_distortion is same Q as n_d 
		   p_e is Q15 */
		L_temp = L_deposit_h(add(*p_distortion++,uhatw_sq));
                for(i=0; i < p; i++) {
		  L_temp = L_mac(L_temp,*p_e++,uhatw[i]);
		}
		/* d_cj is Q15 */
		d_cj = extract_h(L_temp);
                /* determine if d is less than the maximum candidate distortion                   i.e., is the distortion found better than the so-called
                   worst of the best */
		/* increment complexity for if statement */
		compare_nonzero();
                if (d_cj <= n_d[p_max])
                {
                    /* replace the worst with the values just found */
		    /* n_d is now a Q16 */
		    n_d[p_max] = d_cj;    data_move();
                    i = add(shr(extract_l(L_mult(p_max,stages)),1),s);
                    n_indices[i] = j;        data_move();
                    n_parents[p_max] = c;    data_move();
                    /* want to limit the number of times the inner loop
                       is entered (to reduce the *maximum* complexity) */
		    /* increment complexity for if statement */
		    compare_nonzero();
                    if (inner_counter < max_inner)
                    {
                        inner_counter = add(inner_counter,1);
			/* increment complexity for if statement */
			compare_nonzero();
		        if (inner_counter < max_inner)
			{
                            p_max = 0;    data_move();
                            /* find the new maximum */
		            for(i=1; i < ma; i++)
                            {
			        /* increment complexity for if statement */
			        compare_nonzero();
                                if (n_d[i] > n_d[p_max])
                                    p_max = i;    data_move();
		            }
		        }
			else /* inner_counter == max_inner */
                        {
                               /* The inner loop counter now exceeds the
                               maximum, and the inner loop will now not
                               be entered.  Instead of quitting the search
                               or doing something drastic, we simply keep
                               track of the best candidate (rather than the
                               M best) by setting p_max to equal the index
                               of the minimum distortion
                               i.e. only keep one candidate around
                                the MINIMUM distortion */
		            for(i=1; i < ma; i++)
                            {
			        /* increment complexity for if statement */
			        compare_nonzero();
                                if (n_d[i] < n_d[p_max])
                                    p_max = i;    data_move();
		            }
		        }
		    }
                }
            } /* for c */
        } /* for j */
        /* compute the error vectors for each node */
        for(c=0; c < ma; c++)
        {
            /* get the error from the parent node and subtract off
               the codebook value */
	    (void)v_equ(&n_errors[c*p],&p_errors[n_parents[c]*p],p);
	    (void)v_sub(&n_errors[c*p],
			&cb_currentstage[n_indices[c*stages+s]*p],p);
            /* get the indices that were used for the parent node */
            (void)v_equ(&n_indices[c*stages],
			&p_indices[n_parents[c]*stages],s);
        }
        /* increment complexity for if statement */
        compare_nonzero();  data_move();
        m = (m*levels[s] > ma) ? ma : m*levels[s];
    } /* for s */
    /* find the optimum candidate c */
    for(i=1,c=0; i < ma; i++)
    {
        /* increment complexity for if statement */
        compare_nonzero();
        if (n_d[i] < n_d[c])
	    c = i;    data_move();
    }
    d_opt = n_d[c];    data_move();
    
    if (a_indices)
    {
        (void)v_equ(a_indices,&n_indices[c*stages],stages);
    }
    if (u_hat)
    {
        if (u_est)
	    (void)v_equ(u_hat,u_est,p);
        else
	    (void)v_zap(u_hat,p);
        cb_currentstage = cb;
        for(s=0; s < stages; s++)
        {
	  cb_table = &cb_currentstage[n_indices[c*stages+s]*p];	  
	  for (i=0; i<p; i++)
	    u_hat[i] = add(u_hat[i], shr(cb_table[i],2));   data_move();  
	  cb_currentstage += levels[s]*p;
        }
    }
    
    MEM_FREE(FREE,parents);
    MEM_FREE(FREE,d);
    MEM_FREE(FREE,uhatw);
    MEM_FREE(FREE,errors);
    MEM_FREE(FREE,indices);
    MEM_FREE(FREE,tmp_p_e);
    return(d_opt);
} /* vq_ms4 */
/*
    VQ_MSD2-
        Tree search multi-stage VQ decoder
    Synopsis: vq_msd(cb,u,u_est,a,indices,levels,stages,p,conversion)
        Input:
            cb- one dimensional linear codebook array (codebook is structured 
                as [stages][levels for each stage][p])
            indices- the codebook indices (for each stage) indices[0..stages-1]
            levels- the number of levels (for each stage) levels[0..stages-1]
            u_est- dimension p, the estimated parameters or mean (if NULL, assume
               estimate is the all zero vector) (u_est[0..p-1])
            stages- the number of stages of msvq
            p- the predictor order
            conversion- the conversion constant (see lpc.h, lpc_conv.c)
	    diff_Q- the difference between Q value of codebook and u_est
        Output:
            u- dimension p, the decoded parameters (if NULL, use alternate
               temporary storage) (u[0..p-1])
            a- predictor parameters (a[0..p]), if NULL, don't compute
        Returns:
            pointer to reconstruction vector (u)
        Parameters:
*/
Shortword *vq_msd2(Shortword *cb, Shortword *u, Shortword *u_est, 
		   Shortword *a, Shortword *indices, Shortword *levels, 
		   Shortword stages, Shortword p, Shortword conversion,
		   Shortword diff_Q)
{
    Shortword *u_hat,*cb_currentstage,*cb_table;
    Shortword i,j;
    Longword *L_u_hat,L_temp;
    /* allocate memory (if required) */
    MEM_ALLOC(MALLOC,L_u_hat,p,Longword);
    if (u==(Shortword*)NULL)
    {
        MEM_ALLOC(MALLOC,u_hat,p,Shortword);
    }
    else
        u_hat = u;
    /* add estimate on (if non-null), or clear vector */
    if (u_est)
    {
        (void)v_equ(u_hat,u_est,p);
    }
    else
    {
        (void)v_zap(u_hat,p);
    }
    /* put u_hat to a long buffer */
    for (i = 0; i < p; i++)
      L_u_hat[i] = L_shl(L_deposit_l(u_hat[i]),diff_Q);
    /* add the contribution of each stage */
    cb_currentstage = cb;
    for(i=0; i < stages; i++)
    {
      /* (void)v_add(u_hat,&cb_currentstage[indices[i]*p],p); */
      cb_table = &cb_currentstage[indices[i]*p];	  
      for (j = 0; j < p; j++) {
	L_temp = L_deposit_l(cb_table[j]);
        L_u_hat[j] = L_add(L_u_hat[j],L_temp);
      }
      cb_currentstage += levels[i]*p;
    }
    /* convert long buffer back to u_hat */
    for (i = 0; i < p; i++)
      u_hat[i] = extract_l(L_shr(L_u_hat[i],diff_Q));
    MEM_FREE(FREE,L_u_hat);
    return(u);
}
/* VQ_ENC -
   encode vector with full VQ using unweighted Euclidean distance
    Synopsis: vq_enc(cb, u, levels, p, u_hat, indices)
        Input:
            cb- one dimensional linear codebook array
            u- dimension p, the parameters to be encoded (u[0..p-1])
            levels- the number of levels
            p- the predictor order
        Output:
            u_hat-   the reconstruction vector (if non null)
            a_indices- the codebook indices (for each stage) a_indices[0..stages-1]
        Parameters:
*/
/* Q values:
   cb - Q13
   u - Q13
   u_hat - Q13 */
Longword vq_enc(Shortword *cb,Shortword *u,Shortword levels,Shortword p, 
		Shortword *u_hat,Shortword *indices)
{
    Shortword i,j,index;
    Longword d,dmin;
    Shortword *p_cb;
    Shortword temp;
    
    /* Search codebook for minimum distance */
    index = 0;    data_move();
    dmin = LW_MAX;    L_data_move();
    p_cb = cb;
    for (i = 0; i < levels; i++) {
        d = 0;    data_move();
        for (j = 0; j < p; j++) {
	  temp = sub(u[j],*p_cb);
          d = L_mac(d,temp,temp);
          p_cb++;  
        }
	compare_nonzero();
        if (d < dmin) {
            dmin = d;    data_move();
            index = i;    data_move();
        }
    }
    /* Update index and quantized value, and return minimum distance */
    *indices = index;
    v_equ(u_hat,&cb[p*index],p);
    return(dmin);
} /* vq_enc */
/* VQ_FSW - 
   compute the weights for Euclidean distance of Fourier harmonics  */
/* Q values:
   w_fs - Q14
   pitch - Q9 */
void vq_fsw(Shortword *w_fs, Shortword num_harm, Shortword pitch)
{
    Shortword j;
    Shortword tempw0;
    Shortword temp,denom;
    Longword L_temp;
    /* Calculate fundamental frequency */
    /* w0 = TWOPI/pitch */
    /* tempw0 = w0/(0.25*PI) = 8/pitch */
    tempw0 = divide_s((Shortword)EIGHT_Q11,pitch);   /* tempw0 in Q17 */
    for(j=0; j < num_harm; j++) {
      /* Bark-scale weighting */
      /* w_fs[j] = 117.0 / (25.0 + 75.0*
	           pow(1.0 + 1.4*SQR(w0*(j+1)/(0.25*PI)),0.69)) */
      temp = shl(add(j,1),11);    /* Q11 */
      temp = extract_h(L_shl(L_mult(tempw0,temp),1));    /* Q14 */
      temp = mult(temp,temp);    /* Q13 */
      /* L_temp in Q28 */
      L_temp = L_add((Longword)ONE_Q28,L_mult((Shortword)X14_Q14,temp));
      temp = L_pow_fxp(L_temp,(Shortword)X069_Q15,28,13);     /* Q13 */
      denom = add((Shortword)X25_Q6,mult((Shortword)X75_Q8,temp));    /* Q6 */
      w_fs[j] = divide_s((Shortword)X117_Q5,denom);   /* Q14 */   data_move();
    }
}

						