az_lsp.cpp

实现3GPP的GSM中AMR语音的CODECS。

 None.

------------------------------------------------------------------------------
 REFERENCES

 az_lsp.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001

------------------------------------------------------------------------------
 PSEUDO-CODE

void Az_lsp (
    Word16 a[],         // (i)  : predictor coefficients (MP1)
    Word16 lsp[],       // (o)  : line spectral pairs (M)
    Word16 old_lsp[]    // (i)  : old lsp[] (in case not found 10 roots) (M)
)
{
    Word16 i, j, nf, ip;
    Word16 xlow, ylow, xhigh, yhigh, xmid, ymid, xint;
    Word16 x, y, sign, exp;
    Word16 *coef;
    Word16 f1[M / 2 + 1], f2[M / 2 + 1];
    Word32 t0;

     *-------------------------------------------------------------*
     *  find the sum and diff. pol. F1(z) and F2(z)                *
     *    F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1)  *
     *                                                             *
     * f1[0] = 1.0;                                                *
     * f2[0] = 1.0;                                                *
     *                                                             *
     * for (i = 0; i< NC; i++)                                     *
     * {                                                           *
     *   f1[i+1] = a[i+1] + a[M-i] - f1[i] ;                       *
     *   f2[i+1] = a[i+1] - a[M-i] + f2[i] ;                       *
     * }                                                           *
     *-------------------------------------------------------------*

    f1[0] = 1024; // f1[0] = 1.0
    f2[0] = 1024; // f2[0] = 1.0

// The reference ETSI code uses a global flag for Overflow. However, in the
// actual implementation a pointer to Overflow flag is passed in as a
// parameter to the function. This pointer is passed into all the basic math
// functions invoked

    for (i = 0; i < NC; i++)
    {
        t0 = L_mult (a[i + 1], 8192);   // x = (a[i+1] + a[M-i]) >> 2
        t0 = L_mac (t0, a[M - i], 8192);
        x = extract_h (t0);
        // f1[i+1] = a[i+1] + a[M-i] - f1[i]
        f1[i + 1] = sub (x, f1[i]);

        t0 = L_mult (a[i + 1], 8192);   // x = (a[i+1] - a[M-i]) >> 2
        t0 = L_msu (t0, a[M - i], 8192);
        x = extract_h (t0);
        // f2[i+1] = a[i+1] - a[M-i] + f2[i]
        f2[i + 1] = add (x, f2[i]);
    }

     *-------------------------------------------------------------*
     * find the LSPs using the Chebychev pol. evaluation           *
     *-------------------------------------------------------------*

    nf = 0; // number of found frequencies
    ip = 0; // indicator for f1 or f2

    coef = f1;

    xlow = grid[0];
    ylow = Chebps (xlow, coef, NC);

    j = 0;
    // while ( (nf < M) && (j < grid_points) )
    while ((sub (nf, M) < 0) && (sub (j, grid_points) < 0))
    {
        j++;
        xhigh = xlow;
        yhigh = ylow;
        xlow = grid[j];
        ylow = Chebps (xlow, coef, NC);

        if (L_mult (ylow, yhigh) <= (Word32) 0L)
        {

            // divide 4 times the interval

            for (i = 0; i < 4; i++)
            {
                // xmid = (xlow + xhigh)/2
                xmid = add (shr (xlow, 1), shr (xhigh, 1));
                ymid = Chebps (xmid, coef, NC);

                if (L_mult (ylow, ymid) <= (Word32) 0L)
                {
                    yhigh = ymid;
                    xhigh = xmid;
                }
                else
                {
                    ylow = ymid;
                    xlow = xmid;
                }
            }

             *-------------------------------------------------------------*
             * Linear interpolation                                        *
             *    xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow);            *
             *-------------------------------------------------------------*

            x = sub (xhigh, xlow);
            y = sub (yhigh, ylow);

            if (y == 0)
            {
                xint = xlow;
            }
            else
            {
                sign = y;
                y = abs_s (y);
                exp = norm_s (y);
                y = shl (y, exp);
                y = div_s ((Word16) 16383, y);
                t0 = L_mult (x, y);
                t0 = L_shr (t0, sub (20, exp));
                y = extract_l (t0);     // y= (xhigh-xlow)/(yhigh-ylow)

                if (sign < 0)
                    y = negate (y);

                t0 = L_mult (ylow, y);
                t0 = L_shr (t0, 11);
                xint = sub (xlow, extract_l (t0)); // xint = xlow - ylow*y
            }

            lsp[nf] = xint;
            xlow = xint;
            nf++;

            if (ip == 0)
            {
                ip = 1;
                coef = f2;
            }
            else
            {
                ip = 0;
                coef = f1;
            }
            ylow = Chebps (xlow, coef, NC);

        }
    }

    // Check if M roots found

    if (sub (nf, M) < 0)
    {
        for (i = 0; i < M; i++)
        {
            lsp[i] = old_lsp[i];
        }

    }
    return;
}

------------------------------------------------------------------------------
 RESOURCES USED [optional]

 When the code is written for a specific target processor the
 the resources used should be documented below.

 HEAP MEMORY USED: x bytes

 STACK MEMORY USED: x bytes

 CLOCK CYCLES: (cycle count equation for this function) + (variable
                used to represent cycle count for each subroutine
                called)
     where: (cycle count variable) = cycle count for [subroutine
                                     name]

