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

源码开发语言/平台
当前位置：首页 > 源码/资料 > 通讯/手机编程 > Symbian > 查看源码
enc_lpc.c
资源名称：142_61_thumb_advanced.rar　[点击查看]
上传用户：dangjiwu
上传日期：2013-07-19
资源大小：42019k
文件大小：31k
源码类别：

Symbian

开发平台：

Visual C++

enc_lpc.c：源码内容
  1. /*
  2.  *===================================================================
  3.  *  3GPP AMR Wideband Floating-point Speech Codec
  4.  *===================================================================
  5.  */
  6. #include <math.h>
  7. #include "hlxclib/memory.h"
  9. #include "typedef.h"
  10. #include "enc_util.h"
  12. #define ORDER        16          /* order of linear prediction filter   */
  13. #define ISF_GAP      128         /* 50                                  */
  14. #define M            16
  15. #define M16k         20          /* Order of LP filter                  */
  16. #define MP1          (M+1)
  17. #define NC16k        (M16k / 2)
  18. #define MU           10923       /* Prediction factor (1.0/3.0) in Q15  */
  19. #define F_MU         (1.0 / 3.0) /* prediction factor                   */
  20. #define N_SURV_MAX   4           /* 4 survivors max                     */
  22. #ifndef PI
  23. #define PI           3.141592654
  24. #endif
  26. /* isp_isf_conversion */
  27. #define SCALE1  (6400.0/PI)
  29. /* chebyshev */
  30. #define  NO_ITER   4    /* no of iterations for tracking the root */
  31. #define  NO_POINTS 100
  33. #define SIZE_BK1        256
  34. #define SIZE_BK2        256
  35. #define SIZE_BK21       64
  36. #define SIZE_BK22       128
  37. #define SIZE_BK23       128
  38. #define SIZE_BK24       32
  39. #define SIZE_BK25       32
  40. #define SIZE_BK21_36b   128
  41. #define SIZE_BK22_36b   128
  42. #define SIZE_BK23_36b   64
  43. #define SIZE_BK_NOISE1  64
  44. #define SIZE_BK_NOISE2  64
  45. #define SIZE_BK_NOISE3  64
  46. #define SIZE_BK_NOISE4  32
  47. #define SIZE_BK_NOISE5  32
  49. extern const Word16 E_ROM_mean_isf[];
  50. extern const Word16 E_ROM_cos[];
  52. extern const Float32 E_ROM_dico1_isf[];
  53. extern const Float32 E_ROM_dico2_isf[];
  54. extern const Float32 E_ROM_dico21_isf[];
  55. extern const Float32 E_ROM_dico22_isf[];
  56. extern const Float32 E_ROM_dico23_isf[];
  57. extern const Float32 E_ROM_dico24_isf[];
  58. extern const Float32 E_ROM_dico25_isf[];
  59. extern const Float32 E_ROM_dico21_isf_36b[];
  60. extern const Float32 E_ROM_dico22_isf_36b[];
  61. extern const Float32 E_ROM_dico23_isf_36b[];
  62. extern const Float32 E_ROM_f_mean_isf[];
  63. extern const Float32 E_ROM_lag_window[];
  64. extern const Float32 E_ROM_grid[];
  65. extern const Float32 E_ROM_f_interpol_frac[];
  67. /*
  68.  * E_LPC_isf_reorder
  69.  *
  70.  * Parameters:
  71.  *    isf          I/O: vector of isfs
  72.  *    min_dist       I: quantized ISFs (in frequency domain)
  73.  *    n              I: LPC order
  74.  *
  75.  * Function:
  76.  *    To make sure that the  isfs are properly order and to keep a certain
  77.  *    minimum distance between consecutive isfs.
  78.  *
  79.  * Returns:
  80.  *    void
  81.  */
  82. static void E_LPC_isf_reorder(Word16 *isf, Word16 min_dist, Word16 n)
  83. {
  85.    Word32 i, isf_min;
  87.    isf_min = min_dist;
  89.    for (i = 0; i < n - 1; i++)
  90.    {
  91.       if (isf[i] < isf_min)
  92.       {
  93.          isf[i] = (Word16)isf_min;
  94.       }
  96.       isf_min = isf[i] + min_dist;
  97.    }
  99.    return;
  100. }
  102. /*
  103.  * E_LPC_isp_pol_get
  104.  *
  105.  * Parameters:
  106.  *    isp            I: Immitance spectral pairs (cosine domaine)
  107.  *    f              O: the coefficients of F1 or F2
  108.  *    n              I: no of coefficients (m/2)
  109.  *    k16            I: 16k flag
  110.  *
  111.  * Function:
  112.  *    Find the polynomial F1(z) or F2(z) from the ISPs.
  113.  *    This is performed by expanding the product polynomials:
  114.  *
  115.  *    F1(z) =   product   ( 1 - 2 isp_i z^-1 + z^-2 )
  116.  *            i=0,2,4,6,8
  117.  *    F2(z) =   product   ( 1 - 2 isp_i z^-1 + z^-2 )
  118.  *             i=1,3,5,7
  119.  *
  120.  *    where isp_i are the ISPs in the cosine domain.
  121.  *
  122.  * Returns:
  123.  *    void
  124.  */
  125. static void E_LPC_isp_pol_get(Word16 *isp, Word32 *f, Word32 n, Word16 k16)
  126. {
  127.    Word32 i, j, t0, s1, s2;
  128.    Word16 hi, lo;
  130.    s1 = 8388608;
  131.    s2 = 512;
  133.    if(k16)
  134.    {
  135.       s1 >>= 2;
  136.       s2 >>= 2;
  137.    }
  139.    /* All computation in Q23 */
  140.    f[0] = s1;              /* f[0] = 1.0; in Q23         */
  141.    f[1] = isp[0] * (-s2);  /* f[1] = -2.0*isp[0] in Q23  */
  142.    f += 2;                 /* Advance f pointer          */
  143.    isp += 2;               /* Advance isp pointer        */
  145.    for(i = 2; i <= n; i++)
  146.    {
  147.       *f = f[ - 2];
  149.       for(j = 1; j < i; j++, f--)
  150.       {
  151.          E_UTIL_l_extract(f[- 1], &hi, &lo);
  152.          t0 = E_UTIL_mpy_32_16(hi, lo, *isp);   /* t0 = f[-1] * isp */
  153.          t0 = (t0 << 1);
  154.          *f = (*f - t0);         /* *f -= t0    */
  155.          *f = (*f + f[ - 2]);    /* *f += f[-2] */
  156.       }
  158.       *f = *f - (*isp * s2);     /* *f -= isp << 8 */
  159.       f += i;     /* Advance f pointer   */
  160.       isp += 2;   /* Advance isp pointer */
  161.    }
  163.    return;
  164. }
  166. static void E_LPC_f_isp_pol_get(Float32 isp[], Float32 f[], Word32 n)
  167. {
  168.    Float32 b;
  169.    Word32 i, j;
  171.    f[0] = 1;
  172.    b = (Float32)(-2.0 * *isp);
  173.    f[1] = b;
  175.    for (i = 2; i <= n; i++)
  176.    {
  177.       isp += 2;
  178.       b = (Float32)(-2.0 * *isp);
  179.       f[i] = (Float32)(b * f[i - 1] + 2.0 * f[i - 2]);
  181.       for (j = i - 1; j > 1; j--)
  182.       {
  183.          f[j] += b * f[j - 1] + f[j - 2];
  184.       }
  186.       f[1] += b;
  187.    }
  189.    return;
  190. }
  192. /*
  193.  * E_LPC_isp_a_conversion
  194.  *
  195.  * Parameters:
  196.  *    isp            I: (Q15) Immittance spectral pairs
  197.  *    a              O: (Q12) Predictor coefficients (order = M)
  198.  *    m              I: order of LP filter
  199.  *
  200.  * Function:
  201.  *    Convert ISPs to predictor coefficients a[]
  202.  *
  203.  * Returns:
  204.  *    void
  205.  */
  206. void E_LPC_isp_a_conversion(Word16 isp[], Word16 a[], Word16 m)
  207. {
  208.    Word32 f1[NC16k + 1], f2[NC16k];
  209.    Word32 i, j, nc, t0;
  210.    Word16 hi, lo;
  212.    nc = m >> 1;
  214.    if (nc > 8)
  215.    {
  216.       E_LPC_isp_pol_get(&isp[0], f1, nc, 1);
  218.       for (i = 0; i <= nc; i++)
  219.       {
  220.          f1[i] = (f1[i] << 2);
  221.       }
  222.    }
  223.    else
  224.    {
  225.       E_LPC_isp_pol_get(&isp[0], f1, nc, 0);
  226.    }
  228.    if (nc > 8)
  229.    {
  230.       E_LPC_isp_pol_get(&isp[1], f2, nc - 1, 1);
  232.       for (i = 0; i <= nc - 1; i++)
  233.       {
  234.          f2[i] = (f2[i] << 2);
  235.       }
  236.    }
  237.    else
  238.    {
  239.       E_LPC_isp_pol_get(&isp[1], f2, nc - 1, 0);
  240.    }
  242.    /* Multiply F2(z) by (1 - z^-2)  */
  243.    for (i = (nc - 1); i > 1; i--)
  244.    {
  245.       f2[i] = f2[i] - f2[i - 2];     /* f2[i] -= f2[i-2]; */
  246.    }
  248.    /*  Scale F1(z) by (1+isp[m-1])  and  F2(z) by (1-isp[m-1]) */
  249.    for (i = 0; i < nc; i++)
  250.    {
  251.       /* f1[i] *= (1.0 + isp[M-1]); */
  253.       E_UTIL_l_extract(f1[i], &hi, &lo);
  254.       t0 = E_UTIL_mpy_32_16(hi, lo, isp[m - 1]);
  255.       f1[i] = f1[i] + t0;
  257.       /* f2[i] *= (1.0 - isp[M-1]); */
  259.       E_UTIL_l_extract(f2[i], &hi, &lo);
  260.       t0 = E_UTIL_mpy_32_16(hi, lo, isp[m - 1]);
  261.       f2[i] = f2[i] - t0;
  262.    }
  264.    /*
  265.     * A(z) = (F1(z)+F2(z))/2
  266.     * F1(z) is symmetric and F2(z) is antisymmetric
  267.     */
  269.    /* a[0] = 1.0; */
  270.    a[0] = 4096;
  272.    for (i = 1, j = (m - 1); i < nc; i++, j--)
  273.    {
  274.       /* a[i] = 0.5*(f1[i] + f2[i]); */
  275.       t0 = f1[i] + f2[i];                    /* f1[i] + f2[i]             */
  276.       a[i] = (Word16)((t0 + 0x800) >> 12);   /* from Q23 to Q12 and * 0.5 */
  278.       /* a[j] = 0.5*(f1[i] - f2[i]); */
  279.       t0 = (f1[i] - f2[i]);                  /* f1[i] - f2[i]             */
  280.       a[j] = (Word16)((t0 + 0x800) >> 12);   /* from Q23 to Q12 and * 0.5 */
  282.    }
  284.    /* a[NC] = 0.5*f1[NC]*(1.0 + isp[M-1]); */
  285.    E_UTIL_l_extract(f1[nc], &hi, &lo);
  286.    t0 = E_UTIL_mpy_32_16(hi, lo, isp[m - 1]);
  287.    t0 = (f1[nc] + t0);
  288.    a[nc] = (Word16)((t0 + 0x800) >> 12);    /* from Q23 to Q12 and * 0.5 */
  290.    /* a[m] = isp[m-1]; */
  291.    a[m] = (Word16)((isp[m - 1] + 0x4) >> 3); /* from Q15 to Q12          */
  293.    return;
  294. }
  296. void E_LPC_f_isp_a_conversion(Float32 *isp, Float32 *a, Word32 m)
  297. {
  298.    Float32 f1[(M16k / 2) + 1], f2[M16k / 2];
  299.    Word32 i, j, nc;
  301.    nc = m / 2;
  303.    /*
  304.     *  Find the polynomials F1(z) and F2(z)
  305.     */
  307.    E_LPC_f_isp_pol_get(&isp[0], f1, nc);
  308.    E_LPC_f_isp_pol_get(&isp[1], f2, nc-1);
  310.    /*
  311.     *  Multiply F2(z) by (1 - z^-2)
  312.     */
  313.    for (i = (nc - 1); i > 1; i--)
  314.    {
  315.       f2[i] -= f2[i - 2];
  316.    }
  318.    /*
  319.     *  Scale F1(z) by (1+isp[m-1])  and  F2(z) by (1-isp[m-1])
  320.     */
  322.    for (i = 0; i < nc; i++)
  323.    {
  324.       f1[i] *= (Float32)(1.0 + isp[m - 1]);
  325.       f2[i] *= (Float32)(1.0 - isp[m - 1]);
  326.    }
  328.    /*
  329.     *  A(z) = (F1(z)+F2(z))/2
  330.     *  F1(z) is symmetric and F2(z) is antisymmetric
  331.     */
  333.    a[0] = 1.0;
  335.    for (i = 1, j = m - 1; i < nc; i++, j--)
  336.    {
  337.       a[i] = (Float32)(0.5 * (f1[i] + f2[i]));
  338.       a[j] = (Float32)(0.5 * (f1[i] - f2[i]));
  339.    }
  341.    a[nc] = (Float32)(0.5 * f1[nc] * (1.0 + isp[m - 1]));
  342.    a[m] = isp[m - 1];
  344.    return;
  345. }
  347. /*
  348.  * E_LPC_int_isp_find
  349.  *
  350.  * Parameters:
  351.  *    isp_old           I: isps from past frame
  352.  *    isp_new           I: isps from present frame
  353.  *    frac              I: (Q15) fraction for 3 first subfr
  354.  *    Az                O: LP coefficients in 4 subframes
  355.  *
  356.  * Function:
  357.  *    Find the Word32erpolated ISP parameters for all subframes.
  358.  *
  359.  * Returns:
  360.  *    void
  361.  */
  362. void E_LPC_int_isp_find(Word16 isp_old[], Word16 isp_new[],
  363.                         const Word16 frac[], Word16 Az[])
  364. {
  365.    Word32 i, k, fac_old, fac_new, tmp;
  366.    Word16 isp[M];
  368.    for (k = 0; k < 3; k++)
  369.    {
  370.       fac_new = frac[k];
  371.       fac_old = ((32767 - fac_new) + 1);  /* 1.0 - fac_new */
  373.       for (i = 0; i < M; i++)
  374.       {
  375.          tmp = isp_old[i] * fac_old;
  376.          tmp += (isp_new[i] * fac_new);
  377.          isp[i] = (Word16)((tmp + 0x4000) >> 15);
  378.       }
  380.       E_LPC_isp_a_conversion(isp, Az, M);
  381.       Az += MP1;
  382.    }
  384.    /* 4th subframe: isp_new (frac=1.0) */
  385.    E_LPC_isp_a_conversion(isp_new, Az, M);
  387.    return;
  388. }
  390. void E_LPC_f_int_isp_find(Float32 isp_old[], Float32 isp_new[], Float32 a[],
  391.                           Word32 nb_subfr, Word32 m)
  392. {
  393.   Float32 isp[M], fnew, fold;
  394.   Float32 *p_a;
  395.   Word32 i, k;
  397.   p_a = a;
  399.   for (k = 0; k < nb_subfr; k++)
  400.   {
  401.     fnew = E_ROM_f_interpol_frac[k];
  402.     fold = (Float32)(1.0 - fnew);
  404.     for (i = 0; i < m; i++)
  405.     {
  406.       isp[i] = isp_old[i] * fold + isp_new[i] * fnew;
  407.     }
  409.     E_LPC_f_isp_a_conversion(isp, p_a, m);
  410.     p_a += (m + 1);
  411.   }
  413.   return;
  414. }
  416. /*
  417.  * E_LPC_a_weight
  418.  *
  419.  * Parameters:
  420.  *    a              I: LP filter coefficients
  421.  *    ap             O: weighted LP filter coefficients
  422.  *    gamma          I: weighting factor
  423.  *    m              I: order of LP filter
  424.  *
  425.  * Function:
  426.  *    Weighting of LP filter coefficients, ap[i] = a[i] * (gamma^i).
  427.  *
  428.  * Returns:
  429.  *    void
  430.  */
  431. void E_LPC_a_weight(Float32 *a, Float32 *ap, Float32 gamma, Word32 m)
  432. {
  433.    Float32 f;
  434.    Word32 i;
  436.    ap[0] = a[0];
  437.    f = gamma;
  439.    for (i = 1; i <= m; i++)
  440.    {
  441.       ap[i] = f*a[i];
  442.       f *= gamma;
  443.    }
  445.    return;
  446. }
  448. /*
  449.  * E_LPC_isf_2s3s_decode
  450.  *
  451.  * Parameters:
  452.  *    indice            I: quantisation indices
  453.  *    isf_q             O: quantised ISFs in the cosine domain
  454.  *    past_isfq       I/O: past ISF quantizer
  455.  *
  456.  * Function:
  457.  *    Decoding of ISF parameters.
  458.  *
  459.  * Returns:
  460.  *    void
  461.  */
  462. static void E_LPC_isf_2s3s_decode(Word32 *indice, Word16 *isf_q,
  463.                                   Word16 *past_isfq)
  464. {
  465.    Word32 i;
  466.    Word16 tmp;
  468.    for(i = 0; i < 9; i++)
  469.    {
  470.       isf_q[i] = (Word16)((E_ROM_dico1_isf[indice[0] * 9 + i] * 2.56F) + 0.5F);
  471.    }
  473.    for(i = 0; i < 7; i++)
  474.    {
  475.       isf_q[i + 9] =
  476.          (Word16)((E_ROM_dico2_isf[indice[1] * 7 + i] * 2.56F) + 0.5F);
  477.    }
  479.    for(i = 0; i < 5; i++)
  480.    {
  481.       isf_q[i] = (Word16)(isf_q[i] +
  482.          (Word16)((E_ROM_dico21_isf_36b[indice[2] * 5 + i] * 2.56F) + 0.5F));
  483.    }
  485.    for(i = 0; i < 4; i++)
  486.    {
  487.       isf_q[i + 5] = (Word16)(isf_q[i + 5] +
  488.          (Word32)((E_ROM_dico22_isf_36b[indice[3] * 4 + i] * 2.56F) + 0.5F));
  489.    }
  491.    for(i = 0; i < 7; i++)
  492.    {
  493.       isf_q[i + 9] = (Word16)(isf_q[i + 9] +
  494.          (Word32)((E_ROM_dico23_isf_36b[indice[4] * 7 + i] * 2.56F) + 0.5F));
  495.    }
  497.    for(i = 0; i < ORDER; i++)
  498.    {
  499.       tmp = isf_q[i];
  500.       isf_q[i] = (Word16)(tmp + E_ROM_mean_isf[i]);
  501.       isf_q[i] = (Word16)(isf_q[i] + ((MU * past_isfq[i]) >> 15));
  502.       past_isfq[i] = tmp;
  503.    }
  505.    E_LPC_isf_reorder(isf_q, ISF_GAP, ORDER);
  507.    return;
  508. }
  511. /*
  512.  * E_LPC_isf_2s5s_decode
  513.  *
  514.  * Parameters:
  515.  *    indice            I: quantization indices
  516.  *    isf_q             O: quantized ISFs in the cosine domain
  517.  *    past_isfq       I/O: past ISF quantizer
  518.  *    isfold            I: past quantized ISF
  519.  *    isf_buf           I: isf buffer
  520.  *    bfi               I: Bad frame indicator
  521.  *    enc_dec           I:
  522.  *
  523.  * Function:
  524.  *    Decoding of ISF parameters.
  525.  *
  526.  * Returns:
  527.  *    void
  528.  */
  529. void E_LPC_isf_2s5s_decode(Word32 *indice, Word16 *isf_q, Word16 *past_isfq)
  530. {
  531.    Word32 i;
  532.    Word16 tmp;
  534.    for (i = 0; i < 9; i++)
  535.    {
  536.       isf_q[i] = (Word16)((E_ROM_dico1_isf[indice[0] * 9 + i] * 2.56F) + 0.5F);
  537.    }
  539.    for (i = 0; i < 7; i++)
  540.    {
  541.       isf_q[i + 9] =
  542.          (Word16)((E_ROM_dico2_isf[indice[1] * 7 + i] * 2.56F) + 0.5F);
  543.    }
  545.    for (i = 0; i < 3; i++)
  546.    {
  547.       isf_q[i] = (Word16)(isf_q[i] +
  548.          (Word32)((E_ROM_dico21_isf[indice[2] * 3 + i] * 2.56F) + 0.5F));
  549.    }
  551.    for (i = 0; i < 3; i++)
  552.    {
  553.       isf_q[i + 3] = (Word16)(isf_q[i + 3] +
  554.          (Word32)((E_ROM_dico22_isf[indice[3] * 3 + i] * 2.56F) + 0.5F));
  555.    }
  557.    for (i = 0; i < 3; i++)
  558.    {
  559.       isf_q[i + 6] = (Word16)(isf_q[i + 6] +
  560.          (Word32)((E_ROM_dico23_isf[indice[4] * 3 + i] * 2.56F) + 0.5F));
  561.    }
  563.    for (i = 0; i < 3; i++)
  564.    {
  565.       isf_q[i + 9] = (Word16)(isf_q[i + 9] +
  566.          (Word32)((E_ROM_dico24_isf[indice[5] * 3 + i] * 2.56F) + 0.5F));
  567.    }
  569.    for (i = 0; i < 4; i++)
  570.    {
  571.       isf_q[i + 12] = (Word16)(isf_q[i + 12] +
  572.          (Word32)((E_ROM_dico25_isf[indice[6] * 4 + i] * 2.56F) + 0.5F));
  573.    }
  575.    for (i = 0; i < ORDER; i++)
  576.    {
  577.       tmp = isf_q[i];
  578.       isf_q[i] = (Word16)(tmp + E_ROM_mean_isf[i]);
  579.       isf_q[i] = (Word16)(isf_q[i] + ((MU * past_isfq[i]) >> 15));
  580.       past_isfq[i] = tmp;
  581.    }
  583.    E_LPC_isf_reorder(isf_q, ISF_GAP, ORDER);
  585.    return;
  586. }
  588. /*
  589.  * E_LPC_isf_isp_conversion
  590.  *
  591.  * Parameters:
  592.  *    isp            O: (Q15) isp[m] (range: -1<=val<1)
  593.  *    isf            I: (Q15) isf[m] normalized (range: 0.0 <= val <= 0.5)
  594.  *    m              I: LPC order
  595.  *
  596.  * Function:
  597.  *    Transformation isf to isp
  598.  *
  599.  *    ISP are immitance spectral pair in cosine domain (-1 to 1).
  600.  *    ISF are immitance spectral pair in frequency domain (0 to 6400).
  601.  * Returns:
  602.  *    void
  603.  */
  604. void E_LPC_isf_isp_conversion(Word16 isf[], Word16 isp[], Word16 m)
  605. {
  606.    Word32 i, ind, offset, tmp;
  608.    for (i = 0; i < m - 1; i++)
  609.    {
  610.       isp[i] = isf[i];
  611.    }
  613.    isp[m - 1] = (Word16)(isf[m - 1] << 1);
  615.    for (i = 0; i < m; i++)
  616.    {
  617.       ind = isp[i] >> 7;         /* ind    = b7-b15 of isf[i] */
  618.       offset = isp[i] & 0x007f;  /* offset = b0-b6  of isf[i] */
  620.       /* isp[i] = table[ind]+ ((table[ind+1]-table[ind])*offset) / 128 */
  621.       tmp = ((E_ROM_cos[ind + 1] - E_ROM_cos[ind]) * offset) << 1;
  622.       isp[i] = (Word16)(E_ROM_cos[ind] + (tmp >> 8));
  623.    }
  625.    return;
  626. }
  628. /*
  629.  * E_LPC_chebyshev
  630.  *
  631.  * Parameters:
  632.  *    x           I: value of evaluation; x=cos(freq)
  633.  *    f           I: coefficients of sum or diff polynomial
  634.  *    n           I: order of polynomial
  635.  *
  636.  * Function:
  637.  *    Evaluates the Chebyshev polynomial series
  638.  *
  639.  *    The polynomial order is n = m/2   (m is the prediction order)
  640.  *    The polynomial is given by
  641.  *    C(x) = f(0)T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2
  642.  *
  643.  * Returns:
  644.  *    the value of the polynomial C(x)
  645.  */
  646. static Float32 E_LPC_chebyshev(Float32 x, Float32 *f, Word32 n)
  647. {
  648.    Float32 b1, b2, b0, x2;
  649.    Word32 i;                        /* for the special case of 10th order  */
  650.                                     /*       filter (n=5)                  */
  651.    x2 = 2.0F * x;                   /* x2 = 2.0*x;                         */
  652.    b2 = f[0];                       /* b2 = f[0];                          */
  653.    b1 = x2 * b2 + f[1];             /* b1 = x2*b2 + f[1];                  */
  655.    for (i = 2; i < n; i++)
  656.    {
  657.       b0 = x2 * b1 - b2 + f[i];     /* b0 = x2 * b1 - 1. + f[2];           */
  658.       b2 = b1;                      /* b2 = x2 * b0 - b1 + f[3];           */
  659.       b1 = b0;                      /* b1 = x2 * b2 - b0 + f[4];           */
  660.    }
  662.    return (x * b1 - b2 + 0.5F * f[n]); /* return (x*b1 - b2 + 0.5*f[5]);   */
  663. }
  664. /*
  665.  * E_LPC_isf_sub_vq
  666.  *
  667.  * Parameters:
  668.  *    x            I/O: unquantised / quantised ISF
  669.  *    dico           I: codebook
  670.  *    dim            I: dimension of vector
  671.  *    dico_size      I: codebook size
  672.  *    distance       O: minimum distance
  673.  *
  674.  * Function:
  675.  *    Quantization of a subvector of size 2 in Split-VQ of ISFs
  676.  *
  677.  * Returns:
  678.  *    quantisation index
  679.  */
  680. Word16 E_LPC_isf_sub_vq(Float32 *x, const Float32 *E_ROM_dico, Word32 dim,
  681.                         Word32 E_ROM_dico_size, Float32 *distance)
  682. {
  683.    Float32 dist_min, dist, temp;
  684.    const Float32 *p_E_ROM_dico;
  685.    Word32 i, j, index = 0;
  687.    dist_min = 1.0e30f;
  688.    p_E_ROM_dico = E_ROM_dico;
  690.    for (i = 0; i < E_ROM_dico_size; i++)
  691.    {
  692.       dist = x[0] - *p_E_ROM_dico++;
  693.       dist *= dist;
  695.       for (j = 1; j < dim; j++)
  696.       {
  697.          temp = x[j] - *p_E_ROM_dico++;
  698.          dist += temp * temp;
  699.       }
  701.       if (dist < dist_min)
  702.       {
  703.          dist_min = dist;
  704.          index = i;
  705.       }
  706.    }
  708.    *distance = dist_min;
  710.    /* Reading the selected vector */
  711.    memcpy(x, &E_ROM_dico[index * dim], dim * sizeof(Float32));
  713.    return (Word16)index;
  714. }
  716. /*
  717.  * E_LPC_lag_wind
  718.  *
  719.  * Parameters:
  720.  *    r         I/O: autocorrelations vector
  721.  *    m           I: window lenght
  722.  *
  723.  * Function:
  724.  *    Lag windowing of the autocorrelations.
  725.  *
  726.  * Returns:
  727.  *    void
  728.  */
  729. void E_LPC_lag_wind(Float32 r[], Word32 m)
  730. {
  731.    Word32 i;
  733.    for (i = 0; i < m; i++)
  734.    {
  735.       r[i] *= E_ROM_lag_window[i];
  736.    }
  738.    return;
  739. }
  741. /*
  742.  * E_LPC_lev_dur
  743.  *
  744.  * Parameters:
  745.  *    r_h         I: vector of autocorrelations (msb)
  746.  *    r_l         I: vector of autocorrelations (lsb)
  747.  *    A           O: LP coefficients (a[0] = 1.0) (m = 16)
  748. *     rc          O: reflection coefficients
  749.  *    mem       I/O: static memory
  750.  *
  751.  * Function:
  752.  *    Wiener-Levinson-Durbin algorithm to compute
  753.  *    the LPC parameters from the autocorrelations of speech.
  754.  *
  755.  * Returns:
  756.  *    void
  757.  */
  758. void E_LPC_lev_dur(Float32 *a, Float32 *r, Word32 m)
  759. {
  760.    Float32 buf[M];
  761.    Float32 *rc;    /* reflection coefficients  0,...,m-1 */
  762.    Float32 s, at, err;
  763.    Word32 i, j, l;
  765.    rc = &buf[0];
  767.    rc[0] = (-r[1]) / r[0];
  768.    a[0] = 1.0F;
  769.    a[1] = rc[0];
  770.    err = r[0] + r[1] * rc[0];
  772.    for (i = 2; i <= m; i++)
  773.    {
  774.       s = 0.0F;
  776.       for (j = 0; j < i; j++)
  777.       {
  778.          s += r[i - j] * a[j];
  779.       }
  781.       rc[i - 1]= (-s) / (err);
  783.       for (j = 1; j <= (i >> 1); j++)
  784.       {
  785.          l = i - j;
  786.          at = a[j] + rc[i - 1] * a[l];
  787.          a[l] += rc[i - 1] * a[j];
  788.          a[j] = at;
  789.       }
  791.       a[i] = rc[i - 1];
  792.       err += rc[i - 1] * s;
  794.       if (err <= 0.0F)
  795.       {
  796.          err = 0.01F;
  797.       }
  798.    }
  800.    return;
  801. }
  803. /*
  804.  * E_LPC_a_isp_conversion
  805.  *
  806.  * Parameters:
  807.  *    a           I: LP filter coefficients
  808.  *    isp         O: Immittance spectral pairs
  809.  *    old_isp     I: ISP vector from past frame
  810.  *
  811.  * Function:
  812.  *    Compute the ISPs from the LPC coefficients a[] using Chebyshev
  813.  *    polynomials. The found ISPs are in the cosine domain with values
  814.  *    in the range from 1 down to -1.
  815.  *    The table E_ROM_grid[] contains the polongs (in the cosine domain) at
  816.  *    which the polynomials are evaluated.
  817.  *
  818.  *    The ISPs are the roots of the two polynomials F1(z) and F2(z)
  819.  *    defined as
  820.  *               F1(z) = A(z) + z^-m A(z^-1)
  821.  *    and        F2(z) = A(z) - z^-m A(z^-1)
  822.  *
  823.  *    for a even order m=2n, F1(z) has 5 conjugate roots on the unit circle
  824.  *    and F2(z) has 4 conjugate roots on the unit circle in addition to two
  825.  *    roots at 0 and pi.
  826.  *
  827.  *    For a 10th order LP analysis, F1(z) and F2(z) can be written as
  828.  *
  829.  *    F1(z) = (1 + a[10])   PRODUCT  (1 - 2 cos(w_i) z^-1 + z^-2 )
  830.  *                       i=0,2,4,6,8
  831.  *
  832.  *    F2(z) = (1 - a[10]) (1 - z^-2) PRODUCT  (1 - 2 cos(w_i) z^-1 + z^-2 )
  833.  *                                  i=1,3,5,7
  834.  *
  835.  *    The ISPs are the M-1 frequencies w_i, i=0...M-2 plus the last
  836.  *    predictor coefficient a[M].
  837.  *
  838.  * Returns:
  839.  *    void
  840.  */
  841. void E_LPC_a_isp_conversion(Float32 *a, Float32 *isp, Float32 *old_isp,
  842.                             Word32 m)
  843. {
  844.    Float32 f1[(M >> 1) + 1], f2[M >> 1];
  845.    Float32 *pf;
  846.    Float32 xlow, ylow, xhigh, yhigh, xmid, ymid, xint;
  847.    Word32 j, i, nf, ip, order, nc;
  849.    nc = m >> 1;
  851.    /*
  852.     * find the sum and diff polynomials F1(z) and F2(z)
  853.     *      F1(z) = [A(z) + z^m A(z^-1)]
  854.     *      F2(z) = [A(z) - z^m A(z^-1)]/(1-z^-2)
  855.     */
  857.    for (i=0; i < nc; i++)
  858.    {
  859.       f1[i] = a[i] + a[m - i];
  860.       f2[i] = a[i] - a[m - i];
  861.    }
  863.    f1[nc] = 2.0F * a[nc];
  865.    /* divide by (1 - z^-2) */
  866.    for (i = 2; i < nc; i++)
  867.    {
  868.       f2[i] += f2[i - 2];
  869.    }
  871.    /*
  872.     * Find the ISPs (roots of F1(z) and F2(z) ) using the
  873.     * Chebyshev polynomial evaluation.
  874.     * The roots of F1(z) and F2(z) are alternatively searched.
  875.     * We start by finding the first root of F1(z) then we switch
  876.     * to F2(z) then back to F1(z) and so on until all roots are found.
  877.     *
  878.     *  - Evaluate Chebyshev pol. at E_ROM_grid polongs and check for sign change.
  879.     *  - If sign change track the root by subdividing the Word32erval
  880.     *    4 times and ckecking sign change.
  881.     */
  883.    nf=0;      /* number of found frequencies */
  884.    ip=0;      /* flag to first polynomial   */
  886.    pf = f1;  /* start with F1(z) */
  887.    order = nc;
  889.    xlow = E_ROM_grid[0];
  890.    ylow = E_LPC_chebyshev(xlow, pf, nc);
  892.    j = 0;
  894.    while ( (nf < m - 1) && (j < NO_POINTS) )
  895.    {
  896.       j++;
  897.       xhigh = xlow;
  898.       yhigh = ylow;
  899.       xlow = E_ROM_grid[j];
  900.       ylow = E_LPC_chebyshev(xlow, pf, order);
  902.       if (ylow * yhigh <= 0.0F)  /* if sign change new root exists */
  903.       {
  904.          j--;
  906.          /* divide the Word32erval of sign change by NO_ITER */
  908.          for (i = 0; i < NO_ITER; i++)
  909.          {
  910.             xmid = 0.5F * (xlow + xhigh);
  911.             ymid = E_LPC_chebyshev(xmid, pf, order);
  913.             if (ylow * ymid <= 0.0F)
  914.             {
  915.                yhigh = ymid;
  916.                xhigh = xmid;
  917.             }
  918.             else
  919.             {
  920.                ylow = ymid;
  921.                xlow = xmid;
  922.             }
  923.          }
  925.          /* linear interpolation for evaluating the root */
  927.          xint = xlow - ylow * (xhigh - xlow) / (yhigh - ylow);
  929.          isp[nf] = xint;    /* new root */
  930.          nf++;
  932.          ip = 1 - ip;                  /* flag to other polynomial    */
  933.          pf = ip ? f2 : f1;            /* pointer to other polynomial */
  934.          order = ip ? (nc - 1) : nc;   /* order of other polynomial   */
  936.          xlow = xint;
  937.          ylow = E_LPC_chebyshev(xlow, pf, order);
  938.       }
  939.    }
  941.    isp[m - 1] = a[m];
  943.    /*
  944.     * Check if m-1 roots found
  945.     * if not use the ISPs from previous frame
  946.     */
  948.    if (nf < m - 1)
  949.    {
  950.       for(i = 0; i < m; i++)
  951.       {
  952.          isp[i] = old_isp[i];
  953.       }
  954.    }
  956.    return;
  957. }
  959. /*
  960.  * E_LPC_isp_isf_conversion
  961.  *
  962.  * Parameters:
  963.  *    isp            I: isp[m] (range: -1 <= val < 1) (Q15)
  964.  *    isf            O: isf[m] normalized (range: 0 <= val <= 6400)
  965.  *    m              I: LPC order
  966.  *
  967.  * Function:
  968.  *    Transformation isp to isf
  969.  *
  970.  *    ISP are immitance spectral pair in cosine domain (-1 to 1).
  971.  *    ISF are immitance spectral pair in frequency domain (0 to 6400).
  972.  * Returns:
  973.  *    energy of prediction error
  974.  */
  975. void E_LPC_isp_isf_conversion(Float32 isp[], Float32 isf[], Word32 m)
  976. {
  977.    Word32 i;
  979.    /* convert ISPs to frequency domain 0..6400 */
  980.    for(i = 0; i < (m - 1); i++)
  981.    {
  982.       isf[i] = (Float32)(acos(isp[i]) * SCALE1);
  983.    }
  985.    isf[m - 1] = (Float32)(acos(isp[m - 1]) * SCALE1 * 0.5F);
  987.    return;
  988. }
  990. /*
  991.  * E_LPC_stage1_isf_vq
  992.  *
  993.  * Parameters:
  994.  *    x              I: ISF residual vector
  995.  *    dico           I: quantisation codebook
  996.  *    dim            I: dimension of vector
  997.  *    dico_size      I: size of quantization codebook
  998.  *    index          O: indices of survivors
  999.  *    surv           I: number of survivor
  1000.  *
  1001.  * Function:
  1002.  *    1st stage VQ with split-by-2.
  1003.  *
  1004.  * Returns:
  1005.  *    void
  1006.  */
  1007. static void E_LPC_stage1_isf_vq(Float32 *x, const Float32 *E_ROM_dico,
  1008.                                 Word32 dim, Word32 E_ROM_dico_size,
  1009.                                 Word32 *index, Word32 surv)
  1010. {
  1011.    Float32 dist_min[N_SURV_MAX];
  1012.    Float32 dist, temp1, temp2;
  1013.    const Float32 *p_E_ROM_dico;
  1014.    Word32 i, j, k, l;
  1016.    for (i = 0; i < surv; i++)
  1017.    {
  1018.       dist_min[i] = 1.0e30F;
  1019.    }
  1021.    for (i = 0; i < surv; i++)
  1022.    {
  1023.       index[i] = i;
  1024.    }
  1026.    p_E_ROM_dico = E_ROM_dico;
  1028.    for (i = 0; i < E_ROM_dico_size; i++)
  1029.    {
  1030.       dist = x[0] - *p_E_ROM_dico++;
  1031.       dist *= dist;
  1033.       for (j = 1; j < dim; j += 2)
  1034.       {
  1035.          temp1 = x[j] - *p_E_ROM_dico++;
  1036.          temp2 = x[j + 1] - *p_E_ROM_dico++;
  1037.          dist += temp1 * temp1 + temp2 * temp2;
  1038.       }
  1040.       for (k = 0; k < surv; k++)
  1041.       {
  1042.          if (dist < dist_min[k])
  1043.          {
  1044.             for (l = surv - 1; l > k; l--)
  1045.             {
  1046.                dist_min[l] = dist_min[l - 1];
  1047.                index[l] = index[l - 1];
  1048.             }
  1050.             dist_min[k] = dist;
  1051.             index[k] = i;
  1053.             break;
  1054.          }
  1055.       }
  1056.    }
  1058.    return;
  1059. }
  1061. /*
  1062.  * E_LPC_isf_2s3s_quantise
  1063.  *
  1064.  * Parameters:
  1065.  *    isf1              I: ISF in the frequency domain (0..6400)
  1066.  *    isf_q             O: quantized ISF
  1067.  *    past_isfq       I/O: past ISF quantizer
  1068.  *    indice            O: quantisation indices (5 words)
  1069.  *    nb_surv           I: number of survivor (1, 2, 3 or 4)
  1070.  *
  1071.  * Function:
  1072.  *    Quantization of isf parameters with prediction. (36 bits)
  1073.  *
  1074.  *    The isf vector is quantized using two-stage VQ with split-by-2 in
  1075.  *    1st stage and split-by-3 in the second stage.
  1076.  * Returns:
  1077.  *    void
  1078.  */
  1079. void E_LPC_isf_2s3s_quantise(Float32 *isf1, Word16 *isf_q, Word16 *past_isfq,
  1080.                              Word32 *indice, Word32 nb_surv)
  1081. {
  1082.    Float32 isf[ORDER], isf_stage2[ORDER];
  1083.    Float32 temp, min_err, distance;
  1084.    Word32 surv1[N_SURV_MAX];     /* indices of survivors from 1st stage */
  1085.    Word32 tmp_ind[5];
  1086.    Word32 i, k;
  1088.    for (i = 0; i < ORDER; i++)
  1089.    {
  1090.       isf[i] = (Float32)((isf1[i] - E_ROM_f_mean_isf[i]) -
  1091.          F_MU * past_isfq[i] * 0.390625F);
  1092.    }
  1094.    E_LPC_stage1_isf_vq(&isf[0], E_ROM_dico1_isf, 9, SIZE_BK1, surv1, nb_surv);
  1096.    distance = 1.0e30F;
  1098.    for (k = 0; k < nb_surv; k++)
  1099.    {
  1100.       for (i = 0; i < 9; i++)
  1101.       {
  1102.          isf_stage2[i] = isf[i] - E_ROM_dico1_isf[i + surv1[k] * 9];
  1103.       }
  1105.       tmp_ind[0] = E_LPC_isf_sub_vq(&isf_stage2[0], E_ROM_dico21_isf_36b, 5,
  1106.          SIZE_BK21_36b, &min_err);
  1107.       temp = min_err;
  1108.       tmp_ind[1] = E_LPC_isf_sub_vq(&isf_stage2[5], E_ROM_dico22_isf_36b, 4,
  1109.          SIZE_BK22_36b, &min_err);
  1110.       temp += min_err;
  1112.       if (temp < distance)
  1113.       {
  1114.          distance = temp;
  1115.          indice[0] = surv1[k];
  1117.          for (i = 0; i < 2; i++)
  1118.          {
  1119.             indice[i + 2] = tmp_ind[i];
  1120.          }
  1121.       }
  1122.    }
  1124.    E_LPC_stage1_isf_vq(&isf[9], E_ROM_dico2_isf, 7, SIZE_BK2, surv1, nb_surv);
  1126.    distance = 1.0e30F;
  1128.    for (k = 0; k < nb_surv; k++)
  1129.    {
  1130.       for (i = 0; i < 7; i++)
  1131.       {
  1132.          isf_stage2[i] = isf[9 + i] - E_ROM_dico2_isf[i + surv1[k] * 7];
  1133.       }
  1135.       tmp_ind[0] = E_LPC_isf_sub_vq(&isf_stage2[0], E_ROM_dico23_isf_36b, 7,
  1136.          SIZE_BK23_36b, &min_err);
  1137.       temp = min_err;
  1139.       if (temp < distance)
  1140.       {
  1141.          distance = temp;
  1142.          indice[1] = surv1[k];
  1143.          indice[4]= tmp_ind[0];
  1144.       }
  1145.    }
  1147.    /* decoding the ISF */
  1149.    E_LPC_isf_2s3s_decode(indice, isf_q, past_isfq);
  1151.    return;
  1152. }
  1154. /*
  1155.  * E_LPC_isf_2s5s_quantise
  1156.  *
  1157.  * Parameters:
  1158.  *    isf1              I: ISF in the frequency domain (0..6400)
  1159.  *    isf_q             O: quantized ISF
  1160.  *    past_isfq       I/O: past ISF quantizer
  1161.  *    indice            O: quantisation indices (5 words)
  1162.  *    nb_surv           I: number of survivor (1, 2, 3 or 4)
  1163.  *
  1164.  * Function:
  1165.  *    Quantization of isf parameters with prediction. (46 bits)
  1166.  *
  1167.  *    The isf vector is quantized using two-stage VQ with split-by-2 in
  1168.  *    1st stage and split-by-5 in the second stage.
  1169.  * Returns:
  1170.  *    void
  1171.  */
  1172. void E_LPC_isf_2s5s_quantise(Float32 *isf1, Word16 *isf_q, Word16 *past_isfq,
  1173.                              Word32 *indice, Word32 nb_surv)
  1174. {
  1175.    Float32 isf[ORDER], isf_stage2[ORDER];
  1176.    Float32 temp, min_err, distance;
  1177.    Word32 surv1[N_SURV_MAX];     /* indices of survivors from 1st stage */
  1178.    Word32 tmp_ind[5];
  1179.    Word32 i, k;
  1181.    for (i=0; i<ORDER; i++)
  1182.    {
  1183.       isf[i] = (Float32)((isf1[i] - E_ROM_f_mean_isf[i]) -
  1184.          F_MU * past_isfq[i] * 0.390625F);
  1185.    }
  1187.    E_LPC_stage1_isf_vq(&isf[0], E_ROM_dico1_isf, 9, SIZE_BK1, surv1, nb_surv);
  1189.    distance = 1.0e30F;
  1191.    for (k = 0; k < nb_surv; k++)
  1192.    {
  1193.       for (i = 0; i < 9; i++)
  1194.       {
  1195.          isf_stage2[i] = isf[i] - E_ROM_dico1_isf[i + surv1[k] * 9];
  1196.       }
  1198.       tmp_ind[0] = E_LPC_isf_sub_vq(&isf_stage2[0], E_ROM_dico21_isf, 3,
  1199.          SIZE_BK21, &min_err);
  1200.       temp = min_err;
  1201.       tmp_ind[1] = E_LPC_isf_sub_vq(&isf_stage2[3], E_ROM_dico22_isf, 3,
  1202.          SIZE_BK22, &min_err);
  1203.       temp += min_err;
  1204.       tmp_ind[2] = E_LPC_isf_sub_vq(&isf_stage2[6], E_ROM_dico23_isf, 3,
  1205.          SIZE_BK23, &min_err);
  1206.       temp += min_err;
  1208.       if (temp < distance)
  1209.       {
  1210.          distance = temp;
  1211.          indice[0] = surv1[k];
  1213.          for (i = 0; i < 3; i++)
  1214.          {
  1215.             indice[i + 2] = tmp_ind[i];
  1216.          }
  1217.       }
  1218.    }
  1220.    E_LPC_stage1_isf_vq(&isf[9], E_ROM_dico2_isf, 7, SIZE_BK2, surv1, nb_surv);
  1222.    distance = 1.0e30F;
  1224.    for (k=0; k<nb_surv; k++)
  1225.    {
  1226.       for (i = 0; i < 7; i++)
  1227.       {
  1228.          isf_stage2[i] = isf[9+i] - E_ROM_dico2_isf[i+surv1[k]*7];
  1229.       }
  1231.       tmp_ind[0] = E_LPC_isf_sub_vq(&isf_stage2[0], E_ROM_dico24_isf, 3,
  1232.          SIZE_BK24, &min_err);
  1233.       temp = min_err;
  1234.       tmp_ind[1] = E_LPC_isf_sub_vq(&isf_stage2[3], E_ROM_dico25_isf, 4,
  1235.          SIZE_BK25, &min_err);
  1236.       temp += min_err;
  1238.       if (temp < distance)
  1239.       {
  1240.          distance = temp;
  1241.          indice[1] = surv1[k];
  1243.          for (i = 0; i < 2; i++)
  1244.          {
  1245.             indice[i + 5]= tmp_ind[i];
  1246.          }
  1247.       }
  1248.    }
  1250.    /* decoding the ISFs */
  1251.    E_LPC_isf_2s5s_decode(indice, isf_q, past_isfq);
  1253.    return;
  1254. }