1
2 /* -----------------------------------------------------------------------------------------------------------
3 Software License for The Fraunhofer FDK AAC Codec Library for Android
4
5 � Copyright 1995 - 2012 Fraunhofer-Gesellschaft zur F�rderung der angewandten Forschung e.V.
6 All rights reserved.
7
8 1. INTRODUCTION
9 The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software that implements
10 the MPEG Advanced Audio Coding ("AAC") encoding and decoding scheme for digital audio.
11 This FDK AAC Codec software is intended to be used on a wide variety of Android devices.
12
13 AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient general perceptual
14 audio codecs. AAC-ELD is considered the best-performing full-bandwidth communications codec by
15 independent studies and is widely deployed. AAC has been standardized by ISO and IEC as part
16 of the MPEG specifications.
17
18 Patent licenses for necessary patent claims for the FDK AAC Codec (including those of Fraunhofer)
19 may be obtained through Via Licensing (www.vialicensing.com) or through the respective patent owners
20 individually for the purpose of encoding or decoding bit streams in products that are compliant with
21 the ISO/IEC MPEG audio standards. Please note that most manufacturers of Android devices already license
22 these patent claims through Via Licensing or directly from the patent owners, and therefore FDK AAC Codec
23 software may already be covered under those patent licenses when it is used for those licensed purposes only.
24
25 Commercially-licensed AAC software libraries, including floating-point versions with enhanced sound quality,
26 are also available from Fraunhofer. Users are encouraged to check the Fraunhofer website for additional
27 applications information and documentation.
28
29 2. COPYRIGHT LICENSE
30
31 Redistribution and use in source and binary forms, with or without modification, are permitted without
32 payment of copyright license fees provided that you satisfy the following conditions:
33
34 You must retain the complete text of this software license in redistributions of the FDK AAC Codec or
35 your modifications thereto in source code form.
36
37 You must retain the complete text of this software license in the documentation and/or other materials
38 provided with redistributions of the FDK AAC Codec or your modifications thereto in binary form.
39 You must make available free of charge copies of the complete source code of the FDK AAC Codec and your
40 modifications thereto to recipients of copies in binary form.
41
42 The name of Fraunhofer may not be used to endorse or promote products derived from this library without
43 prior written permission.
44
45 You may not charge copyright license fees for anyone to use, copy or distribute the FDK AAC Codec
46 software or your modifications thereto.
47
48 Your modified versions of the FDK AAC Codec must carry prominent notices stating that you changed the software
49 and the date of any change. For modified versions of the FDK AAC Codec, the term
50 "Fraunhofer FDK AAC Codec Library for Android" must be replaced by the term
51 "Third-Party Modified Version of the Fraunhofer FDK AAC Codec Library for Android."
52
53 3. NO PATENT LICENSE
54
55 NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without limitation the patents of Fraunhofer,
56 ARE GRANTED BY THIS SOFTWARE LICENSE. Fraunhofer provides no warranty of patent non-infringement with
57 respect to this software.
58
59 You may use this FDK AAC Codec software or modifications thereto only for purposes that are authorized
60 by appropriate patent licenses.
61
62 4. DISCLAIMER
63
64 This FDK AAC Codec software is provided by Fraunhofer on behalf of the copyright holders and contributors
65 "AS IS" and WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES, including but not limited to the implied warranties
66 of merchantability and fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
67 CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary, or consequential damages,
68 including but not limited to procurement of substitute goods or services; loss of use, data, or profits,
69 or business interruption, however caused and on any theory of liability, whether in contract, strict
70 liability, or tort (including negligence), arising in any way out of the use of this software, even if
71 advised of the possibility of such damage.
72
73 5. CONTACT INFORMATION
74
75 Fraunhofer Institute for Integrated Circuits IIS
76 Attention: Audio and Multimedia Departments - FDK AAC LL
77 Am Wolfsmantel 33
78 91058 Erlangen, Germany
79
80 www.iis.fraunhofer.de/amm
81 amm-info@iis.fraunhofer.de
82 ----------------------------------------------------------------------------------------------------------- */
83
84 /*************************** Fraunhofer IIS FDK Tools **********************
85
86 Author(s): M. Gayer
87 Description: Fixed point specific mathematical functions
88
89 ******************************************************************************/
90
91 #include "fixpoint_math.h"
92
93
94 #define MAX_LD_PRECISION 10
95 #define LD_PRECISION 10
96
97 /* Taylor series coeffcients for ln(1-x), centered at 0 (MacLaurin polinomial). */
98 #ifndef LDCOEFF_16BIT
99 LNK_SECTION_CONSTDATA_L1
100 static const FIXP_DBL ldCoeff[MAX_LD_PRECISION] = {
101 FL2FXCONST_DBL(-1.0),
102 FL2FXCONST_DBL(-1.0/2.0),
103 FL2FXCONST_DBL(-1.0/3.0),
104 FL2FXCONST_DBL(-1.0/4.0),
105 FL2FXCONST_DBL(-1.0/5.0),
106 FL2FXCONST_DBL(-1.0/6.0),
107 FL2FXCONST_DBL(-1.0/7.0),
108 FL2FXCONST_DBL(-1.0/8.0),
109 FL2FXCONST_DBL(-1.0/9.0),
110 FL2FXCONST_DBL(-1.0/10.0)
111 };
112 #else
113 LNK_SECTION_CONSTDATA_L1
114 static const FIXP_SGL ldCoeff[MAX_LD_PRECISION] = {
115 FL2FXCONST_SGL(-1.0),
116 FL2FXCONST_SGL(-1.0/2.0),
117 FL2FXCONST_SGL(-1.0/3.0),
118 FL2FXCONST_SGL(-1.0/4.0),
119 FL2FXCONST_SGL(-1.0/5.0),
120 FL2FXCONST_SGL(-1.0/6.0),
121 FL2FXCONST_SGL(-1.0/7.0),
122 FL2FXCONST_SGL(-1.0/8.0),
123 FL2FXCONST_SGL(-1.0/9.0),
124 FL2FXCONST_SGL(-1.0/10.0)
125 };
126 #endif
127
128 /*****************************************************************************
129
130 functionname: CalcLdData
131 description: Delivers the Logarithm Dualis ld(op)/LD_DATA_SCALING with polynomial approximation.
132 input: Input op is assumed to be double precision fractional 0 < op < 1.0
133 This function does not accept negative values.
134 output: For op == 0, the result is saturated to -1.0
135 This function does not return positive values since input values are treated as fractional values.
136 It does not make sense to input an integer value into this function (and expect a positive output value)
137 since input values are treated as fractional values.
138
139 *****************************************************************************/
140
141 LNK_SECTION_CODE_L1
CalcLdData(FIXP_DBL op)142 FIXP_DBL CalcLdData(FIXP_DBL op)
143 {
144 return fLog2(op, 0);
145 }
146
147
148 /*****************************************************************************
149 functionname: LdDataVector
150 *****************************************************************************/
151 LNK_SECTION_CODE_L1
LdDataVector(FIXP_DBL * srcVector,FIXP_DBL * destVector,INT n)152 void LdDataVector( FIXP_DBL *srcVector,
153 FIXP_DBL *destVector,
154 INT n)
155 {
156 INT i;
157 for ( i=0; i<n; i++) {
158 destVector[i] = CalcLdData(srcVector[i]);
159 }
160 }
161
162
163
164 #define MAX_POW2_PRECISION 8
165 #ifndef SINETABLE_16BIT
166 #define POW2_PRECISION MAX_POW2_PRECISION
167 #else
168 #define POW2_PRECISION 5
169 #endif
170
171 /*
172 Taylor series coefficients of the function x^2. The first coefficient is
173 ommited (equal to 1.0).
174
175 pow2Coeff[i-1] = (1/i!) d^i(2^x)/dx^i, i=1..MAX_POW2_PRECISION
176 To evaluate the taylor series around x = 0, the coefficients are: 1/!i * ln(2)^i
177 */
178 #ifndef POW2COEFF_16BIT
179 LNK_SECTION_CONSTDATA_L1
180 static const FIXP_DBL pow2Coeff[MAX_POW2_PRECISION] = {
181 FL2FXCONST_DBL(0.693147180559945309417232121458177), /* ln(2)^1 /1! */
182 FL2FXCONST_DBL(0.240226506959100712333551263163332), /* ln(2)^2 /2! */
183 FL2FXCONST_DBL(0.0555041086648215799531422637686218), /* ln(2)^3 /3! */
184 FL2FXCONST_DBL(0.00961812910762847716197907157365887), /* ln(2)^4 /4! */
185 FL2FXCONST_DBL(0.00133335581464284434234122219879962), /* ln(2)^5 /5! */
186 FL2FXCONST_DBL(1.54035303933816099544370973327423e-4), /* ln(2)^6 /6! */
187 FL2FXCONST_DBL(1.52527338040598402800254390120096e-5), /* ln(2)^7 /7! */
188 FL2FXCONST_DBL(1.32154867901443094884037582282884e-6) /* ln(2)^8 /8! */
189 };
190 #else
191 LNK_SECTION_CONSTDATA_L1
192 static const FIXP_SGL pow2Coeff[MAX_POW2_PRECISION] = {
193 FL2FXCONST_SGL(0.693147180559945309417232121458177), /* ln(2)^1 /1! */
194 FL2FXCONST_SGL(0.240226506959100712333551263163332), /* ln(2)^2 /2! */
195 FL2FXCONST_SGL(0.0555041086648215799531422637686218), /* ln(2)^3 /3! */
196 FL2FXCONST_SGL(0.00961812910762847716197907157365887), /* ln(2)^4 /4! */
197 FL2FXCONST_SGL(0.00133335581464284434234122219879962), /* ln(2)^5 /5! */
198 FL2FXCONST_SGL(1.54035303933816099544370973327423e-4), /* ln(2)^6 /6! */
199 FL2FXCONST_SGL(1.52527338040598402800254390120096e-5), /* ln(2)^7 /7! */
200 FL2FXCONST_SGL(1.32154867901443094884037582282884e-6) /* ln(2)^8 /8! */
201 };
202 #endif
203
204
205
206 /*****************************************************************************
207
208 functionname: mul_dbl_sgl_rnd
209 description: Multiply with round.
210 *****************************************************************************/
211
212 /* for rounding a dfract to fract */
213 #define ACCU_R (LONG) 0x00008000
214
215 LNK_SECTION_CODE_L1
mul_dbl_sgl_rnd(const FIXP_DBL op1,const FIXP_SGL op2)216 FIXP_DBL mul_dbl_sgl_rnd (const FIXP_DBL op1, const FIXP_SGL op2)
217 {
218 FIXP_DBL prod;
219 LONG v = (LONG)(op1);
220 SHORT u = (SHORT)(op2);
221
222 LONG low = u*(v&SGL_MASK);
223 low = (low+(ACCU_R>>1)) >> (FRACT_BITS-1); /* round */
224 LONG high = u * ((v>>FRACT_BITS)<<1);
225
226 prod = (LONG)(high+low);
227
228 return((FIXP_DBL)prod);
229 }
230
231
232 /*****************************************************************************
233
234 functionname: CalcInvLdData
235 description: Delivers the inverse of function CalcLdData().
236 Delivers 2^(op*LD_DATA_SCALING)
237 input: Input op is assumed to be fractional -1.0 < op < 1.0
238 output: For op == 0, the result is MAXVAL_DBL (almost 1.0).
239 For negative input values the output should be treated as a positive fractional value.
240 For positive input values the output should be treated as a positive integer value.
241 This function does not output negative values.
242
243 *****************************************************************************/
244 LNK_SECTION_CODE_L1
CalcInvLdData(FIXP_DBL op)245 FIXP_DBL CalcInvLdData(FIXP_DBL op)
246 {
247 FIXP_DBL result_m;
248
249 if ( op == FL2FXCONST_DBL(0.0f) ) {
250 result_m = (FIXP_DBL)MAXVAL_DBL;
251 }
252 else if ( op < FL2FXCONST_DBL(0.0f) ) {
253 result_m = f2Pow(op, LD_DATA_SHIFT);
254 }
255 else {
256 int result_e;
257
258 result_m = f2Pow(op, LD_DATA_SHIFT, &result_e);
259 result_e = fixMin(fixMax(result_e+1-(DFRACT_BITS-1), -(DFRACT_BITS-1)), (DFRACT_BITS-1)); /* rounding and saturation */
260
261 if ( (result_e>0) && ( result_m > (((FIXP_DBL)MAXVAL_DBL)>>result_e) ) ) {
262 result_m = (FIXP_DBL)MAXVAL_DBL; /* saturate to max representable value */
263 }
264 else {
265 result_m = (scaleValue(result_m, result_e)+(FIXP_DBL)1)>>1; /* descale result + rounding */
266 }
267 }
268 return result_m;
269 }
270
271
272
273
274
275 /*****************************************************************************
276 functionname: InitLdInt and CalcLdInt
277 description: Create and access table with integer LdData (0 to 193)
278 *****************************************************************************/
279
280
281 LNK_SECTION_CONSTDATA_L1
282 static const FIXP_DBL ldIntCoeff[] = {
283 0x80000001, 0x00000000, 0x02000000, 0x032b8034, 0x04000000, 0x04a4d3c2, 0x052b8034, 0x059d5da0,
284 0x06000000, 0x06570069, 0x06a4d3c2, 0x06eb3a9f, 0x072b8034, 0x0766a009, 0x079d5da0, 0x07d053f7,
285 0x08000000, 0x082cc7ee, 0x08570069, 0x087ef05b, 0x08a4d3c2, 0x08c8ddd4, 0x08eb3a9f, 0x090c1050,
286 0x092b8034, 0x0949a785, 0x0966a009, 0x0982809d, 0x099d5da0, 0x09b74949, 0x09d053f7, 0x09e88c6b,
287 0x0a000000, 0x0a16bad3, 0x0a2cc7ee, 0x0a423162, 0x0a570069, 0x0a6b3d79, 0x0a7ef05b, 0x0a92203d,
288 0x0aa4d3c2, 0x0ab7110e, 0x0ac8ddd4, 0x0ada3f60, 0x0aeb3a9f, 0x0afbd42b, 0x0b0c1050, 0x0b1bf312,
289 0x0b2b8034, 0x0b3abb40, 0x0b49a785, 0x0b584822, 0x0b66a009, 0x0b74b1fd, 0x0b82809d, 0x0b900e61,
290 0x0b9d5da0, 0x0baa708f, 0x0bb74949, 0x0bc3e9ca, 0x0bd053f7, 0x0bdc899b, 0x0be88c6b, 0x0bf45e09,
291 0x0c000000, 0x0c0b73cb, 0x0c16bad3, 0x0c21d671, 0x0c2cc7ee, 0x0c379085, 0x0c423162, 0x0c4caba8,
292 0x0c570069, 0x0c6130af, 0x0c6b3d79, 0x0c7527b9, 0x0c7ef05b, 0x0c88983f, 0x0c92203d, 0x0c9b8926,
293 0x0ca4d3c2, 0x0cae00d2, 0x0cb7110e, 0x0cc0052b, 0x0cc8ddd4, 0x0cd19bb0, 0x0cda3f60, 0x0ce2c97d,
294 0x0ceb3a9f, 0x0cf39355, 0x0cfbd42b, 0x0d03fda9, 0x0d0c1050, 0x0d140ca0, 0x0d1bf312, 0x0d23c41d,
295 0x0d2b8034, 0x0d3327c7, 0x0d3abb40, 0x0d423b08, 0x0d49a785, 0x0d510118, 0x0d584822, 0x0d5f7cff,
296 0x0d66a009, 0x0d6db197, 0x0d74b1fd, 0x0d7ba190, 0x0d82809d, 0x0d894f75, 0x0d900e61, 0x0d96bdad,
297 0x0d9d5da0, 0x0da3ee7f, 0x0daa708f, 0x0db0e412, 0x0db74949, 0x0dbda072, 0x0dc3e9ca, 0x0dca258e,
298 0x0dd053f7, 0x0dd6753e, 0x0ddc899b, 0x0de29143, 0x0de88c6b, 0x0dee7b47, 0x0df45e09, 0x0dfa34e1,
299 0x0e000000, 0x0e05bf94, 0x0e0b73cb, 0x0e111cd2, 0x0e16bad3, 0x0e1c4dfb, 0x0e21d671, 0x0e275460,
300 0x0e2cc7ee, 0x0e323143, 0x0e379085, 0x0e3ce5d8, 0x0e423162, 0x0e477346, 0x0e4caba8, 0x0e51daa8,
301 0x0e570069, 0x0e5c1d0b, 0x0e6130af, 0x0e663b74, 0x0e6b3d79, 0x0e7036db, 0x0e7527b9, 0x0e7a1030,
302 0x0e7ef05b, 0x0e83c857, 0x0e88983f, 0x0e8d602e, 0x0e92203d, 0x0e96d888, 0x0e9b8926, 0x0ea03232,
303 0x0ea4d3c2, 0x0ea96df0, 0x0eae00d2, 0x0eb28c7f, 0x0eb7110e, 0x0ebb8e96, 0x0ec0052b, 0x0ec474e4,
304 0x0ec8ddd4, 0x0ecd4012, 0x0ed19bb0, 0x0ed5f0c4, 0x0eda3f60, 0x0ede8797, 0x0ee2c97d, 0x0ee70525,
305 0x0eeb3a9f, 0x0eef69ff, 0x0ef39355, 0x0ef7b6b4, 0x0efbd42b, 0x0effebcd, 0x0f03fda9, 0x0f0809cf,
306 0x0f0c1050, 0x0f10113b, 0x0f140ca0, 0x0f18028d, 0x0f1bf312, 0x0f1fde3d, 0x0f23c41d, 0x0f27a4c0,
307 0x0f2b8034
308 };
309
310
311 LNK_SECTION_INITCODE
InitLdInt()312 void InitLdInt()
313 {
314 /* nothing to do! Use preinitialized logarithm table */
315 }
316
317
318
319 LNK_SECTION_CODE_L1
CalcLdInt(INT i)320 FIXP_DBL CalcLdInt(INT i)
321 {
322 /* calculates ld(op)/LD_DATA_SCALING */
323 /* op is assumed to be an integer value between 1 and 193 */
324
325 FDK_ASSERT((193>0) && ((FIXP_DBL)ldIntCoeff[0]==(FIXP_DBL)0x80000001)); /* tab has to be initialized */
326
327 if ((i>0)&&(i<193))
328 return ldIntCoeff[i];
329 else
330 {
331 return (0);
332 }
333 }
334
335
336 /*****************************************************************************
337
338 functionname: invSqrtNorm2
339 description: delivers 1/sqrt(op) normalized to .5...1 and the shift value of the OUTPUT
340
341 *****************************************************************************/
342 #define SQRT_BITS 7
343 #define SQRT_VALUES 128
344 #define SQRT_BITS_MASK 0x7f
345
346 LNK_SECTION_CONSTDATA_L1
347 static const FIXP_DBL invSqrtTab[SQRT_VALUES] = {
348 0x5a827999, 0x5a287e03, 0x59cf8cbb, 0x5977a0ab, 0x5920b4de, 0x58cac480, 0x5875cade, 0x5821c364,
349 0x57cea99c, 0x577c792f, 0x572b2ddf, 0x56dac38d, 0x568b3631, 0x563c81df, 0x55eea2c3, 0x55a19521,
350 0x55555555, 0x5509dfd0, 0x54bf311a, 0x547545d0, 0x542c1aa3, 0x53e3ac5a, 0x539bf7cc, 0x5354f9e6,
351 0x530eafa4, 0x52c91617, 0x52842a5e, 0x523fe9ab, 0x51fc513f, 0x51b95e6b, 0x51770e8e, 0x51355f19,
352 0x50f44d88, 0x50b3d768, 0x5073fa4f, 0x5034b3e6, 0x4ff601df, 0x4fb7e1f9, 0x4f7a5201, 0x4f3d4fce,
353 0x4f00d943, 0x4ec4ec4e, 0x4e8986e9, 0x4e4ea718, 0x4e144ae8, 0x4dda7072, 0x4da115d9, 0x4d683948,
354 0x4d2fd8f4, 0x4cf7f31b, 0x4cc08604, 0x4c898fff, 0x4c530f64, 0x4c1d0293, 0x4be767f5, 0x4bb23df9,
355 0x4b7d8317, 0x4b4935ce, 0x4b1554a6, 0x4ae1de2a, 0x4aaed0f0, 0x4a7c2b92, 0x4a49ecb3, 0x4a1812fa,
356 0x49e69d16, 0x49b589bb, 0x4984d7a4, 0x49548591, 0x49249249, 0x48f4fc96, 0x48c5c34a, 0x4896e53c,
357 0x48686147, 0x483a364c, 0x480c6331, 0x47dee6e0, 0x47b1c049, 0x4784ee5f, 0x4758701c, 0x472c447c,
358 0x47006a80, 0x46d4e130, 0x46a9a793, 0x467ebcb9, 0x46541fb3, 0x4629cf98, 0x45ffcb80, 0x45d61289,
359 0x45aca3d5, 0x45837e88, 0x455aa1ca, 0x45320cc8, 0x4509beb0, 0x44e1b6b4, 0x44b9f40b, 0x449275ec,
360 0x446b3b95, 0x44444444, 0x441d8f3b, 0x43f71bbe, 0x43d0e917, 0x43aaf68e, 0x43854373, 0x435fcf14,
361 0x433a98c5, 0x43159fdb, 0x42f0e3ae, 0x42cc6397, 0x42a81ef5, 0x42841527, 0x4260458d, 0x423caf8c,
362 0x4219528b, 0x41f62df1, 0x41d3412a, 0x41b08ba1, 0x418e0cc7, 0x416bc40d, 0x4149b0e4, 0x4127d2c3,
363 0x41062920, 0x40e4b374, 0x40c3713a, 0x40a261ef, 0x40818511, 0x4060da21, 0x404060a1, 0x40201814
364 };
365
366 LNK_SECTION_INITCODE
InitInvSqrtTab()367 void InitInvSqrtTab()
368 {
369 /* nothing to do !
370 use preinitialized square root table
371 */
372 }
373
374
375
376 #if !defined(FUNCTION_invSqrtNorm2)
377 /*****************************************************************************
378 delivers 1/sqrt(op) normalized to .5...1 and the shift value of the OUTPUT,
379 i.e. the denormalized result is 1/sqrt(op) = invSqrtNorm(op) * 2^(shift)
380 uses Newton-iteration for approximation
381 Q(n+1) = Q(n) + Q(n) * (0.5 - 2 * V * Q(n)^2)
382 with Q = 0.5* V ^-0.5; 0.5 <= V < 1.0
383 *****************************************************************************/
invSqrtNorm2(FIXP_DBL op,INT * shift)384 FIXP_DBL invSqrtNorm2(FIXP_DBL op, INT *shift)
385 {
386
387 FIXP_DBL val = op ;
388 FIXP_DBL reg1, reg2, regtmp ;
389
390 if (val == FL2FXCONST_DBL(0.0)) {
391 *shift = 1 ;
392 return((LONG)1); /* minimum positive value */
393 }
394
395
396 /* normalize input, calculate shift value */
397 FDK_ASSERT(val > FL2FXCONST_DBL(0.0));
398 *shift = fNormz(val) - 1; /* CountLeadingBits() is not necessary here since test value is always > 0 */
399 val <<=*shift ; /* normalized input V */
400 *shift+=2 ; /* bias for exponent */
401
402 /* Newton iteration of 1/sqrt(V) */
403 reg1 = invSqrtTab[ (INT)(val>>(DFRACT_BITS-1-(SQRT_BITS+1))) & SQRT_BITS_MASK ];
404 reg2 = FL2FXCONST_DBL(0.0625f); /* 0.5 >> 3 */
405
406 regtmp= fPow2Div2(reg1); /* a = Q^2 */
407 regtmp= reg2 - fMultDiv2(regtmp, val); /* b = 0.5 - 2 * V * Q^2 */
408 reg1 += (fMultDiv2(regtmp, reg1)<<4); /* Q = Q + Q*b */
409
410 /* calculate the output exponent = input exp/2 */
411 if (*shift & 0x00000001) { /* odd shift values ? */
412 reg2 = FL2FXCONST_DBL(0.707106781186547524400844362104849f); /* 1/sqrt(2); */
413 reg1 = fMultDiv2(reg1, reg2) << 2;
414 }
415
416 *shift = *shift>>1;
417
418 return(reg1);
419 }
420 #endif /* !defined(FUNCTION_invSqrtNorm2) */
421
422 /*****************************************************************************
423
424 functionname: sqrtFixp
425 description: delivers sqrt(op)
426
427 *****************************************************************************/
sqrtFixp(FIXP_DBL op)428 FIXP_DBL sqrtFixp(FIXP_DBL op)
429 {
430 INT tmp_exp = 0;
431 FIXP_DBL tmp_inv = invSqrtNorm2(op, &tmp_exp);
432
433 FDK_ASSERT(tmp_exp > 0) ;
434 return( (FIXP_DBL) ( fMultDiv2( (op<<(tmp_exp-1)), tmp_inv ) << 2 ));
435 }
436
437
438 #if !defined(FUNCTION_schur_div)
439 /*****************************************************************************
440
441 functionname: schur_div
442 description: delivers op1/op2 with op3-bit accuracy
443
444 *****************************************************************************/
445
446
schur_div(FIXP_DBL num,FIXP_DBL denum,INT count)447 FIXP_DBL schur_div(FIXP_DBL num, FIXP_DBL denum, INT count)
448 {
449 INT L_num = (LONG)num>>1;
450 INT L_denum = (LONG)denum>>1;
451 INT div = 0;
452 INT k = count;
453
454 FDK_ASSERT (num>=(FIXP_DBL)0);
455 FDK_ASSERT (denum>(FIXP_DBL)0);
456 FDK_ASSERT (num <= denum);
457
458 if (L_num != 0)
459 while (--k)
460 {
461 div <<= 1;
462 L_num <<= 1;
463 if (L_num >= L_denum)
464 {
465 L_num -= L_denum;
466 div++;
467 }
468 }
469 return (FIXP_DBL)(div << (DFRACT_BITS - count));
470 }
471
472
473 #endif /* !defined(FUNCTION_schur_div) */
474
475
476 #ifndef FUNCTION_fMultNorm
fMultNorm(FIXP_DBL f1,FIXP_DBL f2,INT * result_e)477 FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2, INT *result_e)
478 {
479 INT product = 0;
480 INT norm_f1, norm_f2;
481
482 if ( (f1 == (FIXP_DBL)0) || (f2 == (FIXP_DBL)0) ) {
483 *result_e = 0;
484 return (FIXP_DBL)0;
485 }
486 norm_f1 = CountLeadingBits(f1);
487 f1 = f1 << norm_f1;
488 norm_f2 = CountLeadingBits(f2);
489 f2 = f2 << norm_f2;
490
491 product = fMult(f1, f2);
492 *result_e = - (norm_f1 + norm_f2);
493
494 return (FIXP_DBL)product;
495 }
496 #endif
497
498 #ifndef FUNCTION_fDivNorm
fDivNorm(FIXP_DBL L_num,FIXP_DBL L_denum,INT * result_e)499 FIXP_DBL fDivNorm(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e)
500 {
501 FIXP_DBL div;
502 INT norm_num, norm_den;
503
504 FDK_ASSERT (L_num >= (FIXP_DBL)0);
505 FDK_ASSERT (L_denum > (FIXP_DBL)0);
506
507 if(L_num == (FIXP_DBL)0)
508 {
509 *result_e = 0;
510 return ((FIXP_DBL)0);
511 }
512
513 norm_num = CountLeadingBits(L_num);
514 L_num = L_num << norm_num;
515 L_num = L_num >> 1;
516 *result_e = - norm_num + 1;
517
518 norm_den = CountLeadingBits(L_denum);
519 L_denum = L_denum << norm_den;
520 *result_e -= - norm_den;
521
522 div = schur_div(L_num, L_denum, FRACT_BITS);
523
524 return div;
525 }
526 #endif /* !FUNCTION_fDivNorm */
527
528 #ifndef FUNCTION_fDivNorm
fDivNorm(FIXP_DBL num,FIXP_DBL denom)529 FIXP_DBL fDivNorm(FIXP_DBL num, FIXP_DBL denom)
530 {
531 INT e;
532 FIXP_DBL res;
533
534 FDK_ASSERT (denom >= num);
535
536 res = fDivNorm(num, denom, &e);
537
538 /* Avoid overflow since we must output a value with exponent 0
539 there is no other choice than saturating to almost 1.0f */
540 if(res == (FIXP_DBL)(1<<(DFRACT_BITS-2)) && e == 1)
541 {
542 res = (FIXP_DBL)MAXVAL_DBL;
543 }
544 else
545 {
546 res = scaleValue(res, e);
547 }
548
549 return res;
550 }
551 #endif /* !FUNCTION_fDivNorm */
552
553 #ifndef FUNCTION_fDivNormHighPrec
fDivNormHighPrec(FIXP_DBL num,FIXP_DBL denom,INT * result_e)554 FIXP_DBL fDivNormHighPrec(FIXP_DBL num, FIXP_DBL denom, INT *result_e)
555 {
556 FIXP_DBL div;
557 INT norm_num, norm_den;
558
559 FDK_ASSERT (num >= (FIXP_DBL)0);
560 FDK_ASSERT (denom > (FIXP_DBL)0);
561
562 if(num == (FIXP_DBL)0)
563 {
564 *result_e = 0;
565 return ((FIXP_DBL)0);
566 }
567
568 norm_num = CountLeadingBits(num);
569 num = num << norm_num;
570 num = num >> 1;
571 *result_e = - norm_num + 1;
572
573 norm_den = CountLeadingBits(denom);
574 denom = denom << norm_den;
575 *result_e -= - norm_den;
576
577 div = schur_div(num, denom, 31);
578 return div;
579 }
580 #endif /* !FUNCTION_fDivNormHighPrec */
581
582
583
CalcLog2(FIXP_DBL base_m,INT base_e,INT * result_e)584 FIXP_DBL CalcLog2(FIXP_DBL base_m, INT base_e, INT *result_e)
585 {
586 return fLog2(base_m, base_e, result_e);
587 }
588
f2Pow(const FIXP_DBL exp_m,const INT exp_e,INT * result_e)589 FIXP_DBL f2Pow(
590 const FIXP_DBL exp_m, const INT exp_e,
591 INT *result_e
592 )
593 {
594 FIXP_DBL frac_part, result_m;
595 INT int_part;
596
597 if (exp_e > 0)
598 {
599 INT exp_bits = DFRACT_BITS-1 - exp_e;
600 int_part = exp_m >> exp_bits;
601 frac_part = exp_m - (FIXP_DBL)(int_part << exp_bits);
602 frac_part = frac_part << exp_e;
603 }
604 else
605 {
606 int_part = 0;
607 frac_part = exp_m >> -exp_e;
608 }
609
610 /* Best accuracy is around 0, so try to get there with the fractional part. */
611 if( frac_part > FL2FXCONST_DBL(0.5f) )
612 {
613 int_part = int_part + 1;
614 frac_part = frac_part + FL2FXCONST_DBL(-1.0f);
615 }
616 if( frac_part < FL2FXCONST_DBL(-0.5f) )
617 {
618 int_part = int_part - 1;
619 frac_part = -(FL2FXCONST_DBL(-1.0f) - frac_part);
620 }
621
622 /* Evaluate taylor polynomial which approximates 2^x */
623 {
624 FIXP_DBL p;
625
626 /* result_m ~= 2^frac_part */
627 p = frac_part;
628 /* First taylor series coefficient a_0 = 1.0, scaled by 0.5 due to fMultDiv2(). */
629 result_m = FL2FXCONST_DBL(1.0f/2.0f);
630 for (INT i = 0; i < POW2_PRECISION; i++) {
631 /* next taylor series term: a_i * x^i, x=0 */
632 result_m = fMultAddDiv2(result_m, pow2Coeff[i], p);
633 p = fMult(p, frac_part);
634 }
635 }
636
637 /* "+ 1" compensates fMultAddDiv2() of the polynomial evaluation above. */
638 *result_e = int_part + 1;
639
640 return result_m;
641 }
642
f2Pow(const FIXP_DBL exp_m,const INT exp_e)643 FIXP_DBL f2Pow(
644 const FIXP_DBL exp_m, const INT exp_e
645 )
646 {
647 FIXP_DBL result_m;
648 INT result_e;
649
650 result_m = f2Pow(exp_m, exp_e, &result_e);
651 result_e = fixMin(DFRACT_BITS-1,fixMax(-(DFRACT_BITS-1),result_e));
652
653 return scaleValue(result_m, result_e);
654 }
655
fPow(FIXP_DBL base_m,INT base_e,FIXP_DBL exp_m,INT exp_e,INT * result_e)656 FIXP_DBL fPow(
657 FIXP_DBL base_m, INT base_e,
658 FIXP_DBL exp_m, INT exp_e,
659 INT *result_e
660 )
661 {
662 INT ans_lg2_e, baselg2_e;
663 FIXP_DBL base_lg2, ans_lg2, result;
664
665 /* Calc log2 of base */
666 base_lg2 = fLog2(base_m, base_e, &baselg2_e);
667
668 /* Prepare exp */
669 {
670 INT leadingBits;
671
672 leadingBits = CountLeadingBits(fAbs(exp_m));
673 exp_m = exp_m << leadingBits;
674 exp_e -= leadingBits;
675 }
676
677 /* Calc base pow exp */
678 ans_lg2 = fMult(base_lg2, exp_m);
679 ans_lg2_e = exp_e + baselg2_e;
680
681 /* Calc antilog */
682 result = f2Pow(ans_lg2, ans_lg2_e, result_e);
683
684 return result;
685 }
686
fLdPow(FIXP_DBL baseLd_m,INT baseLd_e,FIXP_DBL exp_m,INT exp_e,INT * result_e)687 FIXP_DBL fLdPow(
688 FIXP_DBL baseLd_m,
689 INT baseLd_e,
690 FIXP_DBL exp_m, INT exp_e,
691 INT *result_e
692 )
693 {
694 INT ans_lg2_e;
695 FIXP_DBL ans_lg2, result;
696
697 /* Prepare exp */
698 {
699 INT leadingBits;
700
701 leadingBits = CountLeadingBits(fAbs(exp_m));
702 exp_m = exp_m << leadingBits;
703 exp_e -= leadingBits;
704 }
705
706 /* Calc base pow exp */
707 ans_lg2 = fMult(baseLd_m, exp_m);
708 ans_lg2_e = exp_e + baseLd_e;
709
710 /* Calc antilog */
711 result = f2Pow(ans_lg2, ans_lg2_e, result_e);
712
713 return result;
714 }
715
fLdPow(FIXP_DBL baseLd_m,INT baseLd_e,FIXP_DBL exp_m,INT exp_e)716 FIXP_DBL fLdPow(
717 FIXP_DBL baseLd_m, INT baseLd_e,
718 FIXP_DBL exp_m, INT exp_e
719 )
720 {
721 FIXP_DBL result_m;
722 int result_e;
723
724 result_m = fLdPow(baseLd_m, baseLd_e, exp_m, exp_e, &result_e);
725
726 return SATURATE_SHIFT(result_m, -result_e, DFRACT_BITS);
727 }
728
fPowInt(FIXP_DBL base_m,INT base_e,INT exp,INT * pResult_e)729 FIXP_DBL fPowInt(
730 FIXP_DBL base_m, INT base_e,
731 INT exp,
732 INT *pResult_e
733 )
734 {
735 FIXP_DBL result;
736
737 if (exp != 0) {
738 INT result_e = 0;
739
740 if (base_m != (FIXP_DBL)0) {
741 {
742 INT leadingBits;
743 leadingBits = CountLeadingBits( base_m );
744 base_m <<= leadingBits;
745 base_e -= leadingBits;
746 }
747
748 result = base_m;
749
750 {
751 int i;
752 for (i = 1; i < fAbs(exp); i++) {
753 result = fMult(result, base_m);
754 }
755 }
756
757 if (exp < 0) {
758 /* 1.0 / ans */
759 result = fDivNorm( FL2FXCONST_DBL(0.5f), result, &result_e );
760 result_e++;
761 } else {
762 int ansScale = CountLeadingBits( result );
763 result <<= ansScale;
764 result_e -= ansScale;
765 }
766
767 result_e += exp * base_e;
768
769 } else {
770 result = (FIXP_DBL)0;
771 }
772 *pResult_e = result_e;
773 }
774 else {
775 result = FL2FXCONST_DBL(0.5f);
776 *pResult_e = 1;
777 }
778
779 return result;
780 }
781
fLog2(FIXP_DBL x_m,INT x_e,INT * result_e)782 FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e, INT *result_e)
783 {
784 FIXP_DBL result_m;
785
786 /* Short cut for zero and negative numbers. */
787 if ( x_m <= FL2FXCONST_DBL(0.0f) ) {
788 *result_e = DFRACT_BITS-1;
789 return FL2FXCONST_DBL(-1.0f);
790 }
791
792 /* Calculate log2() */
793 {
794 FIXP_DBL px2_m, x2_m;
795
796 /* Move input value x_m * 2^x_e toward 1.0, where the taylor approximation
797 of the function log(1-x) centered at 0 is most accurate. */
798 {
799 INT b_norm;
800
801 b_norm = fNormz(x_m)-1;
802 x2_m = x_m << b_norm;
803 x_e = x_e - b_norm;
804 }
805
806 /* map x from log(x) domain to log(1-x) domain. */
807 x2_m = - (x2_m + FL2FXCONST_DBL(-1.0) );
808
809 /* Taylor polinomial approximation of ln(1-x) */
810 result_m = FL2FXCONST_DBL(0.0);
811 px2_m = x2_m;
812 for (int i=0; i<LD_PRECISION; i++) {
813 result_m = fMultAddDiv2(result_m, ldCoeff[i], px2_m);
814 px2_m = fMult(px2_m, x2_m);
815 }
816 /* Multiply result with 1/ln(2) = 1.0 + 0.442695040888 (get log2(x) from ln(x) result). */
817 result_m = fMultAddDiv2(result_m, result_m, FL2FXCONST_DBL(2.0*0.4426950408889634073599246810019));
818
819 /* Add exponent part. log2(x_m * 2^x_e) = log2(x_m) + x_e */
820 if (x_e != 0)
821 {
822 int enorm;
823
824 enorm = DFRACT_BITS - fNorm((FIXP_DBL)x_e);
825 /* The -1 in the right shift of result_m compensates the fMultDiv2() above in the taylor polinomial evaluation loop.*/
826 result_m = (result_m >> (enorm-1)) + ((FIXP_DBL)x_e << (DFRACT_BITS-1-enorm));
827
828 *result_e = enorm;
829 } else {
830 /* 1 compensates the fMultDiv2() above in the taylor polinomial evaluation loop.*/
831 *result_e = 1;
832 }
833 }
834
835 return result_m;
836 }
837
fLog2(FIXP_DBL x_m,INT x_e)838 FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e)
839 {
840 if ( x_m <= FL2FXCONST_DBL(0.0f) ) {
841 x_m = FL2FXCONST_DBL(-1.0f);
842 }
843 else {
844 INT result_e;
845 x_m = fLog2(x_m, x_e, &result_e);
846 x_m = scaleValue(x_m, result_e-LD_DATA_SHIFT);
847 }
848 return x_m;
849 }
850
851
852
853
854