1
2/* -----------------------------------------------------------------------------------------------------------
3Software License for The Fraunhofer FDK AAC Codec Library for Android
4
5� Copyright  1995 - 2013 Fraunhofer-Gesellschaft zur F�rderung der angewandten Forschung e.V.
6  All rights reserved.
7
8 1.    INTRODUCTION
9The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software that implements
10the MPEG Advanced Audio Coding ("AAC") encoding and decoding scheme for digital audio.
11This FDK AAC Codec software is intended to be used on a wide variety of Android devices.
12
13AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient general perceptual
14audio codecs. AAC-ELD is considered the best-performing full-bandwidth communications codec by
15independent studies and is widely deployed. AAC has been standardized by ISO and IEC as part
16of the MPEG specifications.
17
18Patent licenses for necessary patent claims for the FDK AAC Codec (including those of Fraunhofer)
19may be obtained through Via Licensing (www.vialicensing.com) or through the respective patent owners
20individually for the purpose of encoding or decoding bit streams in products that are compliant with
21the ISO/IEC MPEG audio standards. Please note that most manufacturers of Android devices already license
22these patent claims through Via Licensing or directly from the patent owners, and therefore FDK AAC Codec
23software may already be covered under those patent licenses when it is used for those licensed purposes only.
24
25Commercially-licensed AAC software libraries, including floating-point versions with enhanced sound quality,
26are also available from Fraunhofer. Users are encouraged to check the Fraunhofer website for additional
27applications information and documentation.
28
292.    COPYRIGHT LICENSE
30
31Redistribution and use in source and binary forms, with or without modification, are permitted without
32payment of copyright license fees provided that you satisfy the following conditions:
33
34You must retain the complete text of this software license in redistributions of the FDK AAC Codec or
35your modifications thereto in source code form.
36
37You must retain the complete text of this software license in the documentation and/or other materials
38provided with redistributions of the FDK AAC Codec or your modifications thereto in binary form.
39You must make available free of charge copies of the complete source code of the FDK AAC Codec and your
40modifications thereto to recipients of copies in binary form.
41
42The name of Fraunhofer may not be used to endorse or promote products derived from this library without
43prior written permission.
44
45You may not charge copyright license fees for anyone to use, copy or distribute the FDK AAC Codec
46software or your modifications thereto.
47
48Your modified versions of the FDK AAC Codec must carry prominent notices stating that you changed the software
49and 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
533.    NO PATENT LICENSE
54
55NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without limitation the patents of Fraunhofer,
56ARE GRANTED BY THIS SOFTWARE LICENSE. Fraunhofer provides no warranty of patent non-infringement with
57respect to this software.
58
59You may use this FDK AAC Codec software or modifications thereto only for purposes that are authorized
60by appropriate patent licenses.
61
624.    DISCLAIMER
63
64This 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
66of merchantability and fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
67CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary, or consequential damages,
68including but not limited to procurement of substitute goods or services; loss of use, data, or profits,
69or business interruption, however caused and on any theory of liability, whether in contract, strict
70liability, or tort (including negligence), arising in any way out of the use of this software, even if
71advised of the possibility of such damage.
72
735.    CONTACT INFORMATION
74
75Fraunhofer Institute for Integrated Circuits IIS
76Attention: Audio and Multimedia Departments - FDK AAC LL
77Am Wolfsmantel 33
7891058 Erlangen, Germany
79
80www.iis.fraunhofer.de/amm
81amm-info@iis.fraunhofer.de
82----------------------------------------------------------------------------------------------------------- */
83
84/***************************  Fraunhofer IIS FDK Tools  **********************
85
86   Author(s):   Haricharan Lakshman, Manuel Jander
87   Description: Trigonometric functions fixed point fractional implementation.
88
89******************************************************************************/
90
91#include "FDK_trigFcts.h"
92
93#include "fixpoint_math.h"
94
95
96
97
98#define IMPROVE_ATAN2_ACCURACY  1  // 0 --> 59 dB SNR     1 --> 65 dB SNR
99#define MINSFTAB  7
100#define MAXSFTAB 25
101
102#if IMPROVE_ATAN2_ACCURACY
103static const FIXP_DBL f_atan_expand_range[MAXSFTAB-(MINSFTAB-1)]  =
104{
105  /*****************************************************************************
106   *
107   *  Table holds fixp_atan() output values which are outside of input range
108   *  of fixp_atan() to improve SNR of fixp_atan2().
109   *
110   *  This Table might also be used in fixp_atan() [todo] so there a wider input
111   *  range can be covered, too.
112   *
113   *  Matlab (generate table):
114   *    for scl = 7:25            % MINSFTAB .. MAXSFTAB
115   *      at=atan(0.5 *(2^scl));  % 0.5 because get in 'middle' area of current scale level 'scl'
116   *      at/2                    % div at by ATO_SCALE
117   *    end
118   *
119   *  Table divided by 2=ATO_SCALE  <--  SF=ATO_SF
120   *****************************************************************************/
121   FL2FXCONST_DBL(7.775862990872099e-001), FL2FXCONST_DBL(7.814919928673978e-001), FL2FXCONST_DBL(7.834450483314648e-001),
122   FL2FXCONST_DBL(7.844216021392089e-001), FL2FXCONST_DBL(7.849098823026687e-001), FL2FXCONST_DBL(7.851540227918509e-001),
123   FL2FXCONST_DBL(7.852760930873737e-001), FL2FXCONST_DBL(7.853371282415015e-001), FL2FXCONST_DBL(7.853676458193612e-001),
124   FL2FXCONST_DBL(7.853829046083906e-001), FL2FXCONST_DBL(7.853905340029177e-001), FL2FXCONST_DBL(7.853943487001828e-001),
125   FL2FXCONST_DBL(7.853962560488155e-001), FL2FXCONST_DBL(7.853972097231319e-001), FL2FXCONST_DBL(7.853976865602901e-001),
126   FL2FXCONST_DBL(7.853979249788692e-001), FL2FXCONST_DBL(7.853980441881587e-001), FL2FXCONST_DBL(7.853981037928035e-001),
127   FL2FXCONST_DBL(7.853981335951259e-001)
128   //     pi/4 = 0.785398163397448 = pi/2/ATO_SCALE
129};
130#endif
131
132FIXP_DBL fixp_atan2(FIXP_DBL y, FIXP_DBL x)
133{
134    FIXP_DBL q;
135    FIXP_DBL at;  // atan  out
136    FIXP_DBL at2; // atan2 out
137    FIXP_DBL ret = FL2FXCONST_DBL(-1.0f);
138    INT sf,sfo,stf;
139
140    // --- division
141
142    if      (y > FL2FXCONST_DBL(0.0f))
143    {
144        if      (x > FL2FXCONST_DBL(0.0f)) {
145                                           q =  fDivNormHighPrec( y, x, &sf); // both pos.
146        }
147        else if (x < FL2FXCONST_DBL(0.0f)) {
148                                           q = -fDivNormHighPrec( y,-x, &sf); // x neg.
149        }
150        else {//(x ==FL2FXCONST_DBL(0.0f))
151                                           q =  FL2FXCONST_DBL(+1.0f);  // y/x = pos/zero = +Inf
152                                           sf = 0;
153        }
154    }
155    else if (y < FL2FXCONST_DBL(0.0f))
156    {
157        if      (x > FL2FXCONST_DBL(0.0f)) {
158                                           q = -fDivNormHighPrec(-y, x, &sf); // y neg.
159        }
160        else if (x < FL2FXCONST_DBL(0.0f)) {
161                                           q =  fDivNormHighPrec(-y,-x, &sf); // both neg.
162        }
163        else {//(x ==FL2FXCONST_DBL(0.0f))
164                                           q =  FL2FXCONST_DBL(-1.0f);  // y/x = neg/zero = -Inf
165                                           sf = 0;
166        }
167    }
168    else { // (y ==FL2FXCONST_DBL(0.0f))
169        q = FL2FXCONST_DBL(0.0f);
170        sf = 0;
171    }
172    sfo = sf;
173
174    // --- atan()
175
176    if  ( sfo > ATI_SF ) {
177        // --- could not calc fixp_atan() here bec of input data out of range
178        //     ==> therefore give back boundary values
179
180        #if IMPROVE_ATAN2_ACCURACY
181        if (sfo > MAXSFTAB) sfo = MAXSFTAB;
182        #endif
183
184        if      (  q > FL2FXCONST_DBL(0.0f) ) {
185           #if IMPROVE_ATAN2_ACCURACY
186            at = +f_atan_expand_range[sfo-ATI_SF-1];
187           #else
188            at = FL2FXCONST_DBL( +M_PI/2 / ATO_SCALE);
189           #endif
190        }
191        else if (  q < FL2FXCONST_DBL(0.0f) ) {
192           #if IMPROVE_ATAN2_ACCURACY
193            at = -f_atan_expand_range[sfo-ATI_SF-1];
194           #else
195            at = FL2FXCONST_DBL( -M_PI/2 / ATO_SCALE);
196           #endif
197        }
198        else {  // q== FL2FXCONST_DBL(0.0f)
199            at = FL2FXCONST_DBL( 0.0f );
200        }
201    }else{
202        // --- calc of fixp_atan() is possible; input data within range
203        //     ==> set q on fixed scale level as desired from fixp_atan()
204        stf = sfo - ATI_SF;
205        if (stf > 0)  q = q << (INT)fMin( stf,DFRACT_BITS-1);
206        else          q = q >> (INT)fMin(-stf,DFRACT_BITS-1);
207        at = fixp_atan(q);  // ATO_SF
208    }
209
210    // --- atan2()
211
212    at2 = at >> (AT2O_SF - ATO_SF); // now AT2O_SF for atan2
213    if      (  x > FL2FXCONST_DBL(0.0f) ) {
214        ret = at2;
215    }
216    else if (  x < FL2FXCONST_DBL(0.0f) ) {
217        if (  y >= FL2FXCONST_DBL(0.0f) ) {
218            ret = at2 + FL2FXCONST_DBL( M_PI / AT2O_SCALE);
219        } else {
220            ret = at2 - FL2FXCONST_DBL( M_PI / AT2O_SCALE);
221        }
222    }
223    else {
224        // x == 0
225        if      ( y >  FL2FXCONST_DBL(0.0f) ) {
226            ret = FL2FXCONST_DBL( +M_PI/2 / AT2O_SCALE);
227        }
228        else if ( y <  FL2FXCONST_DBL(0.0f) ) {
229            ret = FL2FXCONST_DBL( -M_PI/2 / AT2O_SCALE);
230        }
231        else if ( y == FL2FXCONST_DBL(0.0f) ) {
232            ret = FL2FXCONST_DBL(0.0f);
233        }
234    }
235    return ret;
236}
237
238
239FIXP_DBL fixp_atan(FIXP_DBL x)
240{
241    INT sign;
242    FIXP_DBL result, temp;
243
244    // SNR of fixp_atan() = 56 dB
245    FIXP_DBL ONEBY3P56  = (FIXP_DBL)0x26800000; // 1.0/3.56 in q31
246    FIXP_DBL P281       = (FIXP_DBL)0x00013000; // 0.281 in q18
247    FIXP_DBL ONEP571    = (FIXP_DBL)0x6487ef00; // 1.571 in q30
248
249    if (x < FIXP_DBL(0)) {
250      sign = 1;
251      x = - x ;
252    } else {
253      sign = 0;
254    }
255
256    /* calc of arctan */
257    if(x < ( Q(Q_ATANINP)-FL2FXCONST_DBL(0.00395)) )
258    {
259        INT res_e;
260
261        temp = fPow2(x);            // q25 * q25 - (DFRACT_BITS-1) = q19
262        temp = fMult(temp, ONEBY3P56);      // q19 * q31 - (DFRACT_BITS-1) = q19
263        temp = temp + Q(19);                // q19 + q19 = q19
264        result = fDivNorm(x, temp, &res_e);
265        result = scaleValue(result, (Q_ATANOUT-Q_ATANINP+19-DFRACT_BITS+1) + res_e  );
266    }
267    else if( x < FL2FXCONST_DBL(1.28/64.0) )
268    {
269        FIXP_DBL delta_fix;
270        FIXP_DBL PI_BY_4 = FL2FXCONST_DBL(3.1415926/4.0) >> 1; /* pi/4 in q30 */
271
272        delta_fix = (x - FL2FXCONST_DBL(1.0/64.0)) << 5; /* q30 */
273        result = PI_BY_4 + (delta_fix >> 1) - (fPow2Div2(delta_fix));
274    }
275    else
276    {
277        INT res_e;
278
279        temp = fPow2Div2(x);        // q25 * q25 - (DFRACT_BITS-1) - 1 = q18
280        temp = temp + P281;                 // q18 + q18 = q18
281        result = fDivNorm(x, temp, &res_e);
282        result = scaleValue(result, (Q_ATANOUT-Q_ATANINP+18-DFRACT_BITS+1) + res_e );
283        result = ONEP571 - result;          // q30 + q30 = q30
284    }
285    if (sign) {
286      result = -result;
287    }
288
289    return(result);
290}
291
292
293
294#include "FDK_tools_rom.h"
295
296FIXP_DBL fixp_cos(FIXP_DBL x, int scale)
297{
298    FIXP_DBL residual, error, sine, cosine;
299
300    residual = fixp_sin_cos_residual_inline(x, scale, &sine, &cosine);
301    error = fMult(sine, residual);
302
303    return cosine - error;
304}
305
306FIXP_DBL fixp_sin(FIXP_DBL x, int scale)
307{
308    FIXP_DBL residual, error, sine, cosine;
309
310    residual = fixp_sin_cos_residual_inline(x, scale, &sine, &cosine);
311    error = fMult(cosine, residual);
312
313    return sine + error;
314}
315
316void fixp_cos_sin (FIXP_DBL x, int scale, FIXP_DBL *cos, FIXP_DBL *sin)
317{
318    FIXP_DBL residual, error0, error1, sine, cosine;
319
320    residual = fixp_sin_cos_residual_inline(x, scale, &sine, &cosine);
321    error0 = fMult(sine, residual);
322    error1 = fMult(cosine, residual);
323    *cos  = cosine - error0;
324    *sin  = sine + error1;
325}
326
327
328
329
330
331