diff options
Diffstat (limited to 'opencore/codecs_v2/audio/gsm_amr/amr_wb/dec/src/pvamrwb_math_op.cpp')
-rw-r--r-- | opencore/codecs_v2/audio/gsm_amr/amr_wb/dec/src/pvamrwb_math_op.cpp | 627 |
1 files changed, 627 insertions, 0 deletions
diff --git a/opencore/codecs_v2/audio/gsm_amr/amr_wb/dec/src/pvamrwb_math_op.cpp b/opencore/codecs_v2/audio/gsm_amr/amr_wb/dec/src/pvamrwb_math_op.cpp new file mode 100644 index 0000000..d1ec790 --- /dev/null +++ b/opencore/codecs_v2/audio/gsm_amr/amr_wb/dec/src/pvamrwb_math_op.cpp @@ -0,0 +1,627 @@ +/* ------------------------------------------------------------------ + * 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.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] = +{ + 32767, 31790, 30894, 30070, 29309, 28602, 27945, 27330, 26755, 26214, + 25705, 25225, 24770, 24339, 23930, 23541, 23170, 22817, 22479, 22155, + 21845, 21548, 21263, 20988, 20724, 20470, 20225, 19988, 19760, 19539, + 19326, 19119, 18919, 18725, 18536, 18354, 18176, 18004, 17837, 17674, + 17515, 17361, 17211, 17064, 16921, 16782, 16646, 16514, 16384 +}; + +void one_ov_sqrt_norm( + int32 * frac, /* (i/o) Q31: normalized value (1.0 < frac <= 0.5) */ + int16 * exp /* (i/o) : exponent (value = frac x 2^exponent) */ +) +{ + int16 i, a, tmp; + + + if (*frac <= (int32) 0) + { + *exp = 0; + *frac = 0x7fffffffL; + return; + } + + if ((*exp & 1) == 1) /* If exponant odd -> shift right */ + *frac >>= 1; + + *exp = negate_int16((*exp - 1) >> 1); + + *frac >>= 9; + i = extract_h(*frac); /* Extract b25-b31 */ + *frac >>= 1; + a = (int16)(*frac); /* Extract b10-b24 */ + a = (int16)(a & (int16) 0x7fff); + + i -= 16; + + *frac = L_deposit_h(table_isqrt[i]); /* table[i] << 16 */ + tmp = table_isqrt[i] - table_isqrt[i + 1]; /* table[i] - table[i+1]) */ + + *frac = msu_16by16_from_int32(*frac, tmp, a); /* frac -= tmp*a*2 */ + + return; +} + +/*---------------------------------------------------------------------------- + + Function Name : power_2() + + L_x = pow(2.0, exponant.fraction) (exponant = interger part) + = pow(2.0, 0.fraction) << exponant + + Algorithm: + + The function power_2(L_x) is approximated by a table and linear + interpolation. + + 1- i = bit10-b15 of fraction, 0 <= i <= 31 + 2- a = bit0-b9 of fraction + 3- L_x = table[i]<<16 - (table[i] - table[i+1]) * a * 2 + 4- L_x = L_x >> (30-exponant) (with rounding) + ----------------------------------------------------------------------------*/ +const int16 table_pow2[33] = +{ + 16384, 16743, 17109, 17484, 17867, 18258, 18658, 19066, 19484, 19911, + 20347, 20792, 21247, 21713, 22188, 22674, 23170, 23678, 24196, 24726, + 25268, 25821, 26386, 26964, 27554, 28158, 28774, 29405, 30048, 30706, + 31379, 32066, 32767 +}; + +int32 power_of_2( /* (o) Q0 : result (range: 0<=val<=0x7fffffff) */ + int16 exponant, /* (i) Q0 : Integer part. (range: 0<=val<=30) */ + int16 fraction /* (i) Q15 : Fractionnal part. (range: 0.0<=val<1.0) */ +) +{ + int16 exp, i, a, tmp; + int32 L_x; + + L_x = fraction << 5; /* L_x = fraction<<6 */ + i = (fraction >> 10); /* Extract b10-b16 of fraction */ + a = (int16)(L_x); /* Extract b0-b9 of fraction */ + a = (int16)(a & (int16) 0x7fff); + + L_x = ((int32)table_pow2[i]) << 15; /* table[i] << 16 */ + tmp = table_pow2[i] - table_pow2[i + 1]; /* table[i] - table[i+1] */ + L_x -= ((int32)tmp * a); /* L_x -= tmp*a*2 */ + + exp = 29 - exponant ; + + if (exp) + { + L_x = ((L_x >> exp) + ((L_x >> (exp - 1)) & 1)); + } + + return (L_x); +} + +/*---------------------------------------------------------------------------- + * + * Function Name : Dot_product12() + * + * Compute scalar product of <x[],y[]> using accumulator. + * + * The result is normalized (in Q31) with exponent (0..30). + * + * Algorithm: + * + * dot_product = sum(x[i]*y[i]) i=0..N-1 + ----------------------------------------------------------------------------*/ + +int32 Dot_product12( /* (o) Q31: normalized result (1 < val <= -1) */ + int16 x[], /* (i) 12bits: x vector */ + int16 y[], /* (i) 12bits: y vector */ + int16 lg, /* (i) : vector length */ + int16 * exp /* (o) : exponent of result (0..+30) */ +) +{ + int16 i, sft; + int32 L_sum; + int16 *pt_x = x; + int16 *pt_y = y; + + L_sum = 1L; + + + for (i = lg >> 3; i != 0; i--) + { + L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++)); + L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++)); + L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++)); + L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++)); + L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++)); + L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++)); + L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++)); + L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++)); + } + + /* Normalize acc in Q31 */ + + sft = normalize_amr_wb(L_sum); + L_sum <<= sft; + + *exp = 30 - sft; /* exponent = 0..30 */ + + return (L_sum); +} + +/* Table for Log2() */ +const int16 Log2_norm_table[33] = +{ + 0, 1455, 2866, 4236, 5568, 6863, 8124, 9352, 10549, 11716, + 12855, 13967, 15054, 16117, 17156, 18172, 19167, 20142, 21097, 22033, + 22951, 23852, 24735, 25603, 26455, 27291, 28113, 28922, 29716, 30497, + 31266, 32023, 32767 +}; + +/*---------------------------------------------------------------------------- + * + * FUNCTION: Lg2_normalized() + * + * PURPOSE: Computes log2(L_x, exp), where L_x is positive and + * normalized, and exp is the normalisation exponent + * If L_x is negative or zero, the result is 0. + * + * DESCRIPTION: + * The function Log2(L_x) is approximated by a table and linear + * interpolation. The following steps are used to compute Log2(L_x) + * + * 1- exponent = 30-norm_exponent + * 2- i = bit25-b31 of L_x; 32<=i<=63 (because of normalization). + * 3- a = bit10-b24 + * 4- i -=32 + * 5- fraction = table[i]<<16 - (table[i] - table[i+1]) * a * 2 + * +----------------------------------------------------------------------------*/ +void Lg2_normalized( + int32 L_x, /* (i) : input value (normalized) */ + int16 exp, /* (i) : norm_l (L_x) */ + int16 *exponent, /* (o) : Integer part of Log2. (range: 0<=val<=30) */ + int16 *fraction /* (o) : Fractional part of Log2. (range: 0<=val<1) */ +) +{ + int16 i, a, tmp; + int32 L_y; + + if (L_x <= (int32) 0) + { + *exponent = 0; + *fraction = 0;; + return; + } + + *exponent = 30 - exp; + + L_x >>= 9; + i = extract_h(L_x); /* Extract b25-b31 */ + L_x >>= 1; + a = (int16)(L_x); /* Extract b10-b24 of fraction */ + a &= 0x7fff; + + i -= 32; + + L_y = L_deposit_h(Log2_norm_table[i]); /* table[i] << 16 */ + tmp = Log2_norm_table[i] - Log2_norm_table[i + 1]; /* table[i] - table[i+1] */ + L_y = msu_16by16_from_int32(L_y, tmp, a); /* L_y -= tmp*a*2 */ + + *fraction = extract_h(L_y); + + return; +} + + + +/*---------------------------------------------------------------------------- + * + * FUNCTION: amrwb_log_2() + * + * PURPOSE: Computes log2(L_x), where L_x is positive. + * If L_x is negative or zero, the result is 0. + * + * DESCRIPTION: + * normalizes L_x and then calls Lg2_normalized(). + * + ----------------------------------------------------------------------------*/ +void amrwb_log_2( + int32 L_x, /* (i) : input value */ + int16 *exponent, /* (o) : Integer part of Log2. (range: 0<=val<=30) */ + int16 *fraction /* (o) : Fractional part of Log2. (range: 0<=val<1) */ +) +{ + int16 exp; + + exp = normalize_amr_wb(L_x); + Lg2_normalized(shl_int32(L_x, exp), exp, exponent, fraction); +} + + +/***************************************************************************** + * + * These operations are not standard double precision operations. * + * They are used where single precision is not enough but 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. * + * * + * The double precision numbers use a special representation: * + * * + * L_32 = hi<<16 + lo<<1 * + * * + * L_32 is a 32 bit integer. * + * hi and lo are 16 bit signed integers. * + * As the low part also contains the sign, this allows fast multiplication. * + * * + * 0x8000 0000 <= L_32 <= 0x7fff fffe. * + * * + * We will use DPF (Double Precision Format )in this file to specify * + * this special format. * + ***************************************************************************** +*/ + + +/*---------------------------------------------------------------------------- + * + * Function int32_to_dpf() + * + * Extract from a 32 bit integer two 16 bit DPF. + * + * Arguments: + * + * L_32 : 32 bit integer. + * 0x8000 0000 <= L_32 <= 0x7fff ffff. + * hi : b16 to b31 of L_32 + * lo : (L_32 - hi<<16)>>1 + * + ----------------------------------------------------------------------------*/ + +void int32_to_dpf(int32 L_32, int16 *hi, int16 *lo) +{ + *hi = (int16)(L_32 >> 16); + *lo = (int16)((L_32 - (*hi << 16)) >> 1); + return; +} + + +/*---------------------------------------------------------------------------- + * Function mpy_dpf_32() + * + * Multiply two 32 bit integers (DPF). The result is divided by 2**31 + * + * L_32 = (hi1*hi2)<<1 + ( (hi1*lo2)>>15 + (lo1*hi2)>>15 )<<1 + * + * This operation can also be viewed as the multiplication of two Q31 + * number and the result is also in Q31. + * + * Arguments: + * + * hi1 hi part of first number + * lo1 lo part of first number + * hi2 hi part of second number + * lo2 lo part of second number + * + ----------------------------------------------------------------------------*/ + +int32 mpy_dpf_32(int16 hi1, int16 lo1, int16 hi2, int16 lo2) +{ + int32 L_32; + + L_32 = mul_16by16_to_int32(hi1, hi2); + L_32 = mac_16by16_to_int32(L_32, mult_int16(hi1, lo2), 1); + L_32 = mac_16by16_to_int32(L_32, mult_int16(lo1, hi2), 1); + + return (L_32); +} + + |