mirror of https://github.com/mstorsjo/fdk-aac.git
922 lines
28 KiB
C
922 lines
28 KiB
C
/* -----------------------------------------------------------------------------
|
|
Software License for The Fraunhofer FDK AAC Codec Library for Android
|
|
|
|
© Copyright 1995 - 2018 Fraunhofer-Gesellschaft zur Förderung der angewandten
|
|
Forschung e.V. All rights reserved.
|
|
|
|
1. INTRODUCTION
|
|
The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software
|
|
that implements the MPEG Advanced Audio Coding ("AAC") encoding and decoding
|
|
scheme for digital audio. This FDK AAC Codec software is intended to be used on
|
|
a wide variety of Android devices.
|
|
|
|
AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient
|
|
general perceptual audio codecs. AAC-ELD is considered the best-performing
|
|
full-bandwidth communications codec by independent studies and is widely
|
|
deployed. AAC has been standardized by ISO and IEC as part of the MPEG
|
|
specifications.
|
|
|
|
Patent licenses for necessary patent claims for the FDK AAC Codec (including
|
|
those of Fraunhofer) may be obtained through Via Licensing
|
|
(www.vialicensing.com) or through the respective patent owners individually for
|
|
the purpose of encoding or decoding bit streams in products that are compliant
|
|
with the ISO/IEC MPEG audio standards. Please note that most manufacturers of
|
|
Android devices already license these patent claims through Via Licensing or
|
|
directly from the patent owners, and therefore FDK AAC Codec software may
|
|
already be covered under those patent licenses when it is used for those
|
|
licensed purposes only.
|
|
|
|
Commercially-licensed AAC software libraries, including floating-point versions
|
|
with enhanced sound quality, are also available from Fraunhofer. Users are
|
|
encouraged to check the Fraunhofer website for additional applications
|
|
information and documentation.
|
|
|
|
2. COPYRIGHT LICENSE
|
|
|
|
Redistribution and use in source and binary forms, with or without modification,
|
|
are permitted without payment of copyright license fees provided that you
|
|
satisfy the following conditions:
|
|
|
|
You must retain the complete text of this software license in redistributions of
|
|
the FDK AAC Codec or your modifications thereto in source code form.
|
|
|
|
You must retain the complete text of this software license in the documentation
|
|
and/or other materials provided with redistributions of the FDK AAC Codec or
|
|
your modifications thereto in binary form. You must make available free of
|
|
charge copies of the complete source code of the FDK AAC Codec and your
|
|
modifications thereto to recipients of copies in binary form.
|
|
|
|
The name of Fraunhofer may not be used to endorse or promote products derived
|
|
from this library without prior written permission.
|
|
|
|
You may not charge copyright license fees for anyone to use, copy or distribute
|
|
the FDK AAC Codec software or your modifications thereto.
|
|
|
|
Your modified versions of the FDK AAC Codec must carry prominent notices stating
|
|
that you changed the software and the date of any change. For modified versions
|
|
of the FDK AAC Codec, the term "Fraunhofer FDK AAC Codec Library for Android"
|
|
must be replaced by the term "Third-Party Modified Version of the Fraunhofer FDK
|
|
AAC Codec Library for Android."
|
|
|
|
3. NO PATENT LICENSE
|
|
|
|
NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without
|
|
limitation the patents of Fraunhofer, ARE GRANTED BY THIS SOFTWARE LICENSE.
|
|
Fraunhofer provides no warranty of patent non-infringement with respect to this
|
|
software.
|
|
|
|
You may use this FDK AAC Codec software or modifications thereto only for
|
|
purposes that are authorized by appropriate patent licenses.
|
|
|
|
4. DISCLAIMER
|
|
|
|
This FDK AAC Codec software is provided by Fraunhofer on behalf of the copyright
|
|
holders and contributors "AS IS" and WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES,
|
|
including but not limited to the implied warranties of merchantability and
|
|
fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
|
|
CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary,
|
|
or consequential damages, including but not limited to procurement of substitute
|
|
goods or services; loss of use, data, or profits, or business interruption,
|
|
however caused and on any theory of liability, whether in contract, strict
|
|
liability, or tort (including negligence), arising in any way out of the use of
|
|
this software, even if advised of the possibility of such damage.
|
|
|
|
5. CONTACT INFORMATION
|
|
|
|
Fraunhofer Institute for Integrated Circuits IIS
|
|
Attention: Audio and Multimedia Departments - FDK AAC LL
|
|
Am Wolfsmantel 33
|
|
91058 Erlangen, Germany
|
|
|
|
www.iis.fraunhofer.de/amm
|
|
amm-info@iis.fraunhofer.de
|
|
----------------------------------------------------------------------------- */
|
|
|
|
/******************* Library for basic calculation routines ********************
|
|
|
|
Author(s): M. Gayer
|
|
|
|
Description: Fixed point specific mathematical functions
|
|
|
|
*******************************************************************************/
|
|
|
|
#ifndef FIXPOINT_MATH_H
|
|
#define FIXPOINT_MATH_H
|
|
|
|
#include "common_fix.h"
|
|
#include "scale.h"
|
|
|
|
/*
|
|
* Data definitions
|
|
*/
|
|
|
|
#define LD_DATA_SCALING (64.0f)
|
|
#define LD_DATA_SHIFT 6 /* pow(2, LD_DATA_SHIFT) = LD_DATA_SCALING */
|
|
|
|
#define MAX_LD_PRECISION 10
|
|
#define LD_PRECISION 10
|
|
|
|
/* Taylor series coefficients for ln(1-x), centered at 0 (MacLaurin polynomial).
|
|
*/
|
|
#ifndef LDCOEFF_16BIT
|
|
LNK_SECTION_CONSTDATA_L1
|
|
static const FIXP_DBL ldCoeff[MAX_LD_PRECISION] = {
|
|
FL2FXCONST_DBL(-1.0), FL2FXCONST_DBL(-1.0 / 2.0),
|
|
FL2FXCONST_DBL(-1.0 / 3.0), FL2FXCONST_DBL(-1.0 / 4.0),
|
|
FL2FXCONST_DBL(-1.0 / 5.0), FL2FXCONST_DBL(-1.0 / 6.0),
|
|
FL2FXCONST_DBL(-1.0 / 7.0), FL2FXCONST_DBL(-1.0 / 8.0),
|
|
FL2FXCONST_DBL(-1.0 / 9.0), FL2FXCONST_DBL(-1.0 / 10.0)};
|
|
#else /* LDCOEFF_16BIT */
|
|
LNK_SECTION_CONSTDATA_L1
|
|
static const FIXP_SGL ldCoeff[MAX_LD_PRECISION] = {
|
|
FL2FXCONST_SGL(-1.0), FL2FXCONST_SGL(-1.0 / 2.0),
|
|
FL2FXCONST_SGL(-1.0 / 3.0), FL2FXCONST_SGL(-1.0 / 4.0),
|
|
FL2FXCONST_SGL(-1.0 / 5.0), FL2FXCONST_SGL(-1.0 / 6.0),
|
|
FL2FXCONST_SGL(-1.0 / 7.0), FL2FXCONST_SGL(-1.0 / 8.0),
|
|
FL2FXCONST_SGL(-1.0 / 9.0), FL2FXCONST_SGL(-1.0 / 10.0)};
|
|
#endif /* LDCOEFF_16BIT */
|
|
|
|
/*****************************************************************************
|
|
|
|
functionname: invSqrtNorm2
|
|
description: delivers 1/sqrt(op) normalized to .5...1 and the shift value
|
|
of the OUTPUT
|
|
|
|
*****************************************************************************/
|
|
#define SQRT_BITS 7
|
|
#define SQRT_VALUES (128 + 2)
|
|
#define SQRT_BITS_MASK 0x7f
|
|
#define SQRT_FRACT_BITS_MASK 0x007FFFFF
|
|
|
|
extern const FIXP_DBL invSqrtTab[SQRT_VALUES];
|
|
|
|
/*
|
|
* Hardware specific implementations
|
|
*/
|
|
|
|
#if defined(__x86__)
|
|
#include "x86/fixpoint_math_x86.h"
|
|
#endif /* target architecture selector */
|
|
|
|
/*
|
|
* Fallback implementations
|
|
*/
|
|
#if !defined(FUNCTION_fIsLessThan)
|
|
/**
|
|
* \brief Compares two fixpoint values incl. scaling.
|
|
* \param a_m mantissa of the first input value.
|
|
* \param a_e exponent of the first input value.
|
|
* \param b_m mantissa of the second input value.
|
|
* \param b_e exponent of the second input value.
|
|
* \return non-zero if (a_m*2^a_e) < (b_m*2^b_e), 0 otherwise
|
|
*/
|
|
FDK_INLINE INT fIsLessThan(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e) {
|
|
if (a_e > b_e) {
|
|
return ((b_m >> fMin(a_e - b_e, DFRACT_BITS - 1)) > a_m);
|
|
} else {
|
|
return ((a_m >> fMin(b_e - a_e, DFRACT_BITS - 1)) < b_m);
|
|
}
|
|
}
|
|
|
|
FDK_INLINE INT fIsLessThan(FIXP_SGL a_m, INT a_e, FIXP_SGL b_m, INT b_e) {
|
|
if (a_e > b_e) {
|
|
return ((b_m >> fMin(a_e - b_e, FRACT_BITS - 1)) > a_m);
|
|
} else {
|
|
return ((a_m >> fMin(b_e - a_e, FRACT_BITS - 1)) < b_m);
|
|
}
|
|
}
|
|
#endif
|
|
|
|
/**
|
|
* \brief deprecated. Use fLog2() instead.
|
|
*/
|
|
#define CalcLdData(op) fLog2(op, 0)
|
|
|
|
void LdDataVector(FIXP_DBL *srcVector, FIXP_DBL *destVector, INT number);
|
|
|
|
extern const UINT exp2_tab_long[32];
|
|
extern const UINT exp2w_tab_long[32];
|
|
extern const UINT exp2x_tab_long[32];
|
|
|
|
LNK_SECTION_CODE_L1
|
|
FDK_INLINE FIXP_DBL CalcInvLdData(const FIXP_DBL x) {
|
|
int set_zero = (x < FL2FXCONST_DBL(-31.0 / 64.0)) ? 0 : 1;
|
|
int set_max = (x >= FL2FXCONST_DBL(31.0 / 64.0)) | (x == FL2FXCONST_DBL(0.0));
|
|
|
|
FIXP_SGL frac = (FIXP_SGL)((LONG)x & 0x3FF);
|
|
UINT index3 = (UINT)(LONG)(x >> 10) & 0x1F;
|
|
UINT index2 = (UINT)(LONG)(x >> 15) & 0x1F;
|
|
UINT index1 = (UINT)(LONG)(x >> 20) & 0x1F;
|
|
int exp = fMin(31, ((x > FL2FXCONST_DBL(0.0f)) ? (31 - (int)(x >> 25))
|
|
: (int)(-(x >> 25))));
|
|
|
|
UINT lookup1 = exp2_tab_long[index1] * set_zero;
|
|
UINT lookup2 = exp2w_tab_long[index2];
|
|
UINT lookup3 = exp2x_tab_long[index3];
|
|
UINT lookup3f =
|
|
lookup3 + (UINT)(LONG)fMultDiv2((FIXP_DBL)(0x0016302F), (FIXP_SGL)frac);
|
|
|
|
UINT lookup12 = (UINT)(LONG)fMult((FIXP_DBL)lookup1, (FIXP_DBL)lookup2);
|
|
UINT lookup = (UINT)(LONG)fMult((FIXP_DBL)lookup12, (FIXP_DBL)lookup3f);
|
|
|
|
FIXP_DBL retVal = (lookup << 3) >> exp;
|
|
|
|
if (set_max) {
|
|
retVal = (FIXP_DBL)MAXVAL_DBL;
|
|
}
|
|
|
|
return retVal;
|
|
}
|
|
|
|
void InitLdInt();
|
|
FIXP_DBL CalcLdInt(INT i);
|
|
|
|
extern const USHORT sqrt_tab[49];
|
|
|
|
inline FIXP_DBL sqrtFixp_lookup(FIXP_DBL x) {
|
|
UINT y = (INT)x;
|
|
UCHAR is_zero = (y == 0);
|
|
INT zeros = fixnormz_D(y) & 0x1e;
|
|
y <<= zeros;
|
|
UINT idx = (y >> 26) - 16;
|
|
USHORT frac = (y >> 10) & 0xffff;
|
|
USHORT nfrac = 0xffff ^ frac;
|
|
UINT t = (UINT)nfrac * sqrt_tab[idx] + (UINT)frac * sqrt_tab[idx + 1];
|
|
t = t >> (zeros >> 1);
|
|
return (is_zero ? 0 : t);
|
|
}
|
|
|
|
inline FIXP_DBL sqrtFixp_lookup(FIXP_DBL x, INT *x_e) {
|
|
UINT y = (INT)x;
|
|
INT e;
|
|
|
|
if (x == (FIXP_DBL)0) {
|
|
return x;
|
|
}
|
|
|
|
/* Normalize */
|
|
e = fixnormz_D(y);
|
|
y <<= e;
|
|
e = *x_e - e + 2;
|
|
|
|
/* Correct odd exponent. */
|
|
if (e & 1) {
|
|
y >>= 1;
|
|
e++;
|
|
}
|
|
/* Get square root */
|
|
UINT idx = (y >> 26) - 16;
|
|
USHORT frac = (y >> 10) & 0xffff;
|
|
USHORT nfrac = 0xffff ^ frac;
|
|
UINT t = (UINT)nfrac * sqrt_tab[idx] + (UINT)frac * sqrt_tab[idx + 1];
|
|
|
|
/* Write back exponent */
|
|
*x_e = e >> 1;
|
|
return (FIXP_DBL)(LONG)(t >> 1);
|
|
}
|
|
|
|
void InitInvSqrtTab();
|
|
|
|
#ifndef FUNCTION_invSqrtNorm2
|
|
/**
|
|
* \brief calculate 1.0/sqrt(op)
|
|
* \param op_m mantissa of input value.
|
|
* \param result_e pointer to return the exponent of the result
|
|
* \return mantissa of the result
|
|
*/
|
|
/*****************************************************************************
|
|
delivers 1/sqrt(op) normalized to .5...1 and the shift value of the OUTPUT,
|
|
i.e. the denormalized result is 1/sqrt(op) = invSqrtNorm(op) * 2^(shift)
|
|
uses Newton-iteration for approximation
|
|
Q(n+1) = Q(n) + Q(n) * (0.5 - 2 * V * Q(n)^2)
|
|
with Q = 0.5* V ^-0.5; 0.5 <= V < 1.0
|
|
*****************************************************************************/
|
|
static FDK_FORCEINLINE FIXP_DBL invSqrtNorm2(FIXP_DBL op, INT *shift) {
|
|
FIXP_DBL val = op;
|
|
FIXP_DBL reg1, reg2;
|
|
|
|
if (val == FL2FXCONST_DBL(0.0)) {
|
|
*shift = 16;
|
|
return ((LONG)MAXVAL_DBL); /* maximum positive value */
|
|
}
|
|
|
|
#define INVSQRTNORM2_LINEAR_INTERPOLATE
|
|
#define INVSQRTNORM2_LINEAR_INTERPOLATE_HQ
|
|
|
|
/* normalize input, calculate shift value */
|
|
FDK_ASSERT(val > FL2FXCONST_DBL(0.0));
|
|
*shift = fNormz(val) - 1; /* CountLeadingBits() is not necessary here since
|
|
test value is always > 0 */
|
|
val <<= *shift; /* normalized input V */
|
|
*shift += 2; /* bias for exponent */
|
|
|
|
#if defined(INVSQRTNORM2_LINEAR_INTERPOLATE)
|
|
INT index =
|
|
(INT)(val >> (DFRACT_BITS - 1 - (SQRT_BITS + 1))) & SQRT_BITS_MASK;
|
|
FIXP_DBL Fract =
|
|
(FIXP_DBL)(((INT)val & SQRT_FRACT_BITS_MASK) << (SQRT_BITS + 1));
|
|
FIXP_DBL diff = invSqrtTab[index + 1] - invSqrtTab[index];
|
|
reg1 = invSqrtTab[index] + (fMultDiv2(diff, Fract) << 1);
|
|
#if defined(INVSQRTNORM2_LINEAR_INTERPOLATE_HQ)
|
|
/* reg1 = t[i] + (t[i+1]-t[i])*fract ... already computed ...
|
|
+ (1-fract)fract*(t[i+2]-t[i+1])/2 */
|
|
if (Fract != (FIXP_DBL)0) {
|
|
/* fract = fract * (1 - fract) */
|
|
Fract = fMultDiv2(Fract, (FIXP_DBL)((ULONG)0x80000000 - (ULONG)Fract)) << 1;
|
|
diff = diff - (invSqrtTab[index + 2] - invSqrtTab[index + 1]);
|
|
reg1 = fMultAddDiv2(reg1, Fract, diff);
|
|
}
|
|
#endif /* INVSQRTNORM2_LINEAR_INTERPOLATE_HQ */
|
|
#else
|
|
#error \
|
|
"Either define INVSQRTNORM2_NEWTON_ITERATE or INVSQRTNORM2_LINEAR_INTERPOLATE"
|
|
#endif
|
|
/* calculate the output exponent = input exp/2 */
|
|
if (*shift & 0x00000001) { /* odd shift values ? */
|
|
/* Note: Do not use rounded value 0x5A82799A to avoid overflow with
|
|
* shift-by-2 */
|
|
reg2 = (FIXP_DBL)0x5A827999;
|
|
/* FL2FXCONST_DBL(0.707106781186547524400844362104849f);*/ /* 1/sqrt(2);
|
|
*/
|
|
reg1 = fMultDiv2(reg1, reg2) << 2;
|
|
}
|
|
|
|
*shift = *shift >> 1;
|
|
|
|
return (reg1);
|
|
}
|
|
#endif /* FUNCTION_invSqrtNorm2 */
|
|
|
|
#ifndef FUNCTION_sqrtFixp
|
|
static FDK_FORCEINLINE FIXP_DBL sqrtFixp(FIXP_DBL op) {
|
|
INT tmp_exp = 0;
|
|
FIXP_DBL tmp_inv = invSqrtNorm2(op, &tmp_exp);
|
|
|
|
FDK_ASSERT(tmp_exp > 0);
|
|
return ((FIXP_DBL)(fMultDiv2((op << (tmp_exp - 1)), tmp_inv) << 2));
|
|
}
|
|
#endif /* FUNCTION_sqrtFixp */
|
|
|
|
#ifndef FUNCTION_invFixp
|
|
/**
|
|
* \brief calculate 1.0/op
|
|
* \param op mantissa of the input value.
|
|
* \return mantissa of the result with implicit exponent of 31
|
|
* \exceptions are provided for op=0,1 setting max. positive value
|
|
*/
|
|
static inline FIXP_DBL invFixp(FIXP_DBL op) {
|
|
if ((op == (FIXP_DBL)0x00000000) || (op == (FIXP_DBL)0x00000001)) {
|
|
return ((LONG)MAXVAL_DBL);
|
|
}
|
|
INT tmp_exp;
|
|
FIXP_DBL tmp_inv = invSqrtNorm2(op, &tmp_exp);
|
|
FDK_ASSERT((31 - (2 * tmp_exp + 1)) >= 0);
|
|
int shift = 31 - (2 * tmp_exp + 1);
|
|
tmp_inv = fPow2Div2(tmp_inv);
|
|
if (shift) {
|
|
tmp_inv = ((tmp_inv >> (shift - 1)) + (FIXP_DBL)1) >> 1;
|
|
}
|
|
return tmp_inv;
|
|
}
|
|
|
|
/**
|
|
* \brief calculate 1.0/(op_m * 2^op_e)
|
|
* \param op_m mantissa of the input value.
|
|
* \param op_e pointer into were the exponent of the input value is stored, and
|
|
* the result will be stored into.
|
|
* \return mantissa of the result
|
|
*/
|
|
static inline FIXP_DBL invFixp(FIXP_DBL op_m, int *op_e) {
|
|
if ((op_m == (FIXP_DBL)0x00000000) || (op_m == (FIXP_DBL)0x00000001)) {
|
|
*op_e = 31 - *op_e;
|
|
return ((LONG)MAXVAL_DBL);
|
|
}
|
|
|
|
INT tmp_exp;
|
|
FIXP_DBL tmp_inv = invSqrtNorm2(op_m, &tmp_exp);
|
|
|
|
*op_e = (tmp_exp << 1) - *op_e + 1;
|
|
return fPow2Div2(tmp_inv);
|
|
}
|
|
#endif /* FUNCTION_invFixp */
|
|
|
|
#ifndef FUNCTION_schur_div
|
|
|
|
/**
|
|
* \brief Divide two FIXP_DBL values with given precision.
|
|
* \param num dividend
|
|
* \param denum divisor
|
|
* \param count amount of significant bits of the result (starting to the MSB)
|
|
* \return num/divisor
|
|
*/
|
|
|
|
FIXP_DBL schur_div(FIXP_DBL num, FIXP_DBL denum, INT count);
|
|
|
|
#endif /* FUNCTION_schur_div */
|
|
|
|
FIXP_DBL mul_dbl_sgl_rnd(const FIXP_DBL op1, const FIXP_SGL op2);
|
|
|
|
#ifndef FUNCTION_fMultNorm
|
|
/**
|
|
* \brief multiply two values with normalization, thus max precision.
|
|
* Author: Robert Weidner
|
|
*
|
|
* \param f1 first factor
|
|
* \param f2 second factor
|
|
* \param result_e pointer to an INT where the exponent of the result is stored
|
|
* into
|
|
* \return mantissa of the product f1*f2
|
|
*/
|
|
FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2, INT *result_e);
|
|
|
|
/**
|
|
* \brief Multiply 2 values using maximum precision. The exponent of the result
|
|
* is 0.
|
|
* \param f1_m mantissa of factor 1
|
|
* \param f2_m mantissa of factor 2
|
|
* \return mantissa of the result with exponent equal to 0
|
|
*/
|
|
inline FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2) {
|
|
FIXP_DBL m;
|
|
INT e;
|
|
|
|
m = fMultNorm(f1, f2, &e);
|
|
|
|
m = scaleValueSaturate(m, e);
|
|
|
|
return m;
|
|
}
|
|
|
|
/**
|
|
* \brief Multiply 2 values with exponent and use given exponent for the
|
|
* mantissa of the result.
|
|
* \param f1_m mantissa of factor 1
|
|
* \param f1_e exponent of factor 1
|
|
* \param f2_m mantissa of factor 2
|
|
* \param f2_e exponent of factor 2
|
|
* \param result_e exponent for the returned mantissa of the result
|
|
* \return mantissa of the result with exponent equal to result_e
|
|
*/
|
|
inline FIXP_DBL fMultNorm(FIXP_DBL f1_m, INT f1_e, FIXP_DBL f2_m, INT f2_e,
|
|
INT result_e) {
|
|
FIXP_DBL m;
|
|
INT e;
|
|
|
|
m = fMultNorm(f1_m, f2_m, &e);
|
|
|
|
m = scaleValueSaturate(m, e + f1_e + f2_e - result_e);
|
|
|
|
return m;
|
|
}
|
|
#endif /* FUNCTION_fMultNorm */
|
|
|
|
#ifndef FUNCTION_fMultI
|
|
/**
|
|
* \brief Multiplies a fractional value and a integer value and performs
|
|
* rounding to nearest
|
|
* \param a fractional value
|
|
* \param b integer value
|
|
* \return integer value
|
|
*/
|
|
inline INT fMultI(FIXP_DBL a, INT b) {
|
|
FIXP_DBL m, mi;
|
|
INT m_e;
|
|
|
|
m = fMultNorm(a, (FIXP_DBL)b, &m_e);
|
|
|
|
if (m_e < (INT)0) {
|
|
if (m_e > (INT)-DFRACT_BITS) {
|
|
m = m >> ((-m_e) - 1);
|
|
mi = (m + (FIXP_DBL)1) >> 1;
|
|
} else {
|
|
mi = (FIXP_DBL)0;
|
|
}
|
|
} else {
|
|
mi = scaleValueSaturate(m, m_e);
|
|
}
|
|
|
|
return ((INT)mi);
|
|
}
|
|
#endif /* FUNCTION_fMultI */
|
|
|
|
#ifndef FUNCTION_fMultIfloor
|
|
/**
|
|
* \brief Multiplies a fractional value and a integer value and performs floor
|
|
* rounding
|
|
* \param a fractional value
|
|
* \param b integer value
|
|
* \return integer value
|
|
*/
|
|
inline INT fMultIfloor(FIXP_DBL a, INT b) {
|
|
FIXP_DBL m, mi;
|
|
INT m_e;
|
|
|
|
m = fMultNorm(a, (FIXP_DBL)b, &m_e);
|
|
|
|
if (m_e < (INT)0) {
|
|
if (m_e > (INT)-DFRACT_BITS) {
|
|
mi = m >> (-m_e);
|
|
} else {
|
|
mi = (FIXP_DBL)0;
|
|
if (m < (FIXP_DBL)0) {
|
|
mi = (FIXP_DBL)-1;
|
|
}
|
|
}
|
|
} else {
|
|
mi = scaleValueSaturate(m, m_e);
|
|
}
|
|
|
|
return ((INT)mi);
|
|
}
|
|
#endif /* FUNCTION_fMultIfloor */
|
|
|
|
#ifndef FUNCTION_fMultIceil
|
|
/**
|
|
* \brief Multiplies a fractional value and a integer value and performs ceil
|
|
* rounding
|
|
* \param a fractional value
|
|
* \param b integer value
|
|
* \return integer value
|
|
*/
|
|
inline INT fMultIceil(FIXP_DBL a, INT b) {
|
|
FIXP_DBL m, mi;
|
|
INT m_e;
|
|
|
|
m = fMultNorm(a, (FIXP_DBL)b, &m_e);
|
|
|
|
if (m_e < (INT)0) {
|
|
if (m_e > (INT)-DFRACT_BITS) {
|
|
mi = (m >> (-m_e));
|
|
if ((LONG)m & ((1 << (-m_e)) - 1)) {
|
|
mi = mi + (FIXP_DBL)1;
|
|
}
|
|
} else {
|
|
mi = (FIXP_DBL)1;
|
|
if (m < (FIXP_DBL)0) {
|
|
mi = (FIXP_DBL)0;
|
|
}
|
|
}
|
|
} else {
|
|
mi = scaleValueSaturate(m, m_e);
|
|
}
|
|
|
|
return ((INT)mi);
|
|
}
|
|
#endif /* FUNCTION_fMultIceil */
|
|
|
|
#ifndef FUNCTION_fDivNorm
|
|
/**
|
|
* \brief Divide 2 FIXP_DBL values with normalization of input values.
|
|
* \param num numerator
|
|
* \param denum denominator
|
|
* \param result_e pointer to an INT where the exponent of the result is stored
|
|
* into
|
|
* \return num/denum with exponent = *result_e
|
|
*/
|
|
FIXP_DBL fDivNorm(FIXP_DBL num, FIXP_DBL denom, INT *result_e);
|
|
|
|
/**
|
|
* \brief Divide 2 positive FIXP_DBL values with normalization of input values.
|
|
* \param num numerator
|
|
* \param denum denominator
|
|
* \return num/denum with exponent = 0
|
|
*/
|
|
FIXP_DBL fDivNorm(FIXP_DBL num, FIXP_DBL denom);
|
|
|
|
/**
|
|
* \brief Divide 2 signed FIXP_DBL values with normalization of input values.
|
|
* \param num numerator
|
|
* \param denum denominator
|
|
* \param result_e pointer to an INT where the exponent of the result is stored
|
|
* into
|
|
* \return num/denum with exponent = *result_e
|
|
*/
|
|
FIXP_DBL fDivNormSigned(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e);
|
|
|
|
/**
|
|
* \brief Divide 2 signed FIXP_DBL values with normalization of input values.
|
|
* \param num numerator
|
|
* \param denum denominator
|
|
* \return num/denum with exponent = 0
|
|
*/
|
|
FIXP_DBL fDivNormSigned(FIXP_DBL num, FIXP_DBL denom);
|
|
#endif /* FUNCTION_fDivNorm */
|
|
|
|
/**
|
|
* \brief Adjust mantissa to exponent -1
|
|
* \param a_m mantissa of value to be adjusted
|
|
* \param pA_e pointer to the exponen of a_m
|
|
* \return adjusted mantissa
|
|
*/
|
|
inline FIXP_DBL fAdjust(FIXP_DBL a_m, INT *pA_e) {
|
|
INT shift;
|
|
|
|
shift = fNorm(a_m) - 1;
|
|
*pA_e -= shift;
|
|
|
|
return scaleValue(a_m, shift);
|
|
}
|
|
|
|
#ifndef FUNCTION_fAddNorm
|
|
/**
|
|
* \brief Add two values with normalization
|
|
* \param a_m mantissa of first summand
|
|
* \param a_e exponent of first summand
|
|
* \param a_m mantissa of second summand
|
|
* \param a_e exponent of second summand
|
|
* \param pResult_e pointer to where the exponent of the result will be stored
|
|
* to.
|
|
* \return mantissa of result
|
|
*/
|
|
inline FIXP_DBL fAddNorm(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e,
|
|
INT *pResult_e) {
|
|
INT result_e;
|
|
FIXP_DBL result_m;
|
|
|
|
/* If one of the summands is zero, return the other.
|
|
This is necessary for the summation of a very small number to zero */
|
|
if (a_m == (FIXP_DBL)0) {
|
|
*pResult_e = b_e;
|
|
return b_m;
|
|
}
|
|
if (b_m == (FIXP_DBL)0) {
|
|
*pResult_e = a_e;
|
|
return a_m;
|
|
}
|
|
|
|
a_m = fAdjust(a_m, &a_e);
|
|
b_m = fAdjust(b_m, &b_e);
|
|
|
|
if (a_e > b_e) {
|
|
result_m = a_m + (b_m >> fMin(a_e - b_e, DFRACT_BITS - 1));
|
|
result_e = a_e;
|
|
} else {
|
|
result_m = (a_m >> fMin(b_e - a_e, DFRACT_BITS - 1)) + b_m;
|
|
result_e = b_e;
|
|
}
|
|
|
|
*pResult_e = result_e;
|
|
return result_m;
|
|
}
|
|
|
|
inline FIXP_DBL fAddNorm(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e,
|
|
INT result_e) {
|
|
FIXP_DBL result_m;
|
|
|
|
a_m = scaleValue(a_m, a_e - result_e);
|
|
b_m = scaleValue(b_m, b_e - result_e);
|
|
|
|
result_m = a_m + b_m;
|
|
|
|
return result_m;
|
|
}
|
|
#endif /* FUNCTION_fAddNorm */
|
|
|
|
/**
|
|
* \brief Divide 2 FIXP_DBL values with normalization of input values.
|
|
* \param num numerator
|
|
* \param denum denomintator
|
|
* \return num/denum with exponent = 0
|
|
*/
|
|
FIXP_DBL fDivNormHighPrec(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e);
|
|
|
|
#ifndef FUNCTION_fPow
|
|
/**
|
|
* \brief return 2 ^ (exp_m * 2^exp_e)
|
|
* \param exp_m mantissa of the exponent to 2.0f
|
|
* \param exp_e exponent of the exponent to 2.0f
|
|
* \param result_e pointer to a INT where the exponent of the result will be
|
|
* stored into
|
|
* \return mantissa of the result
|
|
*/
|
|
FIXP_DBL f2Pow(const FIXP_DBL exp_m, const INT exp_e, INT *result_e);
|
|
|
|
/**
|
|
* \brief return 2 ^ (exp_m * 2^exp_e). This version returns only the mantissa
|
|
* with implicit exponent of zero.
|
|
* \param exp_m mantissa of the exponent to 2.0f
|
|
* \param exp_e exponent of the exponent to 2.0f
|
|
* \return mantissa of the result
|
|
*/
|
|
FIXP_DBL f2Pow(const FIXP_DBL exp_m, const INT exp_e);
|
|
|
|
/**
|
|
* \brief return x ^ (exp_m * 2^exp_e), where log2(x) = baseLd_m * 2^(baseLd_e).
|
|
* This saves the need to compute log2() of constant values (when x is a
|
|
* constant).
|
|
* \param baseLd_m mantissa of log2() of x.
|
|
* \param baseLd_e exponent of log2() of x.
|
|
* \param exp_m mantissa of the exponent to 2.0f
|
|
* \param exp_e exponent of the exponent to 2.0f
|
|
* \param result_e pointer to a INT where the exponent of the result will be
|
|
* stored into
|
|
* \return mantissa of the result
|
|
*/
|
|
FIXP_DBL fLdPow(FIXP_DBL baseLd_m, INT baseLd_e, FIXP_DBL exp_m, INT exp_e,
|
|
INT *result_e);
|
|
|
|
/**
|
|
* \brief return x ^ (exp_m * 2^exp_e), where log2(x) = baseLd_m * 2^(baseLd_e).
|
|
* This saves the need to compute log2() of constant values (when x is a
|
|
* constant). This version does not return an exponent, which is
|
|
* implicitly 0.
|
|
* \param baseLd_m mantissa of log2() of x.
|
|
* \param baseLd_e exponent of log2() of x.
|
|
* \param exp_m mantissa of the exponent to 2.0f
|
|
* \param exp_e exponent of the exponent to 2.0f
|
|
* \return mantissa of the result
|
|
*/
|
|
FIXP_DBL fLdPow(FIXP_DBL baseLd_m, INT baseLd_e, FIXP_DBL exp_m, INT exp_e);
|
|
|
|
/**
|
|
* \brief return (base_m * 2^base_e) ^ (exp * 2^exp_e). Use fLdPow() instead
|
|
* whenever possible.
|
|
* \param base_m mantissa of the base.
|
|
* \param base_e exponent of the base.
|
|
* \param exp_m mantissa of power to be calculated of the base.
|
|
* \param exp_e exponent of power to be calculated of the base.
|
|
* \param result_e pointer to a INT where the exponent of the result will be
|
|
* stored into.
|
|
* \return mantissa of the result.
|
|
*/
|
|
FIXP_DBL fPow(FIXP_DBL base_m, INT base_e, FIXP_DBL exp_m, INT exp_e,
|
|
INT *result_e);
|
|
|
|
/**
|
|
* \brief return (base_m * 2^base_e) ^ N
|
|
* \param base_m mantissa of the base
|
|
* \param base_e exponent of the base
|
|
* \param N power to be calculated of the base
|
|
* \param result_e pointer to a INT where the exponent of the result will be
|
|
* stored into
|
|
* \return mantissa of the result
|
|
*/
|
|
FIXP_DBL fPowInt(FIXP_DBL base_m, INT base_e, INT N, INT *result_e);
|
|
#endif /* #ifndef FUNCTION_fPow */
|
|
|
|
#ifndef FUNCTION_fLog2
|
|
/**
|
|
* \brief Calculate log(argument)/log(2) (logarithm with base 2). deprecated.
|
|
* Use fLog2() instead.
|
|
* \param arg mantissa of the argument
|
|
* \param arg_e exponent of the argument
|
|
* \param result_e pointer to an INT to store the exponent of the result
|
|
* \return the mantissa of the result.
|
|
* \param
|
|
*/
|
|
FIXP_DBL CalcLog2(FIXP_DBL arg, INT arg_e, INT *result_e);
|
|
|
|
/**
|
|
* \brief calculate logarithm of base 2 of x_m * 2^(x_e)
|
|
* \param x_m mantissa of the input value.
|
|
* \param x_e exponent of the input value.
|
|
* \param pointer to an INT where the exponent of the result is returned into.
|
|
* \return mantissa of the result.
|
|
*/
|
|
FDK_INLINE FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e, INT *result_e) {
|
|
FIXP_DBL result_m;
|
|
|
|
/* Short cut for zero and negative numbers. */
|
|
if (x_m <= FL2FXCONST_DBL(0.0f)) {
|
|
*result_e = DFRACT_BITS - 1;
|
|
return FL2FXCONST_DBL(-1.0f);
|
|
}
|
|
|
|
/* Calculate log2() */
|
|
{
|
|
FIXP_DBL x2_m;
|
|
|
|
/* Move input value x_m * 2^x_e toward 1.0, where the taylor approximation
|
|
of the function log(1-x) centered at 0 is most accurate. */
|
|
{
|
|
INT b_norm;
|
|
|
|
b_norm = fNormz(x_m) - 1;
|
|
x2_m = x_m << b_norm;
|
|
x_e = x_e - b_norm;
|
|
}
|
|
|
|
/* map x from log(x) domain to log(1-x) domain. */
|
|
x2_m = -(x2_m + FL2FXCONST_DBL(-1.0));
|
|
|
|
/* Taylor polynomial approximation of ln(1-x) */
|
|
{
|
|
FIXP_DBL px2_m;
|
|
result_m = FL2FXCONST_DBL(0.0);
|
|
px2_m = x2_m;
|
|
for (int i = 0; i < LD_PRECISION; i++) {
|
|
result_m = fMultAddDiv2(result_m, ldCoeff[i], px2_m);
|
|
px2_m = fMult(px2_m, x2_m);
|
|
}
|
|
}
|
|
/* Multiply result with 1/ln(2) = 1.0 + 0.442695040888 (get log2(x) from
|
|
* ln(x) result). */
|
|
result_m =
|
|
fMultAddDiv2(result_m, result_m,
|
|
FL2FXCONST_DBL(2.0 * 0.4426950408889634073599246810019));
|
|
|
|
/* Add exponent part. log2(x_m * 2^x_e) = log2(x_m) + x_e */
|
|
if (x_e != 0) {
|
|
int enorm;
|
|
|
|
enorm = DFRACT_BITS - fNorm((FIXP_DBL)x_e);
|
|
/* The -1 in the right shift of result_m compensates the fMultDiv2() above
|
|
* in the taylor polynomial evaluation loop.*/
|
|
result_m = (result_m >> (enorm - 1)) +
|
|
((FIXP_DBL)x_e << (DFRACT_BITS - 1 - enorm));
|
|
|
|
*result_e = enorm;
|
|
} else {
|
|
/* 1 compensates the fMultDiv2() above in the taylor polynomial evaluation
|
|
* loop.*/
|
|
*result_e = 1;
|
|
}
|
|
}
|
|
|
|
return result_m;
|
|
}
|
|
|
|
/**
|
|
* \brief calculate logarithm of base 2 of x_m * 2^(x_e)
|
|
* \param x_m mantissa of the input value.
|
|
* \param x_e exponent of the input value.
|
|
* \return mantissa of the result with implicit exponent of LD_DATA_SHIFT.
|
|
*/
|
|
FDK_INLINE FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e) {
|
|
if (x_m <= FL2FXCONST_DBL(0.0f)) {
|
|
x_m = FL2FXCONST_DBL(-1.0f);
|
|
} else {
|
|
INT result_e;
|
|
x_m = fLog2(x_m, x_e, &result_e);
|
|
x_m = scaleValue(x_m, result_e - LD_DATA_SHIFT);
|
|
}
|
|
return x_m;
|
|
}
|
|
|
|
#endif /* FUNCTION_fLog2 */
|
|
|
|
#ifndef FUNCTION_fAddSaturate
|
|
/**
|
|
* \brief Add with saturation of the result.
|
|
* \param a first summand
|
|
* \param b second summand
|
|
* \return saturated sum of a and b.
|
|
*/
|
|
inline FIXP_SGL fAddSaturate(const FIXP_SGL a, const FIXP_SGL b) {
|
|
LONG sum;
|
|
|
|
sum = (LONG)(SHORT)a + (LONG)(SHORT)b;
|
|
sum = fMax(fMin((INT)sum, (INT)MAXVAL_SGL), (INT)MINVAL_SGL);
|
|
return (FIXP_SGL)(SHORT)sum;
|
|
}
|
|
|
|
/**
|
|
* \brief Add with saturation of the result.
|
|
* \param a first summand
|
|
* \param b second summand
|
|
* \return saturated sum of a and b.
|
|
*/
|
|
inline FIXP_DBL fAddSaturate(const FIXP_DBL a, const FIXP_DBL b) {
|
|
LONG sum;
|
|
|
|
sum = (LONG)(a >> 1) + (LONG)(b >> 1);
|
|
sum = fMax(fMin((INT)sum, (INT)(MAXVAL_DBL >> 1)), (INT)(MINVAL_DBL >> 1));
|
|
return (FIXP_DBL)(LONG)(sum << 1);
|
|
}
|
|
#endif /* FUNCTION_fAddSaturate */
|
|
|
|
INT fixp_floorToInt(FIXP_DBL f_inp, INT sf);
|
|
FIXP_DBL fixp_floor(FIXP_DBL f_inp, INT sf);
|
|
|
|
INT fixp_ceilToInt(FIXP_DBL f_inp, INT sf);
|
|
FIXP_DBL fixp_ceil(FIXP_DBL f_inp, INT sf);
|
|
|
|
INT fixp_truncateToInt(FIXP_DBL f_inp, INT sf);
|
|
FIXP_DBL fixp_truncate(FIXP_DBL f_inp, INT sf);
|
|
|
|
INT fixp_roundToInt(FIXP_DBL f_inp, INT sf);
|
|
FIXP_DBL fixp_round(FIXP_DBL f_inp, INT sf);
|
|
|
|
/*****************************************************************************
|
|
|
|
array for 1/n, n=1..80
|
|
|
|
****************************************************************************/
|
|
|
|
extern const FIXP_DBL invCount[80];
|
|
|
|
LNK_SECTION_INITCODE
|
|
inline void InitInvInt(void) {}
|
|
|
|
/**
|
|
* \brief Calculate the value of 1/i where i is a integer value. It supports
|
|
* input values from 1 upto (80-1).
|
|
* \param intValue Integer input value.
|
|
* \param FIXP_DBL representation of 1/intValue
|
|
*/
|
|
inline FIXP_DBL GetInvInt(int intValue) {
|
|
return invCount[fMin(fMax(intValue, 0), 80 - 1)];
|
|
}
|
|
|
|
#endif /* FIXPOINT_MATH_H */
|