------------------------------------------------------------------------------
 CAUTION [optional]
 [State any special notes, constraints or cautions for users of this function]

------------------------------------------------------------------------------
*/

void Az_lsp(
    Word16 a[],         /* (i)  : predictor coefficients (MP1)               */
    Word16 lsp[],       /* (o)  : line spectral pairs (M)                    */
    Word16 old_lsp[],   /* (i)  : old lsp[] (in case not found 10 roots) (M) */
    Flag   *pOverflow   /* (i/o): overflow flag                              */
)
{
    register Word16 i;
    register Word16 j;
    register Word16 nf;
    register Word16 ip;
    Word16 xlow;
    Word16 ylow;
    Word16 xhigh;
    Word16 yhigh;
    Word16 xmid;
    Word16 ymid;
    Word16 xint;
    Word16 x;
    Word16 y;
    Word16 sign;
    Word16 exp;
    Word16 *coef;
    Word16 f1[NC + 1];
    Word16 f2[NC + 1];
    Word32 L_temp1;
    Word32 L_temp2;
    Word16 *p_f1 = f1;
    Word16 *p_f2 = f2;

    /*-------------------------------------------------------------*
     *  find the sum and diff. pol. F1(z) and F2(z)                *
     *    F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1)  *
     *                                                             *
     * f1[0] = 1.0;                                                *
     * f2[0] = 1.0;                                                *
     *                                                             *
     * for (i = 0; i< NC; i++)                                     *
     * {                                                           *
     *   f1[i+1] = a[i+1] + a[M-i] - f1[i] ;                       *
     *   f2[i+1] = a[i+1] - a[M-i] + f2[i] ;                       *
     * }                                                           *
     *-------------------------------------------------------------*/

    *p_f1 = 1024;                       /* f1[0] = 1.0 */
    *p_f2 = 1024;                       /* f2[0] = 1.0 */

    for (i = 0; i < NC; i++)
    {
        L_temp1 = (Word32) * (a + i + 1);
        L_temp2 = (Word32) * (a + M - i);
        /* x = (a[i+1] + a[M-i]) >> 2  */
        x = (Word16)((L_temp1 + L_temp2) >> 2);
        /* y = (a[i+1] - a[M-i]) >> 2 */
        y = (Word16)((L_temp1 - L_temp2) >> 2);
        /* f1[i+1] = a[i+1] + a[M-i] - f1[i] */
        x -= *(p_f1++);
        *(p_f1) = x;
        /* f2[i+1] = a[i+1] - a[M-i] + f2[i] */
        y += *(p_f2++);
        *(p_f2) = y;
    }

    /*-------------------------------------------------------------*
     * find the LSPs using the Chebychev pol. evaluation           *
     *-------------------------------------------------------------*/

    nf = 0;                         /* number of found frequencies */
    ip = 0;                         /* indicator for f1 or f2      */

    coef = f1;

    xlow = *(grid);
    ylow = Chebps(xlow, coef, NC, pOverflow);

    j = 0;

    while ((nf < M) && (j < grid_points))
    {
        j++;
        xhigh = xlow;
        yhigh = ylow;
        xlow = *(grid + j);
        ylow = Chebps(xlow, coef, NC, pOverflow);

        if (((Word32)ylow*yhigh) <= 0)
        {
            /* divide 4 times the interval */
            for (i = 4; i != 0; i--)
            {
                /* xmid = (xlow + xhigh)/2 */
                x = xlow >> 1;
                y = xhigh >> 1;
                xmid = x + y;

                ymid = Chebps(xmid, coef, NC, pOverflow);

                if (((Word32)ylow*ymid) <= 0)
                {
                    yhigh = ymid;
                    xhigh = xmid;
                }
                else
                {
                    ylow = ymid;
                    xlow = xmid;
                }
            }

            /*-------------------------------------------------------------*
             * Linear interpolation                                        *
             *    xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow);            *
             *-------------------------------------------------------------*/

            x = xhigh - xlow;
            y = yhigh - ylow;

            if (y == 0)
            {
                xint = xlow;
            }
            else
            {
                sign = y;
                y = abs_s(y);
                exp = norm_s(y);
                y <<= exp;
                y = div_s((Word16) 16383, y);

                y = ((Word32)x * y) >> (19 - exp);

                if (sign < 0)
                {
                    y = -y;
                }

                /* xint = xlow - ylow*y */
                xint = xlow - (((Word32) ylow * y) >> 10);
            }

            *(lsp + nf) = xint;
            xlow = xint;
            nf++;

            if (ip == 0)
            {
                ip = 1;
                coef = f2;
            }
            else
            {
                ip = 0;
                coef = f1;
            }

            ylow = Chebps(xlow, coef, NC, pOverflow);

        }
    }

    /* Check if M roots found */

    if (nf < M)
    {
        for (i = NC; i != 0 ; i--)
        {
            *lsp++ = *old_lsp++;
            *lsp++ = *old_lsp++;
        }
    }

}

Word16 Chebps_Wrapper(Word16 x,
                      Word16 f[], /* (n) */
                      Word16 n,
                      Flag *pOverflow)
{
    return Chebps(x, f, n, pOverflow);
}