pvamrwb_math_op.cpp

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

/* ------------------------------------------------------------------
 * Copyright (C) 2008 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.173
    ANSI-C code for the Adaptive Multi-Rate - Wideband (AMR-WB) speech codec
    Available from http://www.3gpp.org

(C) 2007, 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.
****************************************************************************************/
/*___________________________________________________________________________

    This file contains mathematic operations in fixed point.

    mult_int16_r()     : Same as mult_int16 with rounding
    shr_rnd()          : Same as shr(var1,var2) but with rounding
    div_16by16()       : fractional integer division
    one_ov_sqrt()      : Compute 1/sqrt(L_x)
    one_ov_sqrt_norm() : Compute 1/sqrt(x)
    power_of_2()       : power of 2
    Dot_product12()    : Compute scalar product of <x[],y[]> using accumulator
    Isqrt()            : inverse square root (16 bits precision).
    amrwb_log_2()      : log2 (16 bits precision).

    These operations are not standard double precision operations.
    They are used where low complexity is important and the full 32 bits
    precision is not necessary. For example, the function Div_32() has a
    24 bits precision which is enough for our purposes.

    In this file, the values use theses representations:

    int32 L_32     : standard signed 32 bits format
    int16 hi, lo   : L_32 = hi<<16 + lo<<1  (DPF - Double Precision Format)
    int32 frac, int16 exp : L_32 = frac << exp-31  (normalised format)
    int16 int, frac        : L_32 = int.frac        (fractional format)
 ----------------------------------------------------------------------------*/

#include "pv_amr_wb_type_defs.h"
#include "pvamrwbdecoder_basic_op.h"
#include "pvamrwb_math_op.h"


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

     Function Name : mult_int16_r

     Purpose :

     Same as mult_int16 with rounding, i.e.:
       mult_int16_r(var1,var2) = extract_l(L_shr(((var1 * var2) + 16384),15)) and
       mult_int16_r(-32768,-32768) = 32767.

     Complexity weight : 2

     Inputs :

      var1
               16 bit short signed integer (int16) whose value falls in the
               range : 0xffff 8000 <= var1 <= 0x0000 7fff.

      var2
               16 bit short signed integer (int16) whose value falls in the
               range : 0xffff 8000 <= var1 <= 0x0000 7fff.

     Outputs :

      none

     Return Value :

      var_out
               16 bit short signed integer (int16) whose value falls in the
               range : 0xffff 8000 <= var_out <= 0x0000 7fff.
 ----------------------------------------------------------------------------*/

int16 mult_int16_r(int16 var1, int16 var2)
{
    int32 L_product_arr;

    L_product_arr = (int32) var1 * (int32) var2;      /* product */
    L_product_arr += (int32) 0x00004000L;      /* round */
    L_product_arr >>= 15;       /* shift */
    if ((L_product_arr >> 15) != (L_product_arr >> 31))
    {
        L_product_arr = (L_product_arr >> 31) ^ MAX_16;
    }

    return ((int16)L_product_arr);
}



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

     Function Name : shr_rnd

     Purpose :

     Same as shr(var1,var2) but with rounding. Saturate the result in case of|
     underflows or overflows :
      - If var2 is greater than zero :
            if (sub(shl_int16(shr(var1,var2),1),shr(var1,sub(var2,1))))
            is equal to zero
                       then
                       shr_rnd(var1,var2) = shr(var1,var2)
                       else
                       shr_rnd(var1,var2) = add_int16(shr(var1,var2),1)
      - If var2 is less than or equal to zero :
                       shr_rnd(var1,var2) = shr(var1,var2).

     Complexity weight : 2

     Inputs :

      var1
               16 bit short signed integer (int16) whose value falls in the
               range : 0xffff 8000 <= var1 <= 0x0000 7fff.

      var2
               16 bit short signed integer (int16) whose value falls in the
               range : 0x0000 0000 <= var2 <= 0x0000 7fff.

     Outputs :

      none

     Return Value :

      var_out
               16 bit short signed integer (int16) whose value falls in the
               range : 0xffff 8000 <= var_out <= 0x0000 7fff.
 ----------------------------------------------------------------------------*/

int16 shr_rnd(int16 var1, int16 var2)
{
    int16 var_out;

    var_out = (int16)(var1 >> (var2 & 0xf));
    if (var2)
    {
        if ((var1 & ((int16) 1 << (var2 - 1))) != 0)
        {
            var_out++;
        }
    }
    return (var_out);
}


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

     Function Name : div_16by16

     Purpose :

     Produces a result which is the fractional integer division of var1  by
     var2; var1 and var2 must be positive and var2 must be greater or equal
     to var1; the result is positive (leading bit equal to 0) and truncated
     to 16 bits.
     If var1 = var2 then div(var1,var2) = 32767.

     Complexity weight : 18

     Inputs :

      var1
               16 bit short signed integer (int16) whose value falls in the
               range : 0x0000 0000 <= var1 <= var2 and var2 != 0.

      var2
               16 bit short signed integer (int16) whose value falls in the
               range : var1 <= var2 <= 0x0000 7fff and var2 != 0.

     Outputs :

      none

     Return Value :

      var_out
               16 bit short signed integer (int16) whose value falls in the
               range : 0x0000 0000 <= var_out <= 0x0000 7fff.
               It's a Q15 value (point between b15 and b14).
 ----------------------------------------------------------------------------*/

int16 div_16by16(int16 var1, int16 var2)
{

    int16 var_out = 0;
    register int16 iteration;
    int32 L_num;
    int32 L_denom;
    int32 L_denom_by_2;
    int32 L_denom_by_4;

    if ((var1 > var2) || (var1 < 0))
    {
        return 0; // used to exit(0);
    }
    if (var1)
    {
        if (var1 != var2)
        {

            L_num = (int32) var1;
            L_denom = (int32) var2;
            L_denom_by_2 = (L_denom << 1);
            L_denom_by_4 = (L_denom << 2);
            for (iteration = 5; iteration > 0; iteration--)
            {
                var_out <<= 3;
                L_num   <<= 3;

                if (L_num >= L_denom_by_4)
                {
                    L_num -= L_denom_by_4;
                    var_out |= 4;
                }

                if (L_num >= L_denom_by_2)
                {
                    L_num -= L_denom_by_2;
                    var_out |=  2;
                }

                if (L_num >= (L_denom))
                {
                    L_num -= (L_denom);
                    var_out |=  1;
                }

            }
        }
        else
        {
            var_out = MAX_16;
        }
    }

    return (var_out);

}



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

     Function Name : one_ov_sqrt

         Compute 1/sqrt(L_x).
         if L_x is negative or zero, result is 1 (7fffffff).

  Algorithm:

     1- Normalization of L_x.
     2- call Isqrt_n(L_x, exponant)
     3- L_y = L_x << exponant
 ----------------------------------------------------------------------------*/
int32 one_ov_sqrt(     /* (o) Q31 : output value (range: 0<=val<1)         */
    int32 L_x         /* (i) Q0  : input value  (range: 0<=val<=7fffffff) */
)
{
    int16 exp;
    int32 L_y;

    exp = normalize_amr_wb(L_x);
    L_x <<= exp;                 /* L_x is normalized */
    exp = 31 - exp;

    one_ov_sqrt_norm(&L_x, &exp);

    L_y = shl_int32(L_x, exp);                 /* denormalization   */

    return (L_y);
}

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

     Function Name : one_ov_sqrt_norm

         Compute 1/sqrt(value).
         if value is negative or zero, result is 1 (frac=7fffffff, exp=0).

  Algorithm:

     The function 1/sqrt(value) is approximated by a table and linear
     interpolation.

     1- If exponant is odd then shift fraction right once.
     2- exponant = -((exponant-1)>>1)
     3- i = bit25-b30 of fraction, 16 <= i <= 63 ->because of normalization.
     4- a = bit10-b24
     5- i -=16
     6- fraction = table[i]<<16 - (table[i] - table[i+1]) * a * 2
 ----------------------------------------------------------------------------*/
static const int16 table_isqrt[49] =