1 /* ----------------------------------------------------------------------------- 2 Software License for The Fraunhofer FDK AAC Codec Library for Android 3 4 Copyright 1995 - 2018 Fraunhofer-Gesellschaft zur Frderung der angewandten 5 Forschung e.V. All rights reserved. 6 7 1. INTRODUCTION 8 The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software 9 that implements the MPEG Advanced Audio Coding ("AAC") encoding and decoding 10 scheme for digital audio. This FDK AAC Codec software is intended to be used on 11 a wide variety of Android devices. 12 13 AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient 14 general perceptual audio codecs. AAC-ELD is considered the best-performing 15 full-bandwidth communications codec by independent studies and is widely 16 deployed. AAC has been standardized by ISO and IEC as part of the MPEG 17 specifications. 18 19 Patent licenses for necessary patent claims for the FDK AAC Codec (including 20 those of Fraunhofer) may be obtained through Via Licensing 21 (www.vialicensing.com) or through the respective patent owners individually for 22 the purpose of encoding or decoding bit streams in products that are compliant 23 with the ISO/IEC MPEG audio standards. Please note that most manufacturers of 24 Android devices already license these patent claims through Via Licensing or 25 directly from the patent owners, and therefore FDK AAC Codec software may 26 already be covered under those patent licenses when it is used for those 27 licensed purposes only. 28 29 Commercially-licensed AAC software libraries, including floating-point versions 30 with enhanced sound quality, are also available from Fraunhofer. Users are 31 encouraged to check the Fraunhofer website for additional applications 32 information and documentation. 33 34 2. COPYRIGHT LICENSE 35 36 Redistribution and use in source and binary forms, with or without modification, 37 are permitted without payment of copyright license fees provided that you 38 satisfy the following conditions: 39 40 You must retain the complete text of this software license in redistributions of 41 the FDK AAC Codec or your modifications thereto in source code form. 42 43 You must retain the complete text of this software license in the documentation 44 and/or other materials provided with redistributions of the FDK AAC Codec or 45 your modifications thereto in binary form. You must make available free of 46 charge copies of the complete source code of the FDK AAC Codec and your 47 modifications thereto to recipients of copies in binary form. 48 49 The name of Fraunhofer may not be used to endorse or promote products derived 50 from this library without prior written permission. 51 52 You may not charge copyright license fees for anyone to use, copy or distribute 53 the FDK AAC Codec software or your modifications thereto. 54 55 Your modified versions of the FDK AAC Codec must carry prominent notices stating 56 that you changed the software and the date of any change. For modified versions 57 of the FDK AAC Codec, the term "Fraunhofer FDK AAC Codec Library for Android" 58 must be replaced by the term "Third-Party Modified Version of the Fraunhofer FDK 59 AAC Codec Library for Android." 60 61 3. NO PATENT LICENSE 62 63 NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without 64 limitation the patents of Fraunhofer, ARE GRANTED BY THIS SOFTWARE LICENSE. 65 Fraunhofer provides no warranty of patent non-infringement with respect to this 66 software. 67 68 You may use this FDK AAC Codec software or modifications thereto only for 69 purposes that are authorized by appropriate patent licenses. 70 71 4. DISCLAIMER 72 73 This FDK AAC Codec software is provided by Fraunhofer on behalf of the copyright 74 holders and contributors "AS IS" and WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES, 75 including but not limited to the implied warranties of merchantability and 76 fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR 77 CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary, 78 or consequential damages, including but not limited to procurement of substitute 79 goods or services; loss of use, data, or profits, or business interruption, 80 however caused and on any theory of liability, whether in contract, strict 81 liability, or tort (including negligence), arising in any way out of the use of 82 this software, even if advised of the possibility of such damage. 83 84 5. CONTACT INFORMATION 85 86 Fraunhofer Institute for Integrated Circuits IIS 87 Attention: Audio and Multimedia Departments - FDK AAC LL 88 Am Wolfsmantel 33 89 91058 Erlangen, Germany 90 91 www.iis.fraunhofer.de/amm 92 amm-info (at) iis.fraunhofer.de 93 ----------------------------------------------------------------------------- */ 94 95 /******************* Library for basic calculation routines ******************** 96 97 Author(s): M. Gayer 98 99 Description: Fixed point specific mathematical functions 100 101 *******************************************************************************/ 102 103 #ifndef FIXPOINT_MATH_H 104 #define FIXPOINT_MATH_H 105 106 #include "common_fix.h" 107 #include "scale.h" 108 109 /* 110 * Data definitions 111 */ 112 113 #define LD_DATA_SCALING (64.0f) 114 #define LD_DATA_SHIFT 6 /* pow(2, LD_DATA_SHIFT) = LD_DATA_SCALING */ 115 116 #define MAX_LD_PRECISION 10 117 #define LD_PRECISION 10 118 119 /* Taylor series coefficients for ln(1-x), centered at 0 (MacLaurin polynomial). 120 */ 121 #ifndef LDCOEFF_16BIT 122 LNK_SECTION_CONSTDATA_L1 123 static const FIXP_DBL ldCoeff[MAX_LD_PRECISION] = { 124 FL2FXCONST_DBL(-1.0), FL2FXCONST_DBL(-1.0 / 2.0), 125 FL2FXCONST_DBL(-1.0 / 3.0), FL2FXCONST_DBL(-1.0 / 4.0), 126 FL2FXCONST_DBL(-1.0 / 5.0), FL2FXCONST_DBL(-1.0 / 6.0), 127 FL2FXCONST_DBL(-1.0 / 7.0), FL2FXCONST_DBL(-1.0 / 8.0), 128 FL2FXCONST_DBL(-1.0 / 9.0), FL2FXCONST_DBL(-1.0 / 10.0)}; 129 #else /* LDCOEFF_16BIT */ 130 LNK_SECTION_CONSTDATA_L1 131 static const FIXP_SGL ldCoeff[MAX_LD_PRECISION] = { 132 FL2FXCONST_SGL(-1.0), FL2FXCONST_SGL(-1.0 / 2.0), 133 FL2FXCONST_SGL(-1.0 / 3.0), FL2FXCONST_SGL(-1.0 / 4.0), 134 FL2FXCONST_SGL(-1.0 / 5.0), FL2FXCONST_SGL(-1.0 / 6.0), 135 FL2FXCONST_SGL(-1.0 / 7.0), FL2FXCONST_SGL(-1.0 / 8.0), 136 FL2FXCONST_SGL(-1.0 / 9.0), FL2FXCONST_SGL(-1.0 / 10.0)}; 137 #endif /* LDCOEFF_16BIT */ 138 139 /***************************************************************************** 140 141 functionname: invSqrtNorm2 142 description: delivers 1/sqrt(op) normalized to .5...1 and the shift value 143 of the OUTPUT 144 145 *****************************************************************************/ 146 #define SQRT_BITS 7 147 #define SQRT_VALUES (128 + 2) 148 #define SQRT_BITS_MASK 0x7f 149 #define SQRT_FRACT_BITS_MASK 0x007FFFFF 150 151 extern const FIXP_DBL invSqrtTab[SQRT_VALUES]; 152 153 /* 154 * Hardware specific implementations 155 */ 156 157 #if defined(__x86__) 158 #include "x86/fixpoint_math_x86.h" 159 #endif /* target architecture selector */ 160 161 /* 162 * Fallback implementations 163 */ 164 #if !defined(FUNCTION_fIsLessThan) 165 /** 166 * \brief Compares two fixpoint values incl. scaling. 167 * \param a_m mantissa of the first input value. 168 * \param a_e exponent of the first input value. 169 * \param b_m mantissa of the second input value. 170 * \param b_e exponent of the second input value. 171 * \return non-zero if (a_m*2^a_e) < (b_m*2^b_e), 0 otherwise 172 */ 173 FDK_INLINE INT fIsLessThan(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e) { 174 if (a_e > b_e) { 175 return ((b_m >> fMin(a_e - b_e, DFRACT_BITS - 1)) > a_m); 176 } else { 177 return ((a_m >> fMin(b_e - a_e, DFRACT_BITS - 1)) < b_m); 178 } 179 } 180 181 FDK_INLINE INT fIsLessThan(FIXP_SGL a_m, INT a_e, FIXP_SGL b_m, INT b_e) { 182 if (a_e > b_e) { 183 return ((b_m >> fMin(a_e - b_e, FRACT_BITS - 1)) > a_m); 184 } else { 185 return ((a_m >> fMin(b_e - a_e, FRACT_BITS - 1)) < b_m); 186 } 187 } 188 #endif 189 190 /** 191 * \brief deprecated. Use fLog2() instead. 192 */ 193 #define CalcLdData(op) fLog2(op, 0) 194 195 void LdDataVector(FIXP_DBL *srcVector, FIXP_DBL *destVector, INT number); 196 197 extern const UINT exp2_tab_long[32]; 198 extern const UINT exp2w_tab_long[32]; 199 extern const UINT exp2x_tab_long[32]; 200 201 LNK_SECTION_CODE_L1 202 FDK_INLINE FIXP_DBL CalcInvLdData(const FIXP_DBL x) { 203 int set_zero = (x < FL2FXCONST_DBL(-31.0 / 64.0)) ? 0 : 1; 204 int set_max = (x >= FL2FXCONST_DBL(31.0 / 64.0)) | (x == FL2FXCONST_DBL(0.0)); 205 206 FIXP_SGL frac = (FIXP_SGL)((LONG)x & 0x3FF); 207 UINT index3 = (UINT)(LONG)(x >> 10) & 0x1F; 208 UINT index2 = (UINT)(LONG)(x >> 15) & 0x1F; 209 UINT index1 = (UINT)(LONG)(x >> 20) & 0x1F; 210 int exp = fMin(31, ((x > FL2FXCONST_DBL(0.0f)) ? (31 - (int)(x >> 25)) 211 : (int)(-(x >> 25)))); 212 213 UINT lookup1 = exp2_tab_long[index1] * set_zero; 214 UINT lookup2 = exp2w_tab_long[index2]; 215 UINT lookup3 = exp2x_tab_long[index3]; 216 UINT lookup3f = 217 lookup3 + (UINT)(LONG)fMultDiv2((FIXP_DBL)(0x0016302F), (FIXP_SGL)frac); 218 219 UINT lookup12 = (UINT)(LONG)fMult((FIXP_DBL)lookup1, (FIXP_DBL)lookup2); 220 UINT lookup = (UINT)(LONG)fMult((FIXP_DBL)lookup12, (FIXP_DBL)lookup3f); 221 222 FIXP_DBL retVal = (lookup << 3) >> exp; 223 224 if (set_max) { 225 retVal = (FIXP_DBL)MAXVAL_DBL; 226 } 227 228 return retVal; 229 } 230 231 void InitLdInt(); 232 FIXP_DBL CalcLdInt(INT i); 233 234 extern const USHORT sqrt_tab[49]; 235 236 inline FIXP_DBL sqrtFixp_lookup(FIXP_DBL x) { 237 UINT y = (INT)x; 238 UCHAR is_zero = (y == 0); 239 INT zeros = fixnormz_D(y) & 0x1e; 240 y <<= zeros; 241 UINT idx = (y >> 26) - 16; 242 USHORT frac = (y >> 10) & 0xffff; 243 USHORT nfrac = 0xffff ^ frac; 244 UINT t = (UINT)nfrac * sqrt_tab[idx] + (UINT)frac * sqrt_tab[idx + 1]; 245 t = t >> (zeros >> 1); 246 return (is_zero ? 0 : t); 247 } 248 249 inline FIXP_DBL sqrtFixp_lookup(FIXP_DBL x, INT *x_e) { 250 UINT y = (INT)x; 251 INT e; 252 253 if (x == (FIXP_DBL)0) { 254 return x; 255 } 256 257 /* Normalize */ 258 e = fixnormz_D(y); 259 y <<= e; 260 e = *x_e - e + 2; 261 262 /* Correct odd exponent. */ 263 if (e & 1) { 264 y >>= 1; 265 e++; 266 } 267 /* Get square root */ 268 UINT idx = (y >> 26) - 16; 269 USHORT frac = (y >> 10) & 0xffff; 270 USHORT nfrac = 0xffff ^ frac; 271 UINT t = (UINT)nfrac * sqrt_tab[idx] + (UINT)frac * sqrt_tab[idx + 1]; 272 273 /* Write back exponent */ 274 *x_e = e >> 1; 275 return (FIXP_DBL)(LONG)(t >> 1); 276 } 277 278 void InitInvSqrtTab(); 279 280 #ifndef FUNCTION_invSqrtNorm2 281 /** 282 * \brief calculate 1.0/sqrt(op) 283 * \param op_m mantissa of input value. 284 * \param result_e pointer to return the exponent of the result 285 * \return mantissa of the result 286 */ 287 /***************************************************************************** 288 delivers 1/sqrt(op) normalized to .5...1 and the shift value of the OUTPUT, 289 i.e. the denormalized result is 1/sqrt(op) = invSqrtNorm(op) * 2^(shift) 290 uses Newton-iteration for approximation 291 Q(n+1) = Q(n) + Q(n) * (0.5 - 2 * V * Q(n)^2) 292 with Q = 0.5* V ^-0.5; 0.5 <= V < 1.0 293 *****************************************************************************/ 294 static FDK_FORCEINLINE FIXP_DBL invSqrtNorm2(FIXP_DBL op, INT *shift) { 295 FIXP_DBL val = op; 296 FIXP_DBL reg1, reg2; 297 298 if (val == FL2FXCONST_DBL(0.0)) { 299 *shift = 16; 300 return ((LONG)MAXVAL_DBL); /* maximum positive value */ 301 } 302 303 #define INVSQRTNORM2_LINEAR_INTERPOLATE 304 #define INVSQRTNORM2_LINEAR_INTERPOLATE_HQ 305 306 /* normalize input, calculate shift value */ 307 FDK_ASSERT(val > FL2FXCONST_DBL(0.0)); 308 *shift = fNormz(val) - 1; /* CountLeadingBits() is not necessary here since 309 test value is always > 0 */ 310 val <<= *shift; /* normalized input V */ 311 *shift += 2; /* bias for exponent */ 312 313 #if defined(INVSQRTNORM2_LINEAR_INTERPOLATE) 314 INT index = 315 (INT)(val >> (DFRACT_BITS - 1 - (SQRT_BITS + 1))) & SQRT_BITS_MASK; 316 FIXP_DBL Fract = 317 (FIXP_DBL)(((INT)val & SQRT_FRACT_BITS_MASK) << (SQRT_BITS + 1)); 318 FIXP_DBL diff = invSqrtTab[index + 1] - invSqrtTab[index]; 319 reg1 = invSqrtTab[index] + (fMultDiv2(diff, Fract) << 1); 320 #if defined(INVSQRTNORM2_LINEAR_INTERPOLATE_HQ) 321 /* reg1 = t[i] + (t[i+1]-t[i])*fract ... already computed ... 322 + (1-fract)fract*(t[i+2]-t[i+1])/2 */ 323 if (Fract != (FIXP_DBL)0) { 324 /* fract = fract * (1 - fract) */ 325 Fract = fMultDiv2(Fract, (FIXP_DBL)((ULONG)0x80000000 - (ULONG)Fract)) << 1; 326 diff = diff - (invSqrtTab[index + 2] - invSqrtTab[index + 1]); 327 reg1 = fMultAddDiv2(reg1, Fract, diff); 328 } 329 #endif /* INVSQRTNORM2_LINEAR_INTERPOLATE_HQ */ 330 #else 331 #error \ 332 "Either define INVSQRTNORM2_NEWTON_ITERATE or INVSQRTNORM2_LINEAR_INTERPOLATE" 333 #endif 334 /* calculate the output exponent = input exp/2 */ 335 if (*shift & 0x00000001) { /* odd shift values ? */ 336 /* Note: Do not use rounded value 0x5A82799A to avoid overflow with 337 * shift-by-2 */ 338 reg2 = (FIXP_DBL)0x5A827999; 339 /* FL2FXCONST_DBL(0.707106781186547524400844362104849f);*/ /* 1/sqrt(2); 340 */ 341 reg1 = fMultDiv2(reg1, reg2) << 2; 342 } 343 344 *shift = *shift >> 1; 345 346 return (reg1); 347 } 348 #endif /* FUNCTION_invSqrtNorm2 */ 349 350 #ifndef FUNCTION_sqrtFixp 351 static FDK_FORCEINLINE FIXP_DBL sqrtFixp(FIXP_DBL op) { 352 INT tmp_exp = 0; 353 FIXP_DBL tmp_inv = invSqrtNorm2(op, &tmp_exp); 354 355 FDK_ASSERT(tmp_exp > 0); 356 return ((FIXP_DBL)(fMultDiv2((op << (tmp_exp - 1)), tmp_inv) << 2)); 357 } 358 #endif /* FUNCTION_sqrtFixp */ 359 360 #ifndef FUNCTION_invFixp 361 /** 362 * \brief calculate 1.0/op 363 * \param op mantissa of the input value. 364 * \return mantissa of the result with implicit exponent of 31 365 * \exceptions are provided for op=0,1 setting max. positive value 366 */ 367 static inline FIXP_DBL invFixp(FIXP_DBL op) { 368 if ((op == (FIXP_DBL)0x00000000) || (op == (FIXP_DBL)0x00000001)) { 369 return ((LONG)MAXVAL_DBL); 370 } 371 INT tmp_exp; 372 FIXP_DBL tmp_inv = invSqrtNorm2(op, &tmp_exp); 373 FDK_ASSERT((31 - (2 * tmp_exp + 1)) >= 0); 374 int shift = 31 - (2 * tmp_exp + 1); 375 tmp_inv = fPow2Div2(tmp_inv); 376 if (shift) { 377 tmp_inv = ((tmp_inv >> (shift - 1)) + (FIXP_DBL)1) >> 1; 378 } 379 return tmp_inv; 380 } 381 382 /** 383 * \brief calculate 1.0/(op_m * 2^op_e) 384 * \param op_m mantissa of the input value. 385 * \param op_e pointer into were the exponent of the input value is stored, and 386 * the result will be stored into. 387 * \return mantissa of the result 388 */ 389 static inline FIXP_DBL invFixp(FIXP_DBL op_m, int *op_e) { 390 if ((op_m == (FIXP_DBL)0x00000000) || (op_m == (FIXP_DBL)0x00000001)) { 391 *op_e = 31 - *op_e; 392 return ((LONG)MAXVAL_DBL); 393 } 394 395 INT tmp_exp; 396 FIXP_DBL tmp_inv = invSqrtNorm2(op_m, &tmp_exp); 397 398 *op_e = (tmp_exp << 1) - *op_e + 1; 399 return fPow2Div2(tmp_inv); 400 } 401 #endif /* FUNCTION_invFixp */ 402 403 #ifndef FUNCTION_schur_div 404 405 /** 406 * \brief Divide two FIXP_DBL values with given precision. 407 * \param num dividend 408 * \param denum divisor 409 * \param count amount of significant bits of the result (starting to the MSB) 410 * \return num/divisor 411 */ 412 413 FIXP_DBL schur_div(FIXP_DBL num, FIXP_DBL denum, INT count); 414 415 #endif /* FUNCTION_schur_div */ 416 417 FIXP_DBL mul_dbl_sgl_rnd(const FIXP_DBL op1, const FIXP_SGL op2); 418 419 #ifndef FUNCTION_fMultNorm 420 /** 421 * \brief multiply two values with normalization, thus max precision. 422 * Author: Robert Weidner 423 * 424 * \param f1 first factor 425 * \param f2 second factor 426 * \param result_e pointer to an INT where the exponent of the result is stored 427 * into 428 * \return mantissa of the product f1*f2 429 */ 430 FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2, INT *result_e); 431 432 /** 433 * \brief Multiply 2 values using maximum precision. The exponent of the result 434 * is 0. 435 * \param f1_m mantissa of factor 1 436 * \param f2_m mantissa of factor 2 437 * \return mantissa of the result with exponent equal to 0 438 */ 439 inline FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2) { 440 FIXP_DBL m; 441 INT e; 442 443 m = fMultNorm(f1, f2, &e); 444 445 m = scaleValueSaturate(m, e); 446 447 return m; 448 } 449 450 /** 451 * \brief Multiply 2 values with exponent and use given exponent for the 452 * mantissa of the result. 453 * \param f1_m mantissa of factor 1 454 * \param f1_e exponent of factor 1 455 * \param f2_m mantissa of factor 2 456 * \param f2_e exponent of factor 2 457 * \param result_e exponent for the returned mantissa of the result 458 * \return mantissa of the result with exponent equal to result_e 459 */ 460 inline FIXP_DBL fMultNorm(FIXP_DBL f1_m, INT f1_e, FIXP_DBL f2_m, INT f2_e, 461 INT result_e) { 462 FIXP_DBL m; 463 INT e; 464 465 m = fMultNorm(f1_m, f2_m, &e); 466 467 m = scaleValueSaturate(m, e + f1_e + f2_e - result_e); 468 469 return m; 470 } 471 #endif /* FUNCTION_fMultNorm */ 472 473 #ifndef FUNCTION_fMultI 474 /** 475 * \brief Multiplies a fractional value and a integer value and performs 476 * rounding to nearest 477 * \param a fractional value 478 * \param b integer value 479 * \return integer value 480 */ 481 inline INT fMultI(FIXP_DBL a, INT b) { 482 FIXP_DBL m, mi; 483 INT m_e; 484 485 m = fMultNorm(a, (FIXP_DBL)b, &m_e); 486 487 if (m_e < (INT)0) { 488 if (m_e > (INT)-DFRACT_BITS) { 489 m = m >> ((-m_e) - 1); 490 mi = (m + (FIXP_DBL)1) >> 1; 491 } else { 492 mi = (FIXP_DBL)0; 493 } 494 } else { 495 mi = scaleValueSaturate(m, m_e); 496 } 497 498 return ((INT)mi); 499 } 500 #endif /* FUNCTION_fMultI */ 501 502 #ifndef FUNCTION_fMultIfloor 503 /** 504 * \brief Multiplies a fractional value and a integer value and performs floor 505 * rounding 506 * \param a fractional value 507 * \param b integer value 508 * \return integer value 509 */ 510 inline INT fMultIfloor(FIXP_DBL a, INT b) { 511 FIXP_DBL m, mi; 512 INT m_e; 513 514 m = fMultNorm(a, (FIXP_DBL)b, &m_e); 515 516 if (m_e < (INT)0) { 517 if (m_e > (INT)-DFRACT_BITS) { 518 mi = m >> (-m_e); 519 } else { 520 mi = (FIXP_DBL)0; 521 if (m < (FIXP_DBL)0) { 522 mi = (FIXP_DBL)-1; 523 } 524 } 525 } else { 526 mi = scaleValueSaturate(m, m_e); 527 } 528 529 return ((INT)mi); 530 } 531 #endif /* FUNCTION_fMultIfloor */ 532 533 #ifndef FUNCTION_fMultIceil 534 /** 535 * \brief Multiplies a fractional value and a integer value and performs ceil 536 * rounding 537 * \param a fractional value 538 * \param b integer value 539 * \return integer value 540 */ 541 inline INT fMultIceil(FIXP_DBL a, INT b) { 542 FIXP_DBL m, mi; 543 INT m_e; 544 545 m = fMultNorm(a, (FIXP_DBL)b, &m_e); 546 547 if (m_e < (INT)0) { 548 if (m_e > (INT)-DFRACT_BITS) { 549 mi = (m >> (-m_e)); 550 if ((LONG)m & ((1 << (-m_e)) - 1)) { 551 mi = mi + (FIXP_DBL)1; 552 } 553 } else { 554 mi = (FIXP_DBL)1; 555 if (m < (FIXP_DBL)0) { 556 mi = (FIXP_DBL)0; 557 } 558 } 559 } else { 560 mi = scaleValueSaturate(m, m_e); 561 } 562 563 return ((INT)mi); 564 } 565 #endif /* FUNCTION_fMultIceil */ 566 567 #ifndef FUNCTION_fDivNorm 568 /** 569 * \brief Divide 2 FIXP_DBL values with normalization of input values. 570 * \param num numerator 571 * \param denum denominator 572 * \param result_e pointer to an INT where the exponent of the result is stored 573 * into 574 * \return num/denum with exponent = *result_e 575 */ 576 FIXP_DBL fDivNorm(FIXP_DBL num, FIXP_DBL denom, INT *result_e); 577 578 /** 579 * \brief Divide 2 positive FIXP_DBL values with normalization of input values. 580 * \param num numerator 581 * \param denum denominator 582 * \return num/denum with exponent = 0 583 */ 584 FIXP_DBL fDivNorm(FIXP_DBL num, FIXP_DBL denom); 585 586 /** 587 * \brief Divide 2 signed FIXP_DBL values with normalization of input values. 588 * \param num numerator 589 * \param denum denominator 590 * \param result_e pointer to an INT where the exponent of the result is stored 591 * into 592 * \return num/denum with exponent = *result_e 593 */ 594 FIXP_DBL fDivNormSigned(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e); 595 596 /** 597 * \brief Divide 2 signed FIXP_DBL values with normalization of input values. 598 * \param num numerator 599 * \param denum denominator 600 * \return num/denum with exponent = 0 601 */ 602 FIXP_DBL fDivNormSigned(FIXP_DBL num, FIXP_DBL denom); 603 #endif /* FUNCTION_fDivNorm */ 604 605 /** 606 * \brief Adjust mantissa to exponent -1 607 * \param a_m mantissa of value to be adjusted 608 * \param pA_e pointer to the exponen of a_m 609 * \return adjusted mantissa 610 */ 611 inline FIXP_DBL fAdjust(FIXP_DBL a_m, INT *pA_e) { 612 INT shift; 613 614 shift = fNorm(a_m) - 1; 615 *pA_e -= shift; 616 617 return scaleValue(a_m, shift); 618 } 619 620 #ifndef FUNCTION_fAddNorm 621 /** 622 * \brief Add two values with normalization 623 * \param a_m mantissa of first summand 624 * \param a_e exponent of first summand 625 * \param a_m mantissa of second summand 626 * \param a_e exponent of second summand 627 * \param pResult_e pointer to where the exponent of the result will be stored 628 * to. 629 * \return mantissa of result 630 */ 631 inline FIXP_DBL fAddNorm(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e, 632 INT *pResult_e) { 633 INT result_e; 634 FIXP_DBL result_m; 635 636 /* If one of the summands is zero, return the other. 637 This is necessary for the summation of a very small number to zero */ 638 if (a_m == (FIXP_DBL)0) { 639 *pResult_e = b_e; 640 return b_m; 641 } 642 if (b_m == (FIXP_DBL)0) { 643 *pResult_e = a_e; 644 return a_m; 645 } 646 647 a_m = fAdjust(a_m, &a_e); 648 b_m = fAdjust(b_m, &b_e); 649 650 if (a_e > b_e) { 651 result_m = a_m + (b_m >> fMin(a_e - b_e, DFRACT_BITS - 1)); 652 result_e = a_e; 653 } else { 654 result_m = (a_m >> fMin(b_e - a_e, DFRACT_BITS - 1)) + b_m; 655 result_e = b_e; 656 } 657 658 *pResult_e = result_e; 659 return result_m; 660 } 661 662 inline FIXP_DBL fAddNorm(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e, 663 INT result_e) { 664 FIXP_DBL result_m; 665 666 a_m = scaleValue(a_m, a_e - result_e); 667 b_m = scaleValue(b_m, b_e - result_e); 668 669 result_m = a_m + b_m; 670 671 return result_m; 672 } 673 #endif /* FUNCTION_fAddNorm */ 674 675 /** 676 * \brief Divide 2 FIXP_DBL values with normalization of input values. 677 * \param num numerator 678 * \param denum denomintator 679 * \return num/denum with exponent = 0 680 */ 681 FIXP_DBL fDivNormHighPrec(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e); 682 683 #ifndef FUNCTION_fPow 684 /** 685 * \brief return 2 ^ (exp_m * 2^exp_e) 686 * \param exp_m mantissa of the exponent to 2.0f 687 * \param exp_e exponent of the exponent to 2.0f 688 * \param result_e pointer to a INT where the exponent of the result will be 689 * stored into 690 * \return mantissa of the result 691 */ 692 FIXP_DBL f2Pow(const FIXP_DBL exp_m, const INT exp_e, INT *result_e); 693 694 /** 695 * \brief return 2 ^ (exp_m * 2^exp_e). This version returns only the mantissa 696 * with implicit exponent of zero. 697 * \param exp_m mantissa of the exponent to 2.0f 698 * \param exp_e exponent of the exponent to 2.0f 699 * \return mantissa of the result 700 */ 701 FIXP_DBL f2Pow(const FIXP_DBL exp_m, const INT exp_e); 702 703 /** 704 * \brief return x ^ (exp_m * 2^exp_e), where log2(x) = baseLd_m * 2^(baseLd_e). 705 * This saves the need to compute log2() of constant values (when x is a 706 * constant). 707 * \param baseLd_m mantissa of log2() of x. 708 * \param baseLd_e exponent of log2() of x. 709 * \param exp_m mantissa of the exponent to 2.0f 710 * \param exp_e exponent of the exponent to 2.0f 711 * \param result_e pointer to a INT where the exponent of the result will be 712 * stored into 713 * \return mantissa of the result 714 */ 715 FIXP_DBL fLdPow(FIXP_DBL baseLd_m, INT baseLd_e, FIXP_DBL exp_m, INT exp_e, 716 INT *result_e); 717 718 /** 719 * \brief return x ^ (exp_m * 2^exp_e), where log2(x) = baseLd_m * 2^(baseLd_e). 720 * This saves the need to compute log2() of constant values (when x is a 721 * constant). This version does not return an exponent, which is 722 * implicitly 0. 723 * \param baseLd_m mantissa of log2() of x. 724 * \param baseLd_e exponent of log2() of x. 725 * \param exp_m mantissa of the exponent to 2.0f 726 * \param exp_e exponent of the exponent to 2.0f 727 * \return mantissa of the result 728 */ 729 FIXP_DBL fLdPow(FIXP_DBL baseLd_m, INT baseLd_e, FIXP_DBL exp_m, INT exp_e); 730 731 /** 732 * \brief return (base_m * 2^base_e) ^ (exp * 2^exp_e). Use fLdPow() instead 733 * whenever possible. 734 * \param base_m mantissa of the base. 735 * \param base_e exponent of the base. 736 * \param exp_m mantissa of power to be calculated of the base. 737 * \param exp_e exponent of power to be calculated of the base. 738 * \param result_e pointer to a INT where the exponent of the result will be 739 * stored into. 740 * \return mantissa of the result. 741 */ 742 FIXP_DBL fPow(FIXP_DBL base_m, INT base_e, FIXP_DBL exp_m, INT exp_e, 743 INT *result_e); 744 745 /** 746 * \brief return (base_m * 2^base_e) ^ N 747 * \param base_m mantissa of the base 748 * \param base_e exponent of the base 749 * \param N power to be calculated of the base 750 * \param result_e pointer to a INT where the exponent of the result will be 751 * stored into 752 * \return mantissa of the result 753 */ 754 FIXP_DBL fPowInt(FIXP_DBL base_m, INT base_e, INT N, INT *result_e); 755 #endif /* #ifndef FUNCTION_fPow */ 756 757 #ifndef FUNCTION_fLog2 758 /** 759 * \brief Calculate log(argument)/log(2) (logarithm with base 2). deprecated. 760 * Use fLog2() instead. 761 * \param arg mantissa of the argument 762 * \param arg_e exponent of the argument 763 * \param result_e pointer to an INT to store the exponent of the result 764 * \return the mantissa of the result. 765 * \param 766 */ 767 FIXP_DBL CalcLog2(FIXP_DBL arg, INT arg_e, INT *result_e); 768 769 /** 770 * \brief calculate logarithm of base 2 of x_m * 2^(x_e) 771 * \param x_m mantissa of the input value. 772 * \param x_e exponent of the input value. 773 * \param pointer to an INT where the exponent of the result is returned into. 774 * \return mantissa of the result. 775 */ 776 FDK_INLINE FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e, INT *result_e) { 777 FIXP_DBL result_m; 778 779 /* Short cut for zero and negative numbers. */ 780 if (x_m <= FL2FXCONST_DBL(0.0f)) { 781 *result_e = DFRACT_BITS - 1; 782 return FL2FXCONST_DBL(-1.0f); 783 } 784 785 /* Calculate log2() */ 786 { 787 FIXP_DBL x2_m; 788 789 /* Move input value x_m * 2^x_e toward 1.0, where the taylor approximation 790 of the function log(1-x) centered at 0 is most accurate. */ 791 { 792 INT b_norm; 793 794 b_norm = fNormz(x_m) - 1; 795 x2_m = x_m << b_norm; 796 x_e = x_e - b_norm; 797 } 798 799 /* map x from log(x) domain to log(1-x) domain. */ 800 x2_m = -(x2_m + FL2FXCONST_DBL(-1.0)); 801 802 /* Taylor polynomial approximation of ln(1-x) */ 803 { 804 FIXP_DBL px2_m; 805 result_m = FL2FXCONST_DBL(0.0); 806 px2_m = x2_m; 807 for (int i = 0; i < LD_PRECISION; i++) { 808 result_m = fMultAddDiv2(result_m, ldCoeff[i], px2_m); 809 px2_m = fMult(px2_m, x2_m); 810 } 811 } 812 /* Multiply result with 1/ln(2) = 1.0 + 0.442695040888 (get log2(x) from 813 * ln(x) result). */ 814 result_m = 815 fMultAddDiv2(result_m, result_m, 816 FL2FXCONST_DBL(2.0 * 0.4426950408889634073599246810019)); 817 818 /* Add exponent part. log2(x_m * 2^x_e) = log2(x_m) + x_e */ 819 if (x_e != 0) { 820 int enorm; 821 822 enorm = DFRACT_BITS - fNorm((FIXP_DBL)x_e); 823 /* The -1 in the right shift of result_m compensates the fMultDiv2() above 824 * in the taylor polynomial evaluation loop.*/ 825 result_m = (result_m >> (enorm - 1)) + 826 ((FIXP_DBL)x_e << (DFRACT_BITS - 1 - enorm)); 827 828 *result_e = enorm; 829 } else { 830 /* 1 compensates the fMultDiv2() above in the taylor polynomial evaluation 831 * loop.*/ 832 *result_e = 1; 833 } 834 } 835 836 return result_m; 837 } 838 839 /** 840 * \brief calculate logarithm of base 2 of x_m * 2^(x_e) 841 * \param x_m mantissa of the input value. 842 * \param x_e exponent of the input value. 843 * \return mantissa of the result with implicit exponent of LD_DATA_SHIFT. 844 */ 845 FDK_INLINE FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e) { 846 if (x_m <= FL2FXCONST_DBL(0.0f)) { 847 x_m = FL2FXCONST_DBL(-1.0f); 848 } else { 849 INT result_e; 850 x_m = fLog2(x_m, x_e, &result_e); 851 x_m = scaleValue(x_m, result_e - LD_DATA_SHIFT); 852 } 853 return x_m; 854 } 855 856 #endif /* FUNCTION_fLog2 */ 857 858 #ifndef FUNCTION_fAddSaturate 859 /** 860 * \brief Add with saturation of the result. 861 * \param a first summand 862 * \param b second summand 863 * \return saturated sum of a and b. 864 */ 865 inline FIXP_SGL fAddSaturate(const FIXP_SGL a, const FIXP_SGL b) { 866 LONG sum; 867 868 sum = (LONG)(SHORT)a + (LONG)(SHORT)b; 869 sum = fMax(fMin((INT)sum, (INT)MAXVAL_SGL), (INT)MINVAL_SGL); 870 return (FIXP_SGL)(SHORT)sum; 871 } 872 873 /** 874 * \brief Add with saturation of the result. 875 * \param a first summand 876 * \param b second summand 877 * \return saturated sum of a and b. 878 */ 879 inline FIXP_DBL fAddSaturate(const FIXP_DBL a, const FIXP_DBL b) { 880 LONG sum; 881 882 sum = (LONG)(a >> 1) + (LONG)(b >> 1); 883 sum = fMax(fMin((INT)sum, (INT)(MAXVAL_DBL >> 1)), (INT)(MINVAL_DBL >> 1)); 884 return (FIXP_DBL)(LONG)(sum << 1); 885 } 886 #endif /* FUNCTION_fAddSaturate */ 887 888 INT fixp_floorToInt(FIXP_DBL f_inp, INT sf); 889 FIXP_DBL fixp_floor(FIXP_DBL f_inp, INT sf); 890 891 INT fixp_ceilToInt(FIXP_DBL f_inp, INT sf); 892 FIXP_DBL fixp_ceil(FIXP_DBL f_inp, INT sf); 893 894 INT fixp_truncateToInt(FIXP_DBL f_inp, INT sf); 895 FIXP_DBL fixp_truncate(FIXP_DBL f_inp, INT sf); 896 897 INT fixp_roundToInt(FIXP_DBL f_inp, INT sf); 898 FIXP_DBL fixp_round(FIXP_DBL f_inp, INT sf); 899 900 /***************************************************************************** 901 902 array for 1/n, n=1..80 903 904 ****************************************************************************/ 905 906 extern const FIXP_DBL invCount[80]; 907 908 LNK_SECTION_INITCODE 909 inline void InitInvInt(void) {} 910 911 /** 912 * \brief Calculate the value of 1/i where i is a integer value. It supports 913 * input values from 1 upto (80-1). 914 * \param intValue Integer input value. 915 * \param FIXP_DBL representation of 1/intValue 916 */ 917 inline FIXP_DBL GetInvInt(int intValue) { 918 return invCount[fMin(fMax(intValue, 0), 80 - 1)]; 919 } 920 921 #endif /* FIXPOINT_MATH_H */ 922
