/* ------------------------------------------------------------------
 * Copyright (C) 1998-2009 PacketVideo
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
 * express or implied.
 * See the License for the specific language governing permissions
 * and limitations under the License.
 * -------------------------------------------------------------------
 */
/****************************************************************************************
Portions of this file are derived from the following 3GPP standard:

    3GPP TS 26.073
    ANSI-C code for the Adaptive Multi-Rate (AMR) speech codec
    Available from http://www.3gpp.org

(C) 2004, 3GPP Organizational Partners (ARIB, ATIS, CCSA, ETSI, TTA, TTC)
Permission to distribute, modify and use this file under the standard license
terms listed above has been obtained from the copyright holder.
****************************************************************************************/
/*

 Pathname: ./audio/gsm-amr/c/src/az_lsp.c
 Funtions: Chebps
           Chebps_Wrapper
           Az_lsp

------------------------------------------------------------------------------
 REVISION HISTORY

 Description: Finished first pass of optimization.

 Description: Made changes based on review comments.

 Description: Made input to Chebps_Wrapper consistent with that of Chebps.

 Description: Replaced current Pseudo-code with the UMTS code version 3.2.0.
              Updated coding template.

 Description: Replaced basic_op.h and oper_32b.h with the header files of the
              math functions used by the file.

 Description: Made the following changes per comments from Phase 2/3 review:
              1. Used "-" operator instead of calling sub function in the
                 az_lsp() code.
              2. Copied detailed function description of az_lsp from the
                 header file.
              3. Modified local variable definition to one per line.
              4. Used NC in the definition of f1 and f2 arrays.
              5. Added curly brackets in the IF statement.

 Description: Changed function interface to pass in a pointer to overflow
              flag into the function instead of using a global flag. Removed
              inclusion of unneeded header files.

 Description:  For Chebps() and Az_lsp()
              1. Eliminated unused include files.
              2. Replaced array addressing by pointers
              3. Eliminated math operations that unnecessary checked for
                 saturation.
              4. Eliminated not needed variables
              5. Eliminated if-else checks for saturation
              6. Deleted unused function cheps_wraper

 Description:  Added casting to eliminate warnings


 Who:                           Date:
 Description:

------------------------------------------------------------------------------
 MODULE DESCRIPTION

 These modules compute the LSPs from the LP coefficients.

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


/*----------------------------------------------------------------------------
; INCLUDES
----------------------------------------------------------------------------*/
#include "az_lsp.h"
#include "cnst.h"
#include "basic_op.h"

/*----------------------------------------------------------------------------
; MACROS
; Define module specific macros here
----------------------------------------------------------------------------*/


/*----------------------------------------------------------------------------
; DEFINES
; Include all pre-processor statements here. Include conditional
; compile variables also.
----------------------------------------------------------------------------*/
#define NC   (M/2)                  /* M = LPC order, NC = M/2 */

/*----------------------------------------------------------------------------
; LOCAL FUNCTION DEFINITIONS
; Function Prototype declaration
----------------------------------------------------------------------------*/


/*----------------------------------------------------------------------------
; LOCAL VARIABLE DEFINITIONS
; Variable declaration - defined here and used outside this module
----------------------------------------------------------------------------*/

/*
------------------------------------------------------------------------------
 FUNCTION NAME: Chebps
------------------------------------------------------------------------------
 INPUT AND OUTPUT DEFINITIONS

 Inputs:
    x = input value (Word16)
    f = polynomial (Word16)
    n = polynomial order (Word16)

    pOverflow = pointer to overflow (Flag)

 Outputs:
    pOverflow -> 1 if the operations in the function resulted in saturation.

 Returns:
    cheb = Chebyshev polynomial for the input value x.(Word16)

 Global Variables Used:
    None.

 Local Variables Needed:
    None.

------------------------------------------------------------------------------
 FUNCTION DESCRIPTION

 This module evaluates the Chebyshev polynomial series.
    - The polynomial order is   n = m/2 = 5
    - The polynomial F(z) (F1(z) or F2(z)) is given by
        F(w) = 2 exp(-j5w) C(x)
        where
        C(x) = T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2
        and T_m(x) = cos(mw) is the mth order Chebyshev
        polynomial ( x=cos(w) )
    - C(x) for the input x is returned.

------------------------------------------------------------------------------
 REQUIREMENTS

 None.

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

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

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

static Word16 Chebps (Word16 x,
                      Word16 f[], // (n)
                      Word16 n)
{
    Word16 i, cheb;
    Word16 b0_h, b0_l, b1_h, b1_l, b2_h, b2_l;
    Word32 t0;

// 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

    b2_h = 256; // b2 = 1.0
    b2_l = 0;

    t0 = L_mult (x, 512);          // 2*x
    t0 = L_mac (t0, f[1], 8192);   // + f[1]
    L_Extract (t0, &b1_h, &b1_l);  // b1 = 2*x + f[1]

    for (i = 2; i < n; i++)
    {
        t0 = Mpy_32_16 (b1_h, b1_l, x);         // t0 = 2.0*x*b1
        t0 = L_shl (t0, 1);
        t0 = L_mac (t0, b2_h, (Word16) 0x8000); // t0 = 2.0*x*b1 - b2
        t0 = L_msu (t0, b2_l, 1);
        t0 = L_mac (t0, f[i], 8192);            // t0 = 2.0*x*b1 - b2 + f[i]

        L_Extract (t0, &b0_h, &b0_l);           // b0 = 2.0*x*b1 - b2 + f[i]

        b2_l = b1_l; // b2 = b1;
        b2_h = b1_h;
        b1_l = b0_l; // b1 = b0;
        b1_h = b0_h;
    }

    t0 = Mpy_32_16 (b1_h, b1_l, x);             // t0 = x*b1;
    t0 = L_mac (t0, b2_h, (Word16) 0x8000);     // t0 = x*b1 - b2
    t0 = L_msu (t0, b2_l, 1);
    t0 = L_mac (t0, f[i], 4096);                // t0 = x*b1 - b2 + f[i]/2

    t0 = L_shl (t0, 6);

    cheb = extract_h (t0);

    return (cheb);
}

------------------------------------------------------------------------------
 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]

------------------------------------------------------------------------------
*/
static Word16 Chebps(Word16 x,
                     Word16 f[], /* (n) */
                     Word16 n,
                     Flag *pOverflow)
{
    Word16 i;
    Word16 cheb;
    Word16 b1_h;
    Word16 b1_l;
    Word32 t0;
    Word32 L_temp;
    Word16 *p_f = &f[1];

    OSCL_UNUSED_ARG(pOverflow);

    /* L_temp = 1.0 */

    L_temp = 0x01000000L;

    t0 = ((Word32) x << 10) + ((Word32) * (p_f++) << 14);

    /* b1 = t0 = 2*x + f[1]  */

    b1_h = (Word16)(t0 >> 16);
    b1_l = (Word16)((t0 >> 1) - (b1_h << 15));


    for (i = 2; i < n; i++)
    {
        /* t0 = 2.0*x*b1    */
        t0  = ((Word32) b1_h * x);
        t0 += ((Word32) b1_l * x) >> 15;
        t0 <<= 2;

        /* t0 = 2.0*x*b1 - b2   */
        t0 -= L_temp;

        /* t0 = 2.0*x*b1 - b2 + f[i] */
        t0 += (Word32) * (p_f++) << 14;

        L_temp = ((Word32) b1_h << 16) + ((Word32) b1_l << 1);

        /* b0 = 2.0*x*b1 - b2 + f[i]*/
        b1_h = (Word16)(t0 >> 16);
        b1_l = (Word16)((t0 >> 1) - (b1_h << 15));

    }

    /* t0 = x*b1; */
    t0  = ((Word32) b1_h * x);
    t0 += ((Word32) b1_l * x) >> 15;
    t0 <<= 1;


    /* t0 = x*b1 - b2   */
    t0 -= L_temp;

    /* t0 = x*b1 - b2 + f[i]/2 */
    t0 += (Word32) * (p_f) << 13;


    if ((UWord32)(t0 - 0xfe000000L) < (UWord32)0x03ffffffL)
    {
        cheb = (Word16)(t0 >> 10);
    }
    else
    {
        if (t0 > (Word32) 0x01ffffffL)
        {
            cheb = MAX_16;

        }
        else
        {
            cheb = MIN_16;
        }
    }

    return (cheb);
}


/*
------------------------------------------------------------------------------
 FUNCTION NAME: Az_lsp
------------------------------------------------------------------------------
 INPUT AND OUTPUT DEFINITIONS FOR Az_lsp

 Inputs:
    a = predictor coefficients (Word16)
    lsp = line spectral pairs (Word16)
    old_lsp = old line spectral pairs (Word16)

    pOverflow = pointer to overflow (Flag)

 Outputs:
    pOverflow -> 1 if the operations in the function resulted in saturation.

 Returns:
    None.

 Global Variables Used:
    None.

 Local Variables Needed:
    None.

------------------------------------------------------------------------------
 FUNCTION DESCRIPTION

 This function computes the LSPs from the LP coefficients.

 The sum and difference filters are computed and divided by 1+z^{-1} and
 1-z^{-1}, respectively.

     f1[i] = a[i] + a[11-i] - f1[i-1] ;   i=1,...,5
     f2[i] = a[i] - a[11-i] + f2[i-1] ;   i=1,...,5

 The roots of F1(z) and F2(z) are found using Chebyshev polynomial evaluation.
 The polynomials are evaluated at 60 points regularly spaced in the
 frequency domain. The sign change interval is subdivided 4 times to better
 track the root. The LSPs are found in the cosine domain [1,-1].

 If less than 10 roots are found, the LSPs from the past frame are used.

------------------------------------------------------------------------------
 REQUIREMENTS

 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                              */
)
{
    Word16 i;
    Word16 j;
    Word16 nf;
    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);
}