开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 通讯/手机编程 > 语音压缩 > 查看源码
lpc_lib.c
资源名称：melpfix.rar　[点击查看]
上传用户：cxx_68
上传日期：2021-02-21
资源大小：161k
文件大小：29k
源码类别：

语音压缩

开发平台：

Visual C++

lpc_lib.c：源码内容
  1. /*
  3. 2.4 kbps MELP Proposed Federal Standard speech coder
  5. Fixed-point C code, version 1.0
  7. Copyright (c) 1998, Texas Instruments, Inc.  
  9. Texas Instruments has intellectual property rights on the MELP
  10. algorithm.  The Texas Instruments contact for licensing issues for
  11. commercial and non-government use is William Gordon, Director,
  12. Government Contracts, Texas Instruments Incorporated, Semiconductor
  13. Group (phone 972 480 7442).
  15. The fixed-point version of the voice codec Mixed Excitation Linear
  16. Prediction (MELP) is based on specifications on the C-language software
  17. simulation contained in GSM 06.06 which is protected by copyright and
  18. is the property of the European Telecommunications Standards Institute
  19. (ETSI). This standard is available from the ETSI publication office
  20. tel. +33 (0)4 92 94 42 58. ETSI has granted a license to United States
  21. Department of Defense to use the C-language software simulation contained
  22. in GSM 06.06 for the purposes of the development of a fixed-point
  23. version of the voice codec Mixed Excitation Linear Prediction (MELP).
  24. Requests for authorization to make other use of the GSM 06.06 or
  25. otherwise distribute or modify them need to be addressed to the ETSI
  26. Secretariat fax: +33 493 65 47 16.
  28. */
  29. /*  
  31.   lpc_lib.c: LPC subroutines.
  33. */
  35. /*  compiler include files  */
  36. #include        <stdio.h>
  37. #include        <stdlib.h>
  38. #include        <math.h>
  39. #include "spbstd.h"
  40. #include "mathhalf.h"
  41. #include "mathdp31.h"
  42. #include "math_lib.h"
  43. #include "wmops.h"
  44. #include "mat.h"
  45. #include "lpc.h"
  46. #include "constant.h"
  48. extern int complexity;
  49. extern Shortword frames;
  50. extern FILE *fp_outdat;
  52. #define PRINT 1
  54. #define ALMOST_ONE_Q14   16382    /* ((1<<14)-2) */
  56. /* Lag window coefficients */
  57. Shortword lagw_cof[10] = {
  58. 32756,32721,32663,32582,32478,32351,32201,32030,31837,31622};
  60. /*
  61.     LPC_ACOR
  62.         Compute autocorrelations based on windowed speech frame
  64.     Synopsis: lpc_acor(s,window,r,hf_correction,lag_win,p,nx)
  65.         Input:
  66.             input- input vector (npts samples, s[0..npts-1])
  67.             win_cof- window vector (npts samples, s[0..npts-1])
  68.             hf_correction- high frequency correction value
  69.             lag_win- lag window value (for bandwidth expansion
  70.                      factor is BeqT, where T is the sampling rate)
  71.      NOT used in this function
  72.             order- order of lpc filter
  73.             npts- number of elements in window
  74.         Output:
  75.             r- output autocorrelation vector (order+1 samples, r[0..order])
  77.     Q values:
  78.     input - Q0
  79.     win_cof - Q15
  80.     hf_correction - Q15
  81. */
  83. void lpc_acor(Shortword input[],Shortword win_cof[],
  84.       Shortword r[],Shortword hf_correction,
  85.       Shortword lag_win,Shortword order,Shortword npts)
  86. {
  87.     Shortword i,j;
  88.     Longword L_temp;
  89.     Shortword *inputw;
  90.     Shortword norm_var, scale_fact, temp;
  92.     /* window optimized for speed and readability
  93.        does windowing and autocorrelation sequentially
  94.        and in the usual manner
  95.     */ 
  97.     MEM_ALLOC(MALLOC,inputw,npts,Shortword);
  99.     if (win_cof) {
  100.       for(i=0; i < npts; i++) {
  101. inputw[i] = mult(win_cof[i],shr(input[i],4));   data_move();
  102.       }
  103.     }
  104.     else {
  105.       v_equ_shr(inputw,input,4,npts);
  106.     }
  107.     
  108.     /* Find scaling factor */
  109.     L_temp = L_v_magsq(inputw,npts,0,1);
  110.     compare_zero();
  111.     if (L_temp) {
  112.       norm_var = norm_l(L_temp);
  113.       norm_var = sub(4,shr(norm_var,1));
  114.       if (norm_var < 0)
  115. norm_var = 0;
  116.     }
  117.     else
  118.       norm_var = 0;
  120.     /* increment complexity for if statement */
  121.     compare_zero();
  122.     if (win_cof) {
  123.       for(i=0; i < npts; i++) {
  124. inputw[i] = shr(mult(win_cof[i],input[i]),norm_var);    data_move();
  125.       }
  126.     }
  127.     else
  128.       v_equ_shr(inputw,input,norm_var,npts);
  130.     /* Compute r[0] */
  131.     L_temp = L_v_magsq(inputw,npts,0,1);
  132.     compare_zero();
  133.     if (L_temp > 0) {
  134.       /* normalize with 1 bit of headroom */
  135.       norm_var = norm_l(L_temp);
  136.       norm_var = sub(norm_var,1);
  137.       L_temp = L_shl(L_temp,norm_var);
  139.       /* High frequency correction */
  140.       L_temp = L_add(L_temp,L_mpy_ls(L_temp,hf_correction));
  142.       /* normalize result */
  143.       temp = norm_s(extract_h(L_temp));
  144.       L_temp = L_shl(L_temp,temp);
  145.       norm_var = add(norm_var,temp);
  146.       r[0] = round(L_temp);
  148.       /* Multiply by 1/r[0] for full normalization */
  149.       scale_fact = divide_s(ALMOST_ONE_Q14,r[0]);
  150.       L_temp = L_shl(L_mpy_ls(L_temp,scale_fact),1);
  151.       r[0] = round(L_temp);
  153.     }
  154.     else {
  155.       norm_var = 0;    data_move();
  156.       r[0] = ONE_Q15;   /* 1 in Q15 */    data_move();
  157.       scale_fact = 0;      data_move();
  158.     }  
  160.     /* Compute remaining autocorrelation terms */
  161.     for(j=1; j <= order; j++)
  162.     {
  163. L_temp = 0;    data_move();
  164.         for(i=j; i < npts; i++) {
  165.     L_temp = L_mac(L_temp,inputw[i],inputw[i-j]);
  166. }
  167. L_temp = L_shl(L_temp,norm_var);
  169. /* Scaling */
  170. L_temp = L_shl(L_mpy_ls(L_temp,scale_fact),1);
  172. /* Lag windowing */
  173. L_temp = L_mpy_ls(L_temp,lagw_cof[j-1]);
  175. r[j] = round(L_temp);
  176.     }
  177.     
  178.     MEM_FREE(FREE,inputw);
  180. } /* lpc_acor */
  182. /* 
  183.     Name: lpc_aejw- Compute square of A(z) evaluated at exp(jw)
  184.     Description:
  185.         Compute the magnitude squared of the z-transform of 
  187.         A(z) = 1+a(1)z^-1 + ... +a(p)z^-p
  189.         evaluated at z=exp(jw)
  190.      Inputs:
  191.         a- LPC filter (a[0] is undefined, a[1..p])
  192.         w- radian frequency
  193.         p- predictor order
  194.      Returns:
  195.         |A(exp(jw))|^2
  196.      See_Also: cos(3), sin(3)
  197.      Includes:
  198.         spbstd.h
  199.         lpc.h
  200.      Systems and Info. Science Lab
  201.      Copyright (c) 1995 by Texas Instruments, Inc.  All rights reserved.
  203. */
  205. /* Q values
  206.    --------
  207.    a - Q12
  208.    w - Q15
  209.    return - Q19
  210. */
  212. /* lower limit of return value to prevent overflow of weighting function */
  213. #define LOW_LIMIT 54
  215. Longword lpc_aejw(Shortword *a,Shortword w, Shortword p)
  216. {
  217.     Shortword i;
  218.     Shortword c_re,c_im;
  219.     Shortword cs,sn,temp;
  220.     Shortword temp1,temp2;
  221.     Longword L_temp;
  222.     
  224.     if (p==0)
  225.         return((Longword)ONE_Q19);
  227.     /* use horners method
  228.     A(exp(jw)) = 1+ e(-jw)[a(1)+e(-jw)[a(2)+e(-jw)[a(3)+..
  229.             ...[a(p-1)+e(-jw)a(p)]]]]
  230.     */
  232.     cs = cos_fxp(w);     data_move();    /* temp in Q15 */
  233.     sn = negate(sin_fxp(w));     data_move();    /* temp in Q15 */
  235.     c_re = shr(mult(cs,a[p]),3);     data_move();  /* -> Q9 */
  236.     c_im = shr(mult(sn,a[p]),3);     data_move();  /* -> Q9 */
  238.     for(i=p-1; i > 0; i--)
  239.     {
  240.         /* add a[i] */
  241.         temp = shr(a[i],3);
  242.         c_re = add(c_re,temp);     data_move();
  244.         /* multiply by exp(-jw) */
  245. temp=c_im;
  246. temp1 = mult(cs,temp);     data_move();   /* temp1 in Q9 */
  247. temp2 = mult(sn,c_re);     data_move();   /* temp2 in Q9 */
  248.         c_im = add(temp1, temp2);     data_move();
  249. temp1 = mult(cs,c_re);     data_move();   /* temp1 in Q9 */
  250. temp2 = mult(sn,temp);     data_move();   /* temp2 in Q9 */
  251.         c_re = sub(temp1,temp2);     data_move();
  252.     }
  254.     /* add one */
  255.     c_re = add(c_re, ONE_Q9);     data_move();
  257.     /* L_temp in Q19 */
  258.     L_temp = L_add(L_mult(c_re,c_re),L_mult(c_im,c_im));
  259.     if (L_temp<LOW_LIMIT)
  260.       L_temp=(Longword)LOW_LIMIT;
  261.     
  262.     return(L_temp);
  263. } /* LPC_AEJW */
  265. #undef LOW_LIMIT
  268. /*
  269.     Name: lpc_bwex- Move the zeros of A(z) toward the origin.
  270.     Aliases: lpc_bw_expand
  271.     Description:
  272.         Expand the zeros of the LPC filter by gamma, which
  273.         moves each zero radially into the origin.
  274. <nf>
  275.         for j = 1 to p
  276.             aw[j] = a[j]*gamma^j
  277. </nf>
  278.         (Can also be used to perform an exponential windowing procedure).
  279.     Inputs:
  280.         a- lpc vector (order p, a[1..p])
  281.         gamma- the bandwidth expansion factor
  282.         p- order of lpc filter
  283.     Outputs:
  284.         aw- the bandwidth expanded LPC filter
  285.     Returns: NULL
  286.     See_Also: lpc_lagw(3l)
  287.     Includes:
  288.         spbstd.h
  289.         lpc.h
  291.     Systems and Info. Science Lab
  292.     Copyright (c) 1995 by Texas Instruments, Inc.  All rights reserved.
  294.     Q values:
  295.     *a,*aw - Q12
  296.     gamma - Q15
  297.     gk - Q15
  298. */
  300. Shortword lpc_bwex(Shortword *a, Shortword *aw, Shortword gamma, Shortword p)
  301. {
  302.     Shortword i;
  303.     Shortword gk;
  305.     gk = gamma;      data_move();
  306.     for(i=1; i <= p; i++) {
  307.         aw[i] = mult(a[i],gk);
  308. gk = mult(gk,gamma);
  309.     }
  310.     return((Shortword)0);
  311. }
  314. /*
  315.     Name: lpc_clmp- Sort and ensure minimum separation in LSPs.
  316.     Aliases: lpc_clamp
  317.     Description:
  318.         Ensure that all LSPs are ordered and separated
  319.         by at least delta.  The algorithm isn't guarenteed
  320.         to work, so it prints an error message when it fails
  321.         to sort the LSPs properly.
  322.     Inputs:
  323.         w- lsp vector (order p, w[1..p])
  324.         delta- the clamping factor
  325.         p- order of lpc filter
  326.     Outputs:
  327.         w- the sorted and clamped lsps
  328.     Returns: NULL
  329.     See_Also:
  330.     Includes:
  331.         spbstd.h
  332.         lpc.h
  333.     Bugs: 
  334.         Currently only supports 10 loops, which is too
  335.         complex and perhaps unneccesary.
  337.     Systems and Info. Science Lab
  338.     Copyright (c) 1995 by Texas Instruments, Inc.  All rights reserved.
  340. */
  341. /*  Q values:
  342.     w - Q15
  343.     delta - Q15 */
  345. #define MAX_LOOPS 10
  347. Shortword lpc_clmp(Shortword *w, Shortword delta, Shortword p)
  348. {
  349.     Shortword i,j,unsorted;
  350.     Shortword temp,d,step1,step2;
  352.     /* sort the LSPs for 10 loops */  
  353.     for (j=0,unsorted=TRUE; unsorted && (j < MAX_LOOPS); j++)
  354.     {
  355.         for(i=1,unsorted=FALSE; i < p; i++)
  356.             if (w[i] > w[i+1])
  357.             {
  358.                 temp = w[i+1];
  359.                 w[i+1] = w[i];
  360.                 w[i] = temp;
  361.                 unsorted = TRUE;
  362.             }
  363.     }
  365.     /* ensure minimum separation */
  366.     if (!unsorted) 
  367.     {
  368.         for(j=0; j < MAX_LOOPS; j++)
  369.         {
  370.             for(i=1; i < p; i++)
  371.             {
  372.         /* increment complexity for if statement */
  373.         compare_nonzero();
  374. d = sub(w[i+1],w[i]);     data_move();
  375.                 if (d < delta)
  376.                 {
  377.     data_move();    data_move();
  378.                     step1 = step2 = shr(sub(delta,d),1);
  380.     /* increment complexity for if statement */
  381.     compare_nonzero();    compare_nonzero();
  382.                     if (i==1 && (w[i] < delta)) {
  383.       step1 = shr(w[i],1);    data_move();
  384.     }
  385.                     else {
  386.       /* increment complexity for if statement */
  387.       compare_nonzero();
  388.       if (i > 1)
  389.       {
  390.         /* increment complexity for if statement */
  391.         compare_nonzero();
  392. temp = sub(w[i],w[i-1]);
  393.                         if (temp < delta) {
  394.   step1 = 0;    data_move();
  395. }
  396.                         else {
  397.   /* increment complexity for if statement */
  398.   compare_nonzero();    compare_nonzero();
  399.   if (temp < shl(delta,1)) {
  400.     step1 = shr(sub(temp,delta),1);
  401.   }
  402. }
  403.       }
  404.                     }
  406.     /* increment complexity for if statement */
  407.     compare_nonzero();    compare_nonzero();
  408.                     if (i==(p-1) && (w[i+1] > sub(ONE_Q15,delta))) {
  409.                         step2 = shr(sub(ONE_Q15,w[i+1]),1);    data_move();
  410.     }
  411.                     else {
  412.       /* increment complexity for if statement */
  413.       compare_nonzero();
  414.       if (i < (p-1))
  415.       {
  416.         /* increment complexity for if statement */
  417.         compare_nonzero();
  418. temp = sub(w[i+2],w[i+1]);
  419.                         if (temp < delta) {
  420.   step2 = 0;    data_move();
  421. }
  422.                         else {
  423.   /* increment complexity for if statement */
  424.   compare_nonzero();
  425.   if (temp < shl(delta,1)) {
  426.                             step2 = shr(sub(temp,delta),1);
  427.   }
  428. }
  429.       }
  430.                     }
  431.                     w[i] = sub(w[i],step1);    data_move();
  432.                     w[i+1] = add(w[i+1],step2);    data_move();
  433.                     
  434.                 }
  435.             }
  436.         }
  437.     }  
  438.     
  439.     /* Debug: check if the minimum separation rule was met */
  440.     temp = mult(32440,delta);   /* temp = 0.99*delta */
  441.     for(j=1; j < p; j++)
  442.       if ((w[j+1]-w[j]) < temp)
  443.           (void)fprintf(stderr,"%s: LSPs not separated by enough (line %d)n",
  444.               __FILE__,__LINE__);
  445.     
  446.     if (unsorted)
  447.       (void)fprintf(stderr,"%s: Fxp LSPs still unsorted (line %d)n",
  448.                     __FILE__,__LINE__);
  450.     return((Shortword)0);
  451. } /* LPC_CLMP */
  454. /*  
  455.     Name: lpc_schr- Schur recursion (autocorrelations to refl coef)
  456.     Aliases: lpc_schur
  457.     Description:
  458.         Compute reflection coefficients from autocorrelations
  459.         based on schur recursion.  Will also compute predictor
  460.         parameters by calling lpc_refl2pred(3l) if necessary.
  461.     Inputs:
  462.         r- autocorrelation vector (r[0..p]).
  463.         p- order of lpc filter.
  464.     Outputs:
  465.         a-     predictor parameters    (can be NULL)
  466.         k_tmp- reflection coefficients (can be NULL)
  467.     Returns:
  468.         alphap- the minimum residual energy
  469.     Includes:
  470.         spbstd.h
  471.         lpc.h
  472.     See_Also:
  473.         lpc_refl2pred(3l) in lpc.h or lpc(3l)
  475.     Q values:
  476.     r - Q0
  477.     a - Q12
  478.     k_tmp - Q15
  479. */
  481. Shortword lpc_schr(Shortword *r, Shortword *a, Shortword *k_tmp, Shortword p)
  482. {
  483.     Shortword i,j;
  484.     Shortword shift,alphap,*k;
  485.     Longword L_temp, *y1, *y2;
  486.     Shortword n_temp,n_y2,nr1,nr0;
  487.     
  488.  
  489.     MEM_ALLOC(MALLOC,y1,p+2,Longword);
  490.     MEM_ALLOC(MALLOC,y2,p+2,Longword);
  492.     if (k_tmp == NULL)
  493.   MEM_ALLOC(MALLOC,k,p+1,Shortword);
  494.     else
  495. k = k_tmp;
  497.     nr1 = abs_s(r[1]);    data_move();
  498.     nr0 = abs_s(r[0]);    data_move();
  500.     k[1] = divide_s(nr1,nr0);    data_move();
  502.     /* increment complexity for if statement */
  503.     logic(); compare_zero();
  504.     /* if (((r[1] < 0) && (r[0] < 0)) || ((r[1] > 0) && (r[0] > 0))) */
  505.     if ((r[1] ^ r[0]) >= 0)
  506.     {
  507.       k[1] = negate(k[1]);    data_move();
  508.     }
  509.     alphap = mult(r[0],sub(ONE_Q15,mult(k[1],k[1])));    data_move();
  511.     y2[1] = L_deposit_h(r[1]);
  512.     y2[2] = L_add(L_deposit_h(r[0]),L_mult(k[1],r[1]));
  514.     for(i=2; i <= p; i++)
  515.     {
  516. y1[1] = L_deposit_h(r[i]);
  517. L_temp = L_deposit_h(r[i]);
  519.         for(j=1; j < i; j++)
  520.         {  
  521.     y1[j+1] = L_add(y2[j], L_mpy_ls(L_temp, k[j]));    L_data_move();
  522.     L_temp = L_add (L_temp, L_mpy_ls(y2[j], k[j]));
  523.     y2[j] = y1[j];    L_data_move();
  524. }
  526.         /* k[i] = -temp/y2[i]; */
  527.         if (L_temp > y2[i]) {
  528.   for (j = i; j <= p; j++)
  529.     k[j] = 0;
  530. #if (PRINT)
  531.   printf("nERROR: numerator bigger than denominator in lpc_schrn");
  532.   printf("Reflection coefficients used in frame %d:n",frames);
  533.   for (i = 1; i <= p; i++) {
  534.     printf(" %5.8f ",(float)k[i]/(1<<15));
  535.   }
  536.   printf("n");
  537. #endif
  538.   break;
  539. }
  541. shift = norm_l(y2[i]);
  542. n_temp = abs_s(extract_h(L_shl(L_temp,shift)));    data_move();
  543. n_y2 = abs_s(extract_h(L_shl(y2[i],shift)));    data_move();
  545. k[i] = divide_s(n_temp,n_y2);    data_move();
  546. /* increment complexity for if statement */
  547. L_logic(); compare_zero();
  548. if ((L_temp ^ y2[i]) >= 0)
  549. {
  550.     k[i] = negate(k[i]);    data_move();
  551. }
  553. y2[i+1] = L_add(y2[i], L_mpy_ls(L_temp, k[i]));    L_data_move();
  554. y2[i] = y1[i];    L_data_move();
  556. alphap = mult(alphap,sub(ONE_Q15,mult(k[i],k[i])));    data_move();
  557.     }
  558.     
  559.     if (a != NULL)
  560.     {
  561. lpc_refl2pred(k,a,p);
  562.     }
  564.     if (k_tmp == NULL) 
  565. MEM_FREE(FREE,k);
  567.     MEM_FREE(FREE,y2);
  568.     MEM_FREE(FREE,y1);
  569.     
  570.     return(alphap);   /* alphap in Q15 */
  571. } /* LPC_SCHR */
  574. /* LPC_REFL2PRED
  575.       get predictor coefficients from the reflection coeffs
  576.    Synopsis: lpc_refl2pred(k,a,p)
  578.       Input:
  579.          k- the reflection coeffs
  580.          p- the predictor order
  581.       Output:
  582.          a- the predictor coefficients
  583.    Reference:  Markel and Gray, Linear Prediction of Speech
  585.    Q values:
  586.    k - Q15
  587.    a - Q12
  588. */
  590. Shortword lpc_refl2pred(Shortword *k,Shortword *a,Shortword p)
  591. {
  592.     Shortword i,j;
  593.     Shortword *a1;
  595.     MEM_ALLOC(MALLOC,a1,p+1,Shortword);
  597.     for(i=1; i <= p; i++)
  598.     {
  599.         /* refl to a recursion */
  600.         a[i] = shift_r(k[i],-3);      /* a in Q12 */    data_move();
  601.         for(j=1; j < i; j++) {
  602.   a1[j] = a[j];    data_move();
  603. }
  604.         for(j=1; j < i; j++) {
  605.             a[j] = add(a1[j],mult(k[i],a1[i-j]));    data_move();
  606. }
  607.     }
  609.     MEM_FREE(FREE,a1);
  610.     return((Shortword)0);
  611. } /* LPC_REFL2PRED */
  614. /* minimum LSP separation-- a global variable */
  615. Shortword lsp_delta = 0;
  617. /* LPC_PRED2LSP
  618.       get LSP coeffs from the predictor coeffs
  619.       Input:
  620.          a- the predictor coefficients
  621.          p- the predictor order
  622.       Output:
  623.          w- the lsp coefficients
  625.       This function uses a DFT to evaluate the P and Q polynomials,
  626.       and is hard-coded to work only for 10th order LPC.
  628.    Q values:
  629.    a - Q12
  630.    w - Q15
  631. */
  632. static void lsp_to_freq(Shortword lsp[], Shortword w[], Shortword order);
  634. #define LPC_ORD 10
  636. Shortword lpc_pred2lsp(Shortword *a,Shortword *w,Shortword p)
  637. {
  638.     Shortword i,j;
  639.     Shortword p_cof[LPC_ORD+1], q_cof[LPC_ORD+1],
  640.               p_freq[LPC_ORD+1], q_freq[LPC_ORD+1];
  641.     Longword L_p_cof[LPC_ORD+1], L_q_cof[LPC_ORD+1];
  642.     Longword L_ai,L_api,L_temp;
  644.     /* Generate P' and Q' polynomials   */
  645.     L_p_cof[0] = (Longword)ONE_Q26;
  646.     L_q_cof[0] = (Longword)ONE_Q26;
  647.     for ( i = 1; i <= p; i++ ) {
  648.         /* temp = sub(a[i],a[p+1-i]); */     /* temp in Q12 */
  649.         L_ai = L_shr(L_deposit_h(a[i]),2);     L_data_move();
  650. L_api = L_shr(L_deposit_h(a[p+1-i]),2);     L_data_move();
  651. L_temp = L_sub(L_ai,L_api);       /* L_temp in Q26 */
  652.         L_p_cof[i] = L_add(L_temp,L_p_cof[i-1]);      /* L_p_cof in Q26 */
  653. /* temp = add(a[i],a[p+1-i]); */
  654. L_temp = L_add(L_ai,L_api);
  655.         L_q_cof[i] = L_sub(L_temp,L_q_cof[i-1]);      /* L_q_cof in Q26 */
  656.     }
  658.     /* Convert p_cof and q_cof to short */
  659.     for ( i = 0; i <= p; i++ ) {
  660.       p_cof[i] = round(L_p_cof[i]);     /* p_cof in Q10 */
  661.       q_cof[i] = round(L_q_cof[i]);     /* q_cof in Q10 */
  662.     }
  664.     /* Find root frequencies of LSP polynomials */
  665.     lsp_to_freq(p_cof,p_freq,p);
  666.     lsp_to_freq(q_cof,q_freq,p);
  668.     /* Combine frequencies into single array    */
  669.     w[0] = 0;
  670.     j = 1;
  671.     for (i = 0; i < (shr(p,1)); i++) {
  672.         w[j++] = q_freq[i];      data_move();
  673.         w[j++] = p_freq[i];      data_move();
  674.     }
  676.     return((Shortword)0);
  678. } /* LPC_PRED2LSP */
  680. /*                                                              */
  681. /*      Subroutine LSP_TO_FREQ: Calculate Line Spectrum Pair    */
  682. /*      root frequencies from LSP polynomial.                   */
  683. /*                                                              */
  684. /* Q values:
  685.    lsp - Q10
  686.    w - Q15 
  687. */
  688. #define DFTLENGTH       512
  689. #define DFTLENGTH_D2    (DFTLENGTH/2)
  690. #define DFTLENGTH_D4    (DFTLENGTH/4)
  692. static void lsp_to_freq(Shortword lsp[], Shortword w[], Shortword p)
  694. {
  695.     Shortword i;
  696.     Shortword p2, count, prev_less, j;
  697.     static Shortword firstcall = 1;
  698.     Longword mag0, mag1, mag2;
  699.     Shortword s_mag0, s_mag1, s_mag2;
  700.     Shortword *p_cos, *p_endlsp, *p_lsp, *p_begcos, *p_endcos;
  701.     Shortword temp1, temp2;
  702.     Longword L_temp1, L_temp2;
  703.     static Shortword lsp_cos[DFTLENGTH];    /* cosine table */
  704.     static Shortword default_w, default_w0;
  705.     static Shortword cos_tab[129] =
  706.     { 
  707.         32767, 32766, 32758, 32746, 32729, 32706, 32679, 32647, 32610,
  708.         32568, 32522, 32470, 32413, 32352, 32286, 32214, 32138, 32058,
  709.         31972, 31881, 31786, 31686, 31581, 31471, 31357, 31238, 31114,
  710.         30986, 30853, 30715, 30572, 30425, 30274, 30118, 29957, 29792,
  711.         29622, 29448, 29269, 29086, 28899, 28707, 28511, 28311, 28106,
  712.         27897, 27684, 27467, 27246, 27020, 26791, 26557, 26320, 26078,
  713.         25833, 25583, 25330, 25073, 24812, 24548, 24279, 24008, 23732,
  714.         23453, 23170, 22884, 22595, 22302, 22006, 21706, 21403, 21097,
  715.         20788, 20475, 20160, 19841, 19520, 19195, 18868, 18538, 18205,
  716.         17869, 17531, 17190, 16846, 16500, 16151, 15800, 15447, 15091,
  717.         14733, 14373, 14010, 13646, 13279, 12910, 12540, 12167, 11793,
  718.         11417, 11039, 10660, 10279,  9896,  9512,  9127,  8740,  8351,
  719.         7962,   7571,  7180,  6787,  6393,  5998,  5602,  5205,  4808,
  720.         4410,   4011,  3612,  3212,  2811,  2411,  2009,  1608,  1206,
  721.         804,     402,     0
  722.     };
  724.     if (firstcall) {
  725.         /* Initialization: fill cosine table */
  726.         firstcall = 0;
  728. /* for (i = 0; i < DFTLENGTH; i++ )
  729.      lsp_cos[i] = cos(i*(TWOPI / DFTLENGTH)); */
  730. for (i = 0; i < DFTLENGTH_D4; i++) {
  731.   lsp_cos[i] = cos_tab[i];
  732.   lsp_cos[i+DFTLENGTH_D4] = -cos_tab[DFTLENGTH_D4-i];
  733.   lsp_cos[i+2*DFTLENGTH_D4] = -cos_tab[i];
  734.   lsp_cos[i+3*DFTLENGTH_D4] = cos_tab[DFTLENGTH_D4-i];
  735. }
  737. /* compute default values for w[] */
  738. /* (1./p2) in Q15 */
  739. default_w = divide_s(ONE_Q11,shl(p,10));
  740. /* w[0] = (0.5/p2) in Q15 */
  741. default_w0 = shr(default_w,1);
  742.     }
  744.     prev_less = 1;
  745.     mag0 = (Longword)0x7fffffff;      L_data_move();
  746.     mag1 = (Longword)0x7fffffff;      L_data_move();
  747.     s_mag0 = (Shortword)0x7fff;      data_move();
  748.     s_mag1 = (Shortword)0x7fff;      data_move();
  749.     p2 = shr(p,1);      /* p/2 */
  750.     count = 0;      data_move();
  751.     p_begcos = &lsp_cos[0];
  752.     p_endcos = &lsp_cos[DFTLENGTH-1];
  753.     p_endlsp = &lsp[p2];
  755.     /*  Search all frequencies for minima of Pc(w) */
  756.     for (i = 0; i <= DFTLENGTH_D2; i++ ) {
  757.         p_lsp = p_endlsp;
  758.         p_cos = p_begcos++;
  759. /* mag2 = 0.5 * *(p_lsp--); */
  760. mag2 = L_mult(*(p_lsp--),(Shortword)X05_Q14);     /* mag2 in Q25 */
  761.         for (j = p2 - 1; j >= 0; j-- ) {
  762.     L_temp1 = L_shr(L_mult(*(p_lsp--),*(p_cos)),1);
  763.             mag2 = L_add(mag2,L_temp1);
  764.     complexity -= 4;    /* this can be done with MAC on TMS320C5x */
  765.             p_cos += i;
  766.             if (p_cos > p_endcos)
  767.                 p_cos -= DFTLENGTH;
  768.         }
  769. s_mag2 = extract_h(mag2);
  770.         mag2 = L_abs(mag2);
  772. L_compare_nonzero();
  773.         if (mag2 < mag1) {
  774.     prev_less = 1;      data_move();
  775. }
  776.         else {
  777.             if (prev_less) {
  778.       /* check for sign change */
  779.       logic();    compare_zero();
  780.       if ((s_mag0 ^ s_mag2) < 0) {
  781. /* Minimum frequency found */
  782. /*w[count] = i - 1 + (0.5 *
  783.              (mag0 - mag2) / (mag0 + mag2 - 2*mag1) );
  784.      w[count] *= (2. / DFTLENGTH) ;*/
  785. L_temp1 = L_shr(L_sub(mag0,mag2),1);
  786. L_temp2 = L_sub(mag0,L_shl(mag1,1));
  787. L_temp2 = L_add(L_temp2,mag2);
  788. temp1 = L_divider2(L_temp1,L_temp2,0,0);  /* temp1 in Q15 */
  789. temp1 = shr(temp1,9);    /* Q6 */
  790. temp2 = sub(i,1);
  791. temp1 = add(shl(temp2,6),temp1);
  792. w[count] = divide_s(temp1,shl(DFTLENGTH,5));      data_move();
  793. count = add(count,1);
  794.       }
  795.     }
  796.     prev_less = 0;      data_move();
  797.         }
  798.         mag0 = mag1;      L_data_move();
  799.         mag1 = mag2;      L_data_move();
  800. s_mag0 = s_mag1;      data_move();
  801. s_mag1 = s_mag2;      data_move();
  803.     }
  804.     
  805.     /*  Verify that all roots were found                */
  806.     if (count != p2) {
  807. #if (PRINT)
  808.         printf("ERROR: couldn't find LSP frequencies in frame %d.n",frames);
  809.         printf(" Found were ");
  810.         for (i = 0; i < count; i++) 
  811.             printf(" %10.4f ",(float)w[i]/(1<<15));
  812.         printf("n");
  813. #endif
  814.         /* use default values */
  815.         w[0] = default_w0;
  816.         for (i = 1; i < p2; i++) {
  817.             w[i] = w[i-1] + default_w;
  818.         }
  819.     }
  820. } /* LSP_TO_FREQ */
  821. #undef  DFTLENGTH
  822. #undef  DFTLENGTH_D2
  823. #undef  DFTLENGTH_D4
  826. /* LPC_PRED2REFL
  827.       get refl coeffs from the predictor coeffs
  828.       Input:
  829.          a- the predictor coefficients
  830.          p- the predictor order
  831.       Output:
  832.          k- the reflection coefficients
  833.       Returns:
  834.          energy - energy of residual signal
  835.    Reference:  Markel and Gray, Linear Prediction of Speech
  837.    Q values:
  838.    *a - Q12
  839.    *k - Q15
  840.    energy - Q15
  841. */
  843. Shortword lpc_pred2refl(Shortword *a,Shortword *k,Shortword p)
  844. {
  845.     Shortword *b,*b1, e;
  846.     Shortword energy = ONE_Q15;
  847.     Shortword i,j;
  848.     Shortword shift;
  849.     Shortword temp, n_temp, n_e;
  851.     MEM_ALLOC(MALLOC,b,p+1,Shortword);
  852.     MEM_ALLOC(MALLOC,b1,p+1,Shortword);
  854.     /* equate temporary variables (b = a) */
  855.     for(i=1; i <= p; i++) {
  856.       b[i] = a[i];      data_move();
  857.     }
  859.     /* compute reflection coefficients */
  860.     for(i=p; i > 0; i--)
  861.     {
  862.         k[i] = shl(b[i],3);     /* Q12 -> Q15 */
  863.         e = sub(ONE_Q15,mult(k[i],k[i]));
  864.         energy = mult(energy,e);
  865.         for(j=1; j < i; j++) {
  866.             b1[j] = b[j];      data_move();
  867. }
  869.         for(j=1; j < i; j++) {
  870.     /* b[j] = (b1[j] - k[i]*b1[i-j])/e; */
  871.     temp = mult(k[i],b1[i-j]);    /* temp in Q12 */
  872.     temp = sub(b1[j],temp);
  874.     /* check signs of temp and e before division */
  875.     n_temp = abs_s(temp);
  876.     n_e = abs_s(e);
  877.     shift = norm_s(n_e);
  878.     n_e = shl(n_e,shift);
  879.     if (n_temp > n_e) {
  880.       n_temp = shr(n_temp,1);
  881.       shift = add(shift,1);
  882.     }
  883.     b[j] = shl(divide_s(n_temp,n_e),shift);    /* b[j] in Q12 */
  885.     if ((temp ^ e) < 0)
  886.       b[j] = negate(b[j]);
  887. }
  888.     }
  890.     MEM_FREE(FREE,b1);
  891.     MEM_FREE(FREE,b);
  892.     return(energy);
  893. } /* LPC_PRED2REFL */
  896. /* LPC_LSP2PRED
  897.       get predictor coefficients from the LSPs
  898.    Synopsis: lpc_lsp2pred(w,a,p)
  899.       Input:
  900.          w- the LSPs
  901.          p- the predictor order
  902.       Output:
  903.          a- the predictor coefficients
  904.    Reference:  Kabal and Ramachandran
  905. */
  906. /* Q values:
  907.    w - Q15
  908.    a - Q12
  909.    f - Q25
  910.    c - Q14 */
  912. Shortword lpc_lsp2pred(Shortword *w,Shortword *a,Shortword p)
  913. {
  914.     Shortword i,j,k,p2;
  915.     Shortword c[2];
  916.     Longword **f;
  917.     Longword L_temp;
  919.     /* ensure minimum separation and sort */
  920.     (void)lpc_clmp(w,lsp_delta,p);
  922.     /* p2 = p/2 */
  923.     p2 = shr(p,1);     data_move();
  925.     MEM_2ALLOC(MALLOC,f,2,p2+1,Longword);
  927.     /* f is Q25 */
  928.     f[0][0] = f[1][0] = (Longword)ONE_Q25;    L_data_move();    L_data_move();
  930.     /* -2.0*cos((double)w[1]*M_PI) */
  931.     f[0][1] = L_shr(L_deposit_h(negate(cos_fxp(w[1]))),5);  L_data_move();
  932.     f[1][1] = L_shr(L_deposit_h(negate(cos_fxp(w[2]))),5);  L_data_move();
  934.     k = 3;
  936.     for(i=2; i <= p2; i++)
  937.       { 
  938.         /* c is Q14 */
  939.         /* multiply by 2 is considered as Q15->Q14 */
  940.         c[0] = negate(cos_fxp(w[k++]));    data_move();
  941.         c[1] = negate(cos_fxp(w[k++]));    data_move();
  942.        
  943.         f[0][i] = f[0][i-2];    L_data_move();
  944.         f[1][i] = f[1][i-2];    L_data_move();
  946.         for(j=i; j >= 2; j--)
  947.           {
  948.     /* f[0][j] += c[0]*f[0][j-1]+f[0][j-2] */
  949.     L_temp = L_mpy_ls(f[0][j-1],c[0]);
  950.     L_temp = L_add(L_shl(L_temp,1),f[0][j-2]);
  951.             f[0][j] = L_add(f[0][j],L_temp);    L_data_move();
  953.             /* f[1][j] += c[1]*f[1][j-1]+f[1][j-2] */
  954.     L_temp = L_mpy_ls(f[1][j-1],c[1]);
  955.     L_temp = L_add(L_shl(L_temp,1),f[1][j-2]);
  956.             f[1][j] = L_add(f[1][j],L_temp);    L_data_move();
  957.           }
  959.         f[0][1] = L_add(f[0][1],L_shl(L_mpy_ls(f[0][0],c[0]),1));
  960.         f[1][1] = L_add(f[1][1],L_shl(L_mpy_ls(f[1][0],c[1]),1));
  961.       }
  963.     for(i=p2; i > 0; i--)
  964.       {
  965.         /* short f (f is a Q14) */
  966.         f[0][i] = L_add(f[0][i],f[0][i-1]);
  967.         f[1][i] = L_sub(f[1][i],f[1][i-1]);
  969.         /* a is Q12 */
  970.         /* a[i] = 0.50*(f[0][i]+f[1][i]) */
  971. /* a[p+1-i] = 0.50*(f[0][i]-f[1][i]) */
  972. /* Q25 -> Q27 -> Q12 */
  973.         a[i] = extract_h(L_shl(L_add(f[0][i],f[1][i]),2));
  974.         a[p+1-i] = extract_h(L_shl(L_sub(f[0][i],f[1][i]),2));
  975.     }
  977.     MEM_2FREE(FREE,f);
  978.     return((Shortword)0);
  979. } /* LPC_LSP2PRED */
  982. /*
  983.     Name: lpc_syn- LPC synthesis filter.
  984.     Aliases: lpc_synthesis
  985.     Description:
  986.         LPC all-pole synthesis filter
  988.         for j = 0 to n-1
  989.             y[j] = x[j] - sum(k=1 to p) y[j-k] a[k]
  991.     Inputs:
  992.         x- input vector (n samples, x[0..n-1])
  993.         a- lpc vector (order p, a[1..p])
  994.         p- order of lpc filter
  995.         n- number of elements in vector which is to be filtered
  996.         y[-p..-1]- filter memory (past outputs)
  997.     Outputs:
  998.         y- output vector (n samples, y[0..n-1])
  999.     Returns: NULL
  1000.     Includes:
  1001.         spbstd.h
  1002.         lpc.h
  1004.     Systems and Info. Science Lab
  1005.     Copyright (c) 1995 by Texas Instruments, Inc.  All rights reserved.
  1007. */
  1009. Shortword lpc_syn(Shortword *x,Shortword *y,Shortword *a,Shortword p,
  1010.   Shortword n)
  1011. {
  1012.     Shortword i,j;
  1013.     Longword accum;
  1015.     for(j=0; j < n; j++) {
  1016.         accum = L_shr(L_deposit_h(x[j]),3);
  1017. for(i=p; i > 0; i--)
  1018.     accum = L_msu(accum,y[j-i],a[i]);
  1019. /* Round off output */
  1020. accum = L_shl(accum,3);
  1021. y[j] = round(accum);
  1022.     }
  1023.     return((Shortword)0);
  1024. }