1 /* -----------------------------------------------------------------------------
2 Software License for The Fraunhofer FDK AAC Codec Library for Android
3
4 © Copyright 1995 - 2018 Fraunhofer-Gesellschaft zur Förderung 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@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 */
fIsLessThan(FIXP_DBL a_m,INT a_e,FIXP_DBL b_m,INT b_e)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
fIsLessThan(FIXP_SGL a_m,INT a_e,FIXP_SGL b_m,INT b_e)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
CalcInvLdData(const FIXP_DBL x)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
sqrtFixp_lookup(FIXP_DBL x)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
sqrtFixp_lookup(FIXP_DBL x,INT * x_e)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 *****************************************************************************/
invSqrtNorm2(FIXP_DBL op,INT * shift)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
sqrtFixp(FIXP_DBL op)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 */
invFixp(FIXP_DBL op)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 */
invFixp(FIXP_DBL op_m,int * op_e)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 */
fMultNorm(FIXP_DBL f1,FIXP_DBL f2)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 */
fMultNorm(FIXP_DBL f1_m,INT f1_e,FIXP_DBL f2_m,INT f2_e,INT result_e)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 */
fMultI(FIXP_DBL a,INT b)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 */
fMultIfloor(FIXP_DBL a,INT b)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 */
fMultIceil(FIXP_DBL a,INT b)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 */
fAdjust(FIXP_DBL a_m,INT * pA_e)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 */
fAddNorm(FIXP_DBL a_m,INT a_e,FIXP_DBL b_m,INT b_e,INT * pResult_e)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
fAddNorm(FIXP_DBL a_m,INT a_e,FIXP_DBL b_m,INT b_e,INT result_e)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 */
fLog2(FIXP_DBL x_m,INT x_e,INT * result_e)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 */
fLog2(FIXP_DBL x_m,INT x_e)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 */
fAddSaturate(const FIXP_SGL a,const FIXP_SGL b)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 */
fAddSaturate(const FIXP_DBL a,const FIXP_DBL b)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
InitInvInt(void)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 */
GetInvInt(int intValue)917 inline FIXP_DBL GetInvInt(int intValue) {
918 return invCount[fMin(fMax(intValue, 0), 80 - 1)];
919 }
920
921 #endif /* FIXPOINT_MATH_H */
922