blob: 784f90907af100428d438e402293c73ae388b955
1 | /* ***** BEGIN LICENSE BLOCK ***** |
2 | * Source last modified: $Id: sbrmath.c,v 1.1 2005/02/26 01:47:35 jrecker Exp $ |
3 | * |
4 | * Portions Copyright (c) 1995-2005 RealNetworks, Inc. All Rights Reserved. |
5 | * |
6 | * The contents of this file, and the files included with this file, |
7 | * are subject to the current version of the RealNetworks Public |
8 | * Source License (the "RPSL") available at |
9 | * http://www.helixcommunity.org/content/rpsl unless you have licensed |
10 | * the file under the current version of the RealNetworks Community |
11 | * Source License (the "RCSL") available at |
12 | * http://www.helixcommunity.org/content/rcsl, in which case the RCSL |
13 | * will apply. You may also obtain the license terms directly from |
14 | * RealNetworks. You may not use this file except in compliance with |
15 | * the RPSL or, if you have a valid RCSL with RealNetworks applicable |
16 | * to this file, the RCSL. Please see the applicable RPSL or RCSL for |
17 | * the rights, obligations and limitations governing use of the |
18 | * contents of the file. |
19 | * |
20 | * This file is part of the Helix DNA Technology. RealNetworks is the |
21 | * developer of the Original Code and owns the copyrights in the |
22 | * portions it created. |
23 | * |
24 | * This file, and the files included with this file, is distributed |
25 | * and made available on an 'AS IS' basis, WITHOUT WARRANTY OF ANY |
26 | * KIND, EITHER EXPRESS OR IMPLIED, AND REALNETWORKS HEREBY DISCLAIMS |
27 | * ALL SUCH WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES |
28 | * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, QUIET |
29 | * ENJOYMENT OR NON-INFRINGEMENT. |
30 | * |
31 | * Technology Compatibility Kit Test Suite(s) Location: |
32 | * http://www.helixcommunity.org/content/tck |
33 | * |
34 | * Contributor(s): |
35 | * |
36 | * ***** END LICENSE BLOCK ***** */ |
37 | |
38 | /************************************************************************************** |
39 | * Fixed-point HE-AAC decoder |
40 | * Jon Recker (jrecker@real.com) |
41 | * February 2005 |
42 | * |
43 | * sbrmath.c - fixed-point math functions for SBR |
44 | **************************************************************************************/ |
45 | |
46 | #include "sbr.h" |
47 | #include "assembly.h" |
48 | |
49 | #define Q28_2 0x20000000 /* Q28: 2.0 */ |
50 | #define Q28_15 0x30000000 /* Q28: 1.5 */ |
51 | |
52 | #define NUM_ITER_IRN 5 |
53 | |
54 | /************************************************************************************** |
55 | * Function: InvRNormalized |
56 | * |
57 | * Description: use Newton's method to solve for x = 1/r |
58 | * |
59 | * Inputs: r = Q31, range = [0.5, 1) (normalize your inputs to this range) |
60 | * |
61 | * Outputs: none |
62 | * |
63 | * Return: x = Q29, range ~= [1.0, 2.0] |
64 | * |
65 | * Notes: guaranteed to converge and not overflow for any r in [0.5, 1) |
66 | * |
67 | * xn+1 = xn - f(xn)/f'(xn) |
68 | * f(x) = 1/r - x = 0 (find root) |
69 | * = 1/x - r |
70 | * f'(x) = -1/x^2 |
71 | * |
72 | * so xn+1 = xn - (1/xn - r) / (-1/xn^2) |
73 | * = xn * (2 - r*xn) |
74 | * |
75 | * NUM_ITER_IRN = 2, maxDiff = 6.2500e-02 (precision of about 4 bits) |
76 | * NUM_ITER_IRN = 3, maxDiff = 3.9063e-03 (precision of about 8 bits) |
77 | * NUM_ITER_IRN = 4, maxDiff = 1.5288e-05 (precision of about 16 bits) |
78 | * NUM_ITER_IRN = 5, maxDiff = 3.0034e-08 (precision of about 24 bits) |
79 | **************************************************************************************/ |
80 | int InvRNormalized(int r) |
81 | { |
82 | int i, xn, t; |
83 | |
84 | /* r = [0.5, 1.0) |
85 | * 1/r = (1.0, 2.0] |
86 | * so use 1.5 as initial guess |
87 | */ |
88 | xn = Q28_15; |
89 | |
90 | /* xn = xn*(2.0 - r*xn) */ |
91 | for (i = NUM_ITER_IRN; i != 0; i--) { |
92 | t = MULSHIFT32(r, xn); /* Q31*Q29 = Q28 */ |
93 | t = Q28_2 - t; /* Q28 */ |
94 | xn = MULSHIFT32(xn, t) << 4; /* Q29*Q28 << 4 = Q29 */ |
95 | } |
96 | |
97 | return xn; |
98 | } |
99 | |
100 | #define NUM_TERMS_RPI 5 |
101 | #define LOG2_EXP_INV 0x58b90bfc /* 1/log2(e), Q31 */ |
102 | |
103 | /* invTab[x] = 1/(x+1), format = Q30 */ |
104 | static const int invTab[NUM_TERMS_RPI] = {0x40000000, 0x20000000, 0x15555555, 0x10000000, 0x0ccccccd}; |
105 | |
106 | /************************************************************************************** |
107 | * Function: RatioPowInv |
108 | * |
109 | * Description: use Taylor (MacLaurin) series expansion to calculate (a/b) ^ (1/c) |
110 | * |
111 | * Inputs: a = [1, 64], b = [1, 64], c = [1, 64], a >= b |
112 | * |
113 | * Outputs: none |
114 | * |
115 | * Return: y = Q24, range ~= [0.015625, 64] |
116 | **************************************************************************************/ |
117 | int RatioPowInv(int a, int b, int c) |
118 | { |
119 | int lna, lnb, i, p, t, y; |
120 | |
121 | if (a < 1 || b < 1 || c < 1 || a > 64 || b > 64 || c > 64 || a < b) { |
122 | return 0; |
123 | } |
124 | |
125 | lna = MULSHIFT32(log2Tab[a], LOG2_EXP_INV) << 1; /* ln(a), Q28 */ |
126 | lnb = MULSHIFT32(log2Tab[b], LOG2_EXP_INV) << 1; /* ln(b), Q28 */ |
127 | p = (lna - lnb) / c; /* Q28 */ |
128 | |
129 | /* sum in Q24 */ |
130 | y = (1 << 24); |
131 | t = p >> 4; /* t = p^1 * 1/1! (Q24)*/ |
132 | y += t; |
133 | |
134 | for (i = 2; i <= NUM_TERMS_RPI; i++) { |
135 | t = MULSHIFT32(invTab[i - 1], t) << 2; |
136 | t = MULSHIFT32(p, t) << 4; /* t = p^i * 1/i! (Q24) */ |
137 | y += t; |
138 | } |
139 | |
140 | return y; |
141 | } |
142 | |
143 | /************************************************************************************** |
144 | * Function: SqrtFix |
145 | * |
146 | * Description: use binary search to calculate sqrt(q) |
147 | * |
148 | * Inputs: q = Q30 |
149 | * number of fraction bits in input |
150 | * |
151 | * Outputs: number of fraction bits in output |
152 | * |
153 | * Return: lo = Q(fBitsOut) |
154 | * |
155 | * Notes: absolute precision varies depending on fBitsIn |
156 | * normalizes input to range [0x200000000, 0x7fffffff] and takes |
157 | * floor(sqrt(input)), and sets fBitsOut appropriately |
158 | **************************************************************************************/ |
159 | int SqrtFix(int q, int fBitsIn, int *fBitsOut) |
160 | { |
161 | int z, lo, hi, mid; |
162 | |
163 | if (q <= 0) { |
164 | *fBitsOut = fBitsIn; |
165 | return 0; |
166 | } |
167 | |
168 | /* force even fBitsIn */ |
169 | z = fBitsIn & 0x01; |
170 | q >>= z; |
171 | fBitsIn -= z; |
172 | |
173 | /* for max precision, normalize to [0x20000000, 0x7fffffff] */ |
174 | z = (CLZ(q) - 1); |
175 | z >>= 1; |
176 | q <<= (2 * z); |
177 | |
178 | /* choose initial bounds */ |
179 | lo = 1; |
180 | if (q >= 0x10000000) { |
181 | lo = 16384; /* (int)sqrt(0x10000000) */ |
182 | } |
183 | hi = 46340; /* (int)sqrt(0x7fffffff) */ |
184 | |
185 | /* do binary search with 32x32->32 multiply test */ |
186 | do { |
187 | mid = (lo + hi) >> 1; |
188 | if (mid * mid > q) { |
189 | hi = mid - 1; |
190 | } else { |
191 | lo = mid + 1; |
192 | } |
193 | } while (hi >= lo); |
194 | lo--; |
195 | |
196 | *fBitsOut = ((fBitsIn + 2 * z) >> 1); |
197 | return lo; |
198 | } |
199 |