/src/aac/libSBRdec/src/HFgen_preFlat.cpp
Line | Count | Source |
1 | | /* ----------------------------------------------------------------------------- |
2 | | Software License for The Fraunhofer FDK AAC Codec Library for Android |
3 | | |
4 | | © Copyright 1995 - 2019 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 | | /**************************** SBR decoder library ****************************** |
96 | | |
97 | | Author(s): Oliver Moser, Manuel Jander, Matthias Hildenbrand |
98 | | |
99 | | Description: QMF frequency pre-whitening for SBR. |
100 | | In the documentation the terms "scale factor" and "exponent" |
101 | | mean the same. Variables containing such information have |
102 | | the suffix "_sf". |
103 | | |
104 | | *******************************************************************************/ |
105 | | |
106 | | #include "HFgen_preFlat.h" |
107 | | |
108 | 12.8M | #define POLY_ORDER 3 |
109 | | #define MAXLOWBANDS 32 |
110 | | #define LOG10FAC 0.752574989159953f /* == 10/log2(10) * 2^-2 */ |
111 | | #define LOG10FAC_INV 0.664385618977472f /* == log2(10)/20 * 2^2 */ |
112 | | |
113 | | #define FIXP_CHB FIXP_SGL /* STB sinus Tab used in transformation */ |
114 | | #define CHC(a) (FX_DBL2FXCONST_SGL(a)) |
115 | 731k | #define FX_CHB2FX_DBL(a) FX_SGL2FX_DBL(a) |
116 | | |
117 | | typedef struct backsubst_data { |
118 | | FIXP_CHB Lnorm1d[3]; /*!< Normalized L matrix */ |
119 | | SCHAR Lnorm1d_sf[3]; |
120 | | FIXP_CHB Lnormii |
121 | | [3]; /*!< The diagonal data points [i][i] of the normalized L matrix */ |
122 | | SCHAR Lnormii_sf[3]; |
123 | | FIXP_CHB Bmul0 |
124 | | [4]; /*!< To normalize L*x=b, Bmul0 is what we need to multiply b with. */ |
125 | | SCHAR Bmul0_sf[4]; |
126 | | FIXP_CHB LnormInv1d[6]; /*!< Normalized inverted L matrix (L') */ |
127 | | SCHAR LnormInv1d_sf[6]; |
128 | | FIXP_CHB |
129 | | Bmul1[4]; /*!< To normalize L'*x=b, Bmul1 is what we need to multiply b |
130 | | with. */ |
131 | | SCHAR Bmul1_sf[4]; |
132 | | } backsubst_data; |
133 | | |
134 | | /* for each element n do, f(n) = trunc(log2(n))+1 */ |
135 | | const UCHAR getLog2[32] = {0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, |
136 | | 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5}; |
137 | | |
138 | | /** \def BSD_IDX_OFFSET |
139 | | * |
140 | | * bsd[] begins at index 0 with data for numBands=5. The correct bsd[] is |
141 | | * indexed like bsd[numBands-BSD_IDX_OFFSET]. |
142 | | */ |
143 | 365k | #define BSD_IDX_OFFSET 5 |
144 | | |
145 | | #define N_NUMBANDS \ |
146 | | MAXLOWBANDS - BSD_IDX_OFFSET + \ |
147 | | 1 /*!< Number of backsubst_data elements in bsd */ |
148 | | |
149 | | const backsubst_data bsd[N_NUMBANDS] = { |
150 | | { |
151 | | /* numBands=5 */ |
152 | | {CHC(0x66c85a52), CHC(0x4278e587), CHC(0x697dcaff)}, |
153 | | {-1, 0, 0}, |
154 | | {CHC(0x66a61789), CHC(0x5253b8e3), CHC(0x5addad81)}, |
155 | | {3, 4, 1}, |
156 | | {CHC(0x7525ee90), CHC(0x6e2a1210), CHC(0x6523bb40), CHC(0x59822ead)}, |
157 | | {-6, -4, -2, 0}, |
158 | | {CHC(0x609e4cad), CHC(0x59c7e312), CHC(0x681eecac), CHC(0x440ea893), |
159 | | CHC(0x4a214bb3), CHC(0x53c345a1)}, |
160 | | {1, 0, -1, -1, -3, -5}, |
161 | | {CHC(0x7525ee90), CHC(0x58587936), CHC(0x410d0b38), CHC(0x7f1519d6)}, |
162 | | {-6, -1, 2, 0}, |
163 | | }, |
164 | | { |
165 | | /* numBands=6 */ |
166 | | {CHC(0x68943285), CHC(0x4841d2c3), CHC(0x6a6214c7)}, |
167 | | {-1, 0, 0}, |
168 | | {CHC(0x63c5923e), CHC(0x4e906e18), CHC(0x6285af8a)}, |
169 | | {3, 4, 1}, |
170 | | {CHC(0x7263940b), CHC(0x424a69a5), CHC(0x4ae8383a), CHC(0x517b7730)}, |
171 | | {-7, -4, -2, 0}, |
172 | | {CHC(0x518aee5f), CHC(0x4823a096), CHC(0x43764a39), CHC(0x6e6faf23), |
173 | | CHC(0x61bba44f), CHC(0x59d8b132)}, |
174 | | {1, 0, -1, -2, -4, -6}, |
175 | | {CHC(0x7263940b), CHC(0x6757bff2), CHC(0x5bf40fe0), CHC(0x7d6f4292)}, |
176 | | {-7, -2, 1, 0}, |
177 | | }, |
178 | | { |
179 | | /* numBands=7 */ |
180 | | {CHC(0x699b4c3c), CHC(0x4b8b702f), CHC(0x6ae51a4f)}, |
181 | | {-1, 0, 0}, |
182 | | {CHC(0x623a7f49), CHC(0x4ccc91fc), CHC(0x68f048dd)}, |
183 | | {3, 4, 1}, |
184 | | {CHC(0x7e6ebe18), CHC(0x5701daf2), CHC(0x74a8198b), CHC(0x4b399aa1)}, |
185 | | {-8, -5, -3, 0}, |
186 | | {CHC(0x464a64a6), CHC(0x78e42633), CHC(0x5ee174ba), CHC(0x5d0008c8), |
187 | | CHC(0x455cff0f), CHC(0x6b9100e7)}, |
188 | | {1, -1, -2, -2, -4, -7}, |
189 | | {CHC(0x7e6ebe18), CHC(0x42c52efe), CHC(0x45fe401f), CHC(0x7b5808ef)}, |
190 | | {-8, -2, 1, 0}, |
191 | | }, |
192 | | { |
193 | | /* numBands=8 */ |
194 | | {CHC(0x6a3fd9b4), CHC(0x4d99823f), CHC(0x6b372a94)}, |
195 | | {-1, 0, 0}, |
196 | | {CHC(0x614c6ef7), CHC(0x4bd06699), CHC(0x6e59cfca)}, |
197 | | {3, 4, 1}, |
198 | | {CHC(0x4c389cc5), CHC(0x79686681), CHC(0x5e2544c2), CHC(0x46305b43)}, |
199 | | {-8, -6, -3, 0}, |
200 | | {CHC(0x7b4ca7c6), CHC(0x68270ac5), CHC(0x467c644c), CHC(0x505c1b0f), |
201 | | CHC(0x67a14778), CHC(0x45801767)}, |
202 | | {0, -1, -2, -2, -5, -7}, |
203 | | {CHC(0x4c389cc5), CHC(0x5c499ceb), CHC(0x6f863c9f), CHC(0x79059bfc)}, |
204 | | {-8, -3, 0, 0}, |
205 | | }, |
206 | | { |
207 | | /* numBands=9 */ |
208 | | {CHC(0x6aad9988), CHC(0x4ef8ac18), CHC(0x6b6df116)}, |
209 | | {-1, 0, 0}, |
210 | | {CHC(0x60b159b0), CHC(0x4b33f772), CHC(0x72f5573d)}, |
211 | | {3, 4, 1}, |
212 | | {CHC(0x6206cb18), CHC(0x58a7d8dc), CHC(0x4e0b2d0b), CHC(0x4207ad84)}, |
213 | | {-9, -6, -3, 0}, |
214 | | {CHC(0x6dadadae), CHC(0x5b8b2cfc), CHC(0x6cf61db2), CHC(0x46c3c90b), |
215 | | CHC(0x506314ea), CHC(0x5f034acd)}, |
216 | | {0, -1, -3, -2, -5, -8}, |
217 | | {CHC(0x6206cb18), CHC(0x42f8b8de), CHC(0x5bb4776f), CHC(0x769acc79)}, |
218 | | {-9, -3, 0, 0}, |
219 | | }, |
220 | | { |
221 | | /* numBands=10 */ |
222 | | {CHC(0x6afa7252), CHC(0x4feed3ed), CHC(0x6b94504d)}, |
223 | | {-1, 0, 0}, |
224 | | {CHC(0x60467899), CHC(0x4acbafba), CHC(0x76eb327f)}, |
225 | | {3, 4, 1}, |
226 | | {CHC(0x42415b15), CHC(0x431080da), CHC(0x420f1c32), CHC(0x7d0c1aeb)}, |
227 | | {-9, -6, -3, -1}, |
228 | | {CHC(0x62b2c7a4), CHC(0x51b040a6), CHC(0x56caddb4), CHC(0x7e74a2c8), |
229 | | CHC(0x4030adf5), CHC(0x43d1dc4f)}, |
230 | | {0, -1, -3, -3, -5, -8}, |
231 | | {CHC(0x42415b15), CHC(0x64e299b3), CHC(0x4d33b5e8), CHC(0x742cee5f)}, |
232 | | {-9, -4, 0, 0}, |
233 | | }, |
234 | | { |
235 | | /* numBands=11 */ |
236 | | {CHC(0x6b3258bb), CHC(0x50a21233), CHC(0x6bb03c19)}, |
237 | | {-1, 0, 0}, |
238 | | {CHC(0x5ff997c6), CHC(0x4a82706e), CHC(0x7a5aae36)}, |
239 | | {3, 4, 1}, |
240 | | {CHC(0x5d2fb4fb), CHC(0x685bddd8), CHC(0x71b5e983), CHC(0x7708c90b)}, |
241 | | {-10, -7, -4, -1}, |
242 | | {CHC(0x59aceea2), CHC(0x49c428a0), CHC(0x46ca5527), CHC(0x724be884), |
243 | | CHC(0x68e586da), CHC(0x643485b6)}, |
244 | | {0, -1, -3, -3, -6, -9}, |
245 | | {CHC(0x5d2fb4fb), CHC(0x4e3fad1a), CHC(0x42310ba2), CHC(0x71c8b3ce)}, |
246 | | {-10, -4, 0, 0}, |
247 | | }, |
248 | | { |
249 | | /* numBands=12 */ |
250 | | {CHC(0x6b5c4726), CHC(0x5128a4a8), CHC(0x6bc52ee1)}, |
251 | | {-1, 0, 0}, |
252 | | {CHC(0x5fc06618), CHC(0x4a4ce559), CHC(0x7d5c16e9)}, |
253 | | {3, 4, 1}, |
254 | | {CHC(0x43af8342), CHC(0x531533d3), CHC(0x633660a6), CHC(0x71ce6052)}, |
255 | | {-10, -7, -4, -1}, |
256 | | {CHC(0x522373d7), CHC(0x434150cb), CHC(0x75b58afc), CHC(0x68474f2d), |
257 | | CHC(0x575348a5), CHC(0x4c20973f)}, |
258 | | {0, -1, -4, -3, -6, -9}, |
259 | | {CHC(0x43af8342), CHC(0x7c4d3d11), CHC(0x732e13db), CHC(0x6f756ac4)}, |
260 | | {-10, -5, -1, 0}, |
261 | | }, |
262 | | { |
263 | | /* numBands=13 */ |
264 | | {CHC(0x6b7c8953), CHC(0x51903fcd), CHC(0x6bd54d2e)}, |
265 | | {-1, 0, 0}, |
266 | | {CHC(0x5f94abf0), CHC(0x4a2480fa), CHC(0x40013553)}, |
267 | | {3, 4, 2}, |
268 | | {CHC(0x6501236e), CHC(0x436b9c4e), CHC(0x578d7881), CHC(0x6d34f92e)}, |
269 | | {-11, -7, -4, -1}, |
270 | | {CHC(0x4bc0e2b2), CHC(0x7b9d12ac), CHC(0x636c1c1b), CHC(0x5fe15c2b), |
271 | | CHC(0x49d54879), CHC(0x7662cfa5)}, |
272 | | {0, -2, -4, -3, -6, -10}, |
273 | | {CHC(0x6501236e), CHC(0x64b059fe), CHC(0x656d8359), CHC(0x6d370900)}, |
274 | | {-11, -5, -1, 0}, |
275 | | }, |
276 | | { |
277 | | /* numBands=14 */ |
278 | | {CHC(0x6b95e276), CHC(0x51e1b637), CHC(0x6be1f7ed)}, |
279 | | {-1, 0, 0}, |
280 | | {CHC(0x5f727a1c), CHC(0x4a053e9c), CHC(0x412e528c)}, |
281 | | {3, 4, 2}, |
282 | | {CHC(0x4d178bd4), CHC(0x6f33b4e8), CHC(0x4e028f7f), CHC(0x691ee104)}, |
283 | | {-11, -8, -4, -1}, |
284 | | {CHC(0x46473d3f), CHC(0x725bd0a6), CHC(0x55199885), CHC(0x58bcc56b), |
285 | | CHC(0x7e7e6288), CHC(0x5ddef6eb)}, |
286 | | {0, -2, -4, -3, -7, -10}, |
287 | | {CHC(0x4d178bd4), CHC(0x52ebd467), CHC(0x5a395a6e), CHC(0x6b0f724f)}, |
288 | | {-11, -5, -1, 0}, |
289 | | }, |
290 | | { |
291 | | /* numBands=15 */ |
292 | | {CHC(0x6baa2a22), CHC(0x5222eb91), CHC(0x6bec1a86)}, |
293 | | {-1, 0, 0}, |
294 | | {CHC(0x5f57393b), CHC(0x49ec8934), CHC(0x423b5b58)}, |
295 | | {3, 4, 2}, |
296 | | {CHC(0x77fd2486), CHC(0x5cfbdf2c), CHC(0x46153bd1), CHC(0x65757ed9)}, |
297 | | {-12, -8, -4, -1}, |
298 | | {CHC(0x41888ee6), CHC(0x6a661db3), CHC(0x49abc8c8), CHC(0x52965848), |
299 | | CHC(0x6d9301b7), CHC(0x4bb04721)}, |
300 | | {0, -2, -4, -3, -7, -10}, |
301 | | {CHC(0x77fd2486), CHC(0x45424c68), CHC(0x50f33cc6), CHC(0x68ff43f0)}, |
302 | | {-12, -5, -1, 0}, |
303 | | }, |
304 | | { |
305 | | /* numBands=16 */ |
306 | | {CHC(0x6bbaa499), CHC(0x5257ed94), CHC(0x6bf456e4)}, |
307 | | {-1, 0, 0}, |
308 | | {CHC(0x5f412594), CHC(0x49d8a766), CHC(0x432d1dbd)}, |
309 | | {3, 4, 2}, |
310 | | {CHC(0x5ef5cfde), CHC(0x4eafcd2d), CHC(0x7ed36893), CHC(0x62274b45)}, |
311 | | {-12, -8, -5, -1}, |
312 | | {CHC(0x7ac438f5), CHC(0x637aab21), CHC(0x4067617a), CHC(0x4d3c6ec7), |
313 | | CHC(0x5fd6e0dd), CHC(0x7bd5f024)}, |
314 | | {-1, -2, -4, -3, -7, -11}, |
315 | | {CHC(0x5ef5cfde), CHC(0x751d0d4f), CHC(0x492b3c41), CHC(0x67065409)}, |
316 | | {-12, -6, -1, 0}, |
317 | | }, |
318 | | { |
319 | | /* numBands=17 */ |
320 | | {CHC(0x6bc836c9), CHC(0x5283997e), CHC(0x6bfb1f5e)}, |
321 | | {-1, 0, 0}, |
322 | | {CHC(0x5f2f02b6), CHC(0x49c868e9), CHC(0x44078151)}, |
323 | | {3, 4, 2}, |
324 | | {CHC(0x4c43b65a), CHC(0x4349dcf6), CHC(0x73799e2d), CHC(0x5f267274)}, |
325 | | {-12, -8, -5, -1}, |
326 | | {CHC(0x73726394), CHC(0x5d68511a), CHC(0x7191bbcc), CHC(0x48898c70), |
327 | | CHC(0x548956e1), CHC(0x66981ce8)}, |
328 | | {-1, -2, -5, -3, -7, -11}, |
329 | | {CHC(0x4c43b65a), CHC(0x64131116), CHC(0x429028e2), CHC(0x65240211)}, |
330 | | {-12, -6, -1, 0}, |
331 | | }, |
332 | | { |
333 | | /* numBands=18 */ |
334 | | {CHC(0x6bd3860d), CHC(0x52a80156), CHC(0x6c00c68d)}, |
335 | | {-1, 0, 0}, |
336 | | {CHC(0x5f1fed86), CHC(0x49baf636), CHC(0x44cdb9dc)}, |
337 | | {3, 4, 2}, |
338 | | {CHC(0x7c189389), CHC(0x742666d8), CHC(0x69b8c776), CHC(0x5c67e27d)}, |
339 | | {-13, -9, -5, -1}, |
340 | | {CHC(0x6cf1ea76), CHC(0x58095703), CHC(0x64e351a9), CHC(0x4460da90), |
341 | | CHC(0x4b1f8083), CHC(0x55f2d3e1)}, |
342 | | {-1, -2, -5, -3, -7, -11}, |
343 | | {CHC(0x7c189389), CHC(0x5651792a), CHC(0x79cb9b3d), CHC(0x635769c0)}, |
344 | | {-13, -6, -2, 0}, |
345 | | }, |
346 | | { |
347 | | /* numBands=19 */ |
348 | | {CHC(0x6bdd0c40), CHC(0x52c6abf6), CHC(0x6c058950)}, |
349 | | {-1, 0, 0}, |
350 | | {CHC(0x5f133f88), CHC(0x49afb305), CHC(0x45826d73)}, |
351 | | {3, 4, 2}, |
352 | | {CHC(0x6621a164), CHC(0x6512528e), CHC(0x61449fc8), CHC(0x59e2a0c0)}, |
353 | | {-13, -9, -5, -1}, |
354 | | {CHC(0x6721cadb), CHC(0x53404cd4), CHC(0x5a389e91), CHC(0x40abcbd2), |
355 | | CHC(0x43332f01), CHC(0x48b82e46)}, |
356 | | {-1, -2, -5, -3, -7, -11}, |
357 | | {CHC(0x6621a164), CHC(0x4b12cc28), CHC(0x6ffd4df8), CHC(0x619f835e)}, |
358 | | {-13, -6, -2, 0}, |
359 | | }, |
360 | | { |
361 | | /* numBands=20 */ |
362 | | {CHC(0x6be524c5), CHC(0x52e0beb3), CHC(0x6c099552)}, |
363 | | {-1, 0, 0}, |
364 | | {CHC(0x5f087c68), CHC(0x49a62bb5), CHC(0x4627d175)}, |
365 | | {3, 4, 2}, |
366 | | {CHC(0x54ec6afe), CHC(0x58991a42), CHC(0x59e23e8c), CHC(0x578f4ef4)}, |
367 | | {-13, -9, -5, -1}, |
368 | | {CHC(0x61e78f6f), CHC(0x4ef5e1e9), CHC(0x5129c3b8), CHC(0x7ab0f7b2), |
369 | | CHC(0x78efb076), CHC(0x7c2567ea)}, |
370 | | {-1, -2, -5, -4, -8, -12}, |
371 | | {CHC(0x54ec6afe), CHC(0x41c7812c), CHC(0x676f6f8d), CHC(0x5ffb383f)}, |
372 | | {-13, -6, -2, 0}, |
373 | | }, |
374 | | { |
375 | | /* numBands=21 */ |
376 | | {CHC(0x6bec1542), CHC(0x52f71929), CHC(0x6c0d0d5e)}, |
377 | | {-1, 0, 0}, |
378 | | {CHC(0x5eff45c5), CHC(0x499e092d), CHC(0x46bfc0c9)}, |
379 | | {3, 4, 2}, |
380 | | {CHC(0x47457a78), CHC(0x4e2d99b3), CHC(0x53637ea5), CHC(0x5567d0e9)}, |
381 | | {-13, -9, -5, -1}, |
382 | | {CHC(0x5d2dc61b), CHC(0x4b1760c8), CHC(0x4967cf39), CHC(0x74b113d8), |
383 | | CHC(0x6d6676b6), CHC(0x6ad114e9)}, |
384 | | {-1, -2, -5, -4, -8, -12}, |
385 | | {CHC(0x47457a78), CHC(0x740accaa), CHC(0x5feb6609), CHC(0x5e696f95)}, |
386 | | {-13, -7, -2, 0}, |
387 | | }, |
388 | | { |
389 | | /* numBands=22 */ |
390 | | {CHC(0x6bf21387), CHC(0x530a683c), CHC(0x6c100c59)}, |
391 | | {-1, 0, 0}, |
392 | | {CHC(0x5ef752ea), CHC(0x499708c6), CHC(0x474bcd1b)}, |
393 | | {3, 4, 2}, |
394 | | {CHC(0x78a21ab7), CHC(0x45658aec), CHC(0x4da3c4fe), CHC(0x5367094b)}, |
395 | | {-14, -9, -5, -1}, |
396 | | {CHC(0x58e2df6a), CHC(0x4795990e), CHC(0x42b5e0f7), CHC(0x6f408c64), |
397 | | CHC(0x6370bebf), CHC(0x5c91ca85)}, |
398 | | {-1, -2, -5, -4, -8, -12}, |
399 | | {CHC(0x78a21ab7), CHC(0x66f951d6), CHC(0x594605bb), CHC(0x5ce91657)}, |
400 | | {-14, -7, -2, 0}, |
401 | | }, |
402 | | { |
403 | | /* numBands=23 */ |
404 | | {CHC(0x6bf749b2), CHC(0x531b3348), CHC(0x6c12a750)}, |
405 | | {-1, 0, 0}, |
406 | | {CHC(0x5ef06b17), CHC(0x4990f6c9), CHC(0x47cd4c5b)}, |
407 | | {3, 4, 2}, |
408 | | {CHC(0x66dede36), CHC(0x7bdf90a9), CHC(0x4885b2b9), CHC(0x5188a6b7)}, |
409 | | {-14, -10, -5, -1}, |
410 | | {CHC(0x54f85812), CHC(0x446414ae), CHC(0x79c8d519), CHC(0x6a4c2f31), |
411 | | CHC(0x5ac8325f), CHC(0x50bf9200)}, |
412 | | {-1, -2, -6, -4, -8, -12}, |
413 | | {CHC(0x66dede36), CHC(0x5be0d90e), CHC(0x535cc453), CHC(0x5b7923f0)}, |
414 | | {-14, -7, -2, 0}, |
415 | | }, |
416 | | { |
417 | | /* numBands=24 */ |
418 | | {CHC(0x6bfbd91d), CHC(0x5329e580), CHC(0x6c14eeed)}, |
419 | | {-1, 0, 0}, |
420 | | {CHC(0x5eea6179), CHC(0x498baa90), CHC(0x4845635d)}, |
421 | | {3, 4, 2}, |
422 | | {CHC(0x58559b7e), CHC(0x6f1b231f), CHC(0x43f1789b), CHC(0x4fc8fcb8)}, |
423 | | {-14, -10, -5, -1}, |
424 | | {CHC(0x51621775), CHC(0x417881a3), CHC(0x6f9ba9b6), CHC(0x65c412b2), |
425 | | CHC(0x53352c61), CHC(0x46db9caf)}, |
426 | | {-1, -2, -6, -4, -8, -12}, |
427 | | {CHC(0x58559b7e), CHC(0x52636003), CHC(0x4e13b316), CHC(0x5a189cdf)}, |
428 | | {-14, -7, -2, 0}, |
429 | | }, |
430 | | { |
431 | | /* numBands=25 */ |
432 | | {CHC(0x6bffdc73), CHC(0x5336d4af), CHC(0x6c16f084)}, |
433 | | {-1, 0, 0}, |
434 | | {CHC(0x5ee51249), CHC(0x498703cc), CHC(0x48b50e4f)}, |
435 | | {3, 4, 2}, |
436 | | {CHC(0x4c5616cf), CHC(0x641b9fad), CHC(0x7fa735e0), CHC(0x4e24e57a)}, |
437 | | {-14, -10, -6, -1}, |
438 | | {CHC(0x4e15f47a), CHC(0x7d9481d6), CHC(0x66a82f8a), CHC(0x619ae971), |
439 | | CHC(0x4c8b2f5f), CHC(0x7d09ec11)}, |
440 | | {-1, -3, -6, -4, -8, -13}, |
441 | | {CHC(0x4c5616cf), CHC(0x4a3770fb), CHC(0x495402de), CHC(0x58c693fa)}, |
442 | | {-14, -7, -2, 0}, |
443 | | }, |
444 | | { |
445 | | /* numBands=26 */ |
446 | | {CHC(0x6c036943), CHC(0x53424625), CHC(0x6c18b6dc)}, |
447 | | {-1, 0, 0}, |
448 | | {CHC(0x5ee060aa), CHC(0x4982e88a), CHC(0x491d277f)}, |
449 | | {3, 4, 2}, |
450 | | {CHC(0x425ada5b), CHC(0x5a9368ac), CHC(0x78380a42), CHC(0x4c99aa05)}, |
451 | | {-14, -10, -6, -1}, |
452 | | {CHC(0x4b0b569c), CHC(0x78a420da), CHC(0x5ebdf203), CHC(0x5dc57e63), |
453 | | CHC(0x46a650ff), CHC(0x6ee13fb8)}, |
454 | | {-1, -3, -6, -4, -8, -13}, |
455 | | {CHC(0x425ada5b), CHC(0x4323073c), CHC(0x450ae92b), CHC(0x57822ad5)}, |
456 | | {-14, -7, -2, 0}, |
457 | | }, |
458 | | { |
459 | | /* numBands=27 */ |
460 | | {CHC(0x6c06911a), CHC(0x534c7261), CHC(0x6c1a4aba)}, |
461 | | {-1, 0, 0}, |
462 | | {CHC(0x5edc3524), CHC(0x497f43c0), CHC(0x497e6cd8)}, |
463 | | {3, 4, 2}, |
464 | | {CHC(0x73fb550e), CHC(0x5244894f), CHC(0x717aad78), CHC(0x4b24ef6c)}, |
465 | | {-15, -10, -6, -1}, |
466 | | {CHC(0x483aebe4), CHC(0x74139116), CHC(0x57b58037), CHC(0x5a3a4f3c), |
467 | | CHC(0x416950fe), CHC(0x62c7f4f2)}, |
468 | | {-1, -3, -6, -4, -8, -13}, |
469 | | {CHC(0x73fb550e), CHC(0x79efb994), CHC(0x4128cab7), CHC(0x564a919a)}, |
470 | | {-15, -8, -2, 0}, |
471 | | }, |
472 | | { |
473 | | /* numBands=28 */ |
474 | | {CHC(0x6c096264), CHC(0x535587cd), CHC(0x6c1bb355)}, |
475 | | {-1, 0, 0}, |
476 | | {CHC(0x5ed87c76), CHC(0x497c0439), CHC(0x49d98452)}, |
477 | | {3, 4, 2}, |
478 | | {CHC(0x65dec5bf), CHC(0x4afd1ba3), CHC(0x6b58b4b3), CHC(0x49c4a7b0)}, |
479 | | {-15, -10, -6, -1}, |
480 | | {CHC(0x459e6eb1), CHC(0x6fd850b7), CHC(0x516e7be9), CHC(0x56f13d05), |
481 | | CHC(0x79785594), CHC(0x58617de7)}, |
482 | | {-1, -3, -6, -4, -9, -13}, |
483 | | {CHC(0x65dec5bf), CHC(0x6f2168aa), CHC(0x7b41310f), CHC(0x551f0692)}, |
484 | | {-15, -8, -3, 0}, |
485 | | }, |
486 | | { |
487 | | /* numBands=29 */ |
488 | | {CHC(0x6c0be913), CHC(0x535dacd5), CHC(0x6c1cf6a3)}, |
489 | | {-1, 0, 0}, |
490 | | {CHC(0x5ed526b4), CHC(0x49791bc5), CHC(0x4a2eff99)}, |
491 | | {3, 4, 2}, |
492 | | {CHC(0x59e44afe), CHC(0x44949ada), CHC(0x65bf36f5), CHC(0x487705a0)}, |
493 | | {-15, -10, -6, -1}, |
494 | | {CHC(0x43307779), CHC(0x6be959c4), CHC(0x4bce2122), CHC(0x53e34d89), |
495 | | CHC(0x7115ff82), CHC(0x4f6421a1)}, |
496 | | {-1, -3, -6, -4, -9, -13}, |
497 | | {CHC(0x59e44afe), CHC(0x659eab7d), CHC(0x74cea459), CHC(0x53fed574)}, |
498 | | {-15, -8, -3, 0}, |
499 | | }, |
500 | | { |
501 | | /* numBands=30 */ |
502 | | {CHC(0x6c0e2f17), CHC(0x53650181), CHC(0x6c1e199d)}, |
503 | | {-1, 0, 0}, |
504 | | {CHC(0x5ed2269f), CHC(0x49767e9e), CHC(0x4a7f5f0b)}, |
505 | | {3, 4, 2}, |
506 | | {CHC(0x4faa4ae6), CHC(0x7dd3bf11), CHC(0x609e2732), CHC(0x473a72e9)}, |
507 | | {-15, -11, -6, -1}, |
508 | | {CHC(0x40ec57c6), CHC(0x683ee147), CHC(0x46be261d), CHC(0x510a7983), |
509 | | CHC(0x698a84cb), CHC(0x4794a927)}, |
510 | | {-1, -3, -6, -4, -9, -13}, |
511 | | {CHC(0x4faa4ae6), CHC(0x5d3615ad), CHC(0x6ee74773), CHC(0x52e956a1)}, |
512 | | {-15, -8, -3, 0}, |
513 | | }, |
514 | | { |
515 | | /* numBands=31 */ |
516 | | {CHC(0x6c103cc9), CHC(0x536ba0ac), CHC(0x6c1f2070)}, |
517 | | {-1, 0, 0}, |
518 | | {CHC(0x5ecf711e), CHC(0x497422ea), CHC(0x4acb1438)}, |
519 | | {3, 4, 2}, |
520 | | {CHC(0x46e322ad), CHC(0x73c32f3c), CHC(0x5be7d172), CHC(0x460d8800)}, |
521 | | {-15, -11, -6, -1}, |
522 | | {CHC(0x7d9bf8ad), CHC(0x64d22351), CHC(0x422bdc81), CHC(0x4e6184aa), |
523 | | CHC(0x62ba2375), CHC(0x40c325de)}, |
524 | | {-2, -3, -6, -4, -9, -13}, |
525 | | {CHC(0x46e322ad), CHC(0x55bef2a3), CHC(0x697b3135), CHC(0x51ddee4d)}, |
526 | | {-15, -8, -3, 0}, |
527 | | }, |
528 | | { |
529 | | // numBands=32 |
530 | | {CHC(0x6c121933), CHC(0x5371a104), CHC(0x6c200ea0)}, |
531 | | {-1, 0, 0}, |
532 | | {CHC(0x5eccfcd3), CHC(0x49720060), CHC(0x4b1283f0)}, |
533 | | {3, 4, 2}, |
534 | | {CHC(0x7ea12a52), CHC(0x6aca3303), CHC(0x579072bf), CHC(0x44ef056e)}, |
535 | | {-16, -11, -6, -1}, |
536 | | {CHC(0x79a3a9ab), CHC(0x619d38fc), CHC(0x7c0f0734), CHC(0x4be3dd5d), |
537 | | CHC(0x5c8d7163), CHC(0x7591065f)}, |
538 | | {-2, -3, -7, -4, -9, -14}, |
539 | | {CHC(0x7ea12a52), CHC(0x4f1782a6), CHC(0x647cbcb2), CHC(0x50dc0bb1)}, |
540 | | {-16, -8, -3, 0}, |
541 | | }, |
542 | | }; |
543 | | |
544 | | /** \def SUM_SAFETY |
545 | | * |
546 | | * SUM_SAFTEY defines the bits needed to right-shift every summand in |
547 | | * order to be overflow-safe. In the two backsubst functions we sum up 4 |
548 | | * values. Since one of which is definitely not MAXVAL_DBL (the L[x][y]), |
549 | | * we spare just 2 safety bits instead of 3. |
550 | | */ |
551 | 877k | #define SUM_SAFETY 2 |
552 | | |
553 | | /** |
554 | | * \brief Solves L*x=b via backsubstitution according to the following |
555 | | * structure: |
556 | | * |
557 | | * x[0] = b[0]; |
558 | | * x[1] = (b[1] - x[0]) / L[1][1]; |
559 | | * x[2] = (b[2] - x[1]*L[2][1] - x[0]) / L[2][2]; |
560 | | * x[3] = (b[3] - x[2]*L[3][2] - x[1]*L[3][1] - x[0]) / L[3][3]; |
561 | | * |
562 | | * \param[in] numBands SBR crossover band index |
563 | | * \param[in] b the b in L*x=b (one-dimensional) |
564 | | * \param[out] x output polynomial coefficients (mantissa) |
565 | | * \param[out] x_sf exponents of x[] |
566 | | */ |
567 | | static void backsubst_fw(const int numBands, const FIXP_DBL *const b, |
568 | 36.5k | FIXP_DBL *RESTRICT x, int *RESTRICT x_sf) { |
569 | 36.5k | int i, k; |
570 | 36.5k | int m; /* the trip counter that indexes incrementally through Lnorm1d[] */ |
571 | | |
572 | 36.5k | const FIXP_CHB *RESTRICT pLnorm1d = bsd[numBands - BSD_IDX_OFFSET].Lnorm1d; |
573 | 36.5k | const SCHAR *RESTRICT pLnorm1d_sf = bsd[numBands - BSD_IDX_OFFSET].Lnorm1d_sf; |
574 | 36.5k | const FIXP_CHB *RESTRICT pLnormii = bsd[numBands - BSD_IDX_OFFSET].Lnormii; |
575 | 36.5k | const SCHAR *RESTRICT pLnormii_sf = bsd[numBands - BSD_IDX_OFFSET].Lnormii_sf; |
576 | | |
577 | 36.5k | x[0] = b[0]; |
578 | | |
579 | 146k | for (i = 1, m = 0; i <= POLY_ORDER; ++i) { |
580 | 109k | FIXP_DBL sum = b[i] >> SUM_SAFETY; |
581 | 109k | int sum_sf = x_sf[i]; |
582 | 219k | for (k = i - 1; k > 0; --k, ++m) { |
583 | 109k | int e; |
584 | 109k | FIXP_DBL mult = fMultNorm(FX_CHB2FX_DBL(pLnorm1d[m]), x[k], &e); |
585 | 109k | int mult_sf = pLnorm1d_sf[m] + x_sf[k] + e; |
586 | | |
587 | | /* check if the new summand mult has a different sf than the sum currently |
588 | | * has */ |
589 | 109k | int diff = mult_sf - sum_sf; |
590 | | |
591 | 109k | if (diff > 0) { |
592 | | /* yes, and it requires the sum to be adjusted (scaled down) */ |
593 | 48.5k | sum >>= diff; |
594 | 48.5k | sum_sf = mult_sf; |
595 | 61.1k | } else if (diff < 0) { |
596 | | /* yes, but here mult needs to be scaled down */ |
597 | 55.0k | mult >>= -diff; |
598 | 55.0k | } |
599 | 109k | sum -= (mult >> SUM_SAFETY); |
600 | 109k | } |
601 | | |
602 | | /* - x[0] */ |
603 | 109k | if (x_sf[0] > sum_sf) { |
604 | 32.2k | sum >>= (x_sf[0] - sum_sf); |
605 | 32.2k | sum_sf = x_sf[0]; |
606 | 32.2k | } |
607 | 109k | sum -= (x[0] >> (sum_sf - x_sf[0] + SUM_SAFETY)); |
608 | | |
609 | | /* instead of the division /L[i][i], we multiply by the inverse */ |
610 | 109k | int e; |
611 | 109k | x[i] = fMultNorm(sum, FX_CHB2FX_DBL(pLnormii[i - 1]), &e); |
612 | 109k | x_sf[i] = sum_sf + pLnormii_sf[i - 1] + e + SUM_SAFETY; |
613 | 109k | } |
614 | 36.5k | } |
615 | | |
616 | | /** |
617 | | * \brief Solves L*x=b via backsubstitution according to the following |
618 | | * structure: |
619 | | * |
620 | | * x[3] = b[3]; |
621 | | * x[2] = b[2] - L[2][3]*x[3]; |
622 | | * x[1] = b[1] - L[1][2]*x[2] - L[1][3]*x[3]; |
623 | | * x[0] = b[0] - L[0][1]*x[1] - L[0][2]*x[2] - L[0][3]*x[3]; |
624 | | * |
625 | | * \param[in] numBands SBR crossover band index |
626 | | * \param[in] b the b in L*x=b (one-dimensional) |
627 | | * \param[out] x solution vector |
628 | | * \param[out] x_sf exponents of x[] |
629 | | */ |
630 | | static void backsubst_bw(const int numBands, const FIXP_DBL *const b, |
631 | 36.5k | FIXP_DBL *RESTRICT x, int *RESTRICT x_sf) { |
632 | 36.5k | int i, k; |
633 | 36.5k | int m; /* the trip counter that indexes incrementally through LnormInv1d[] */ |
634 | | |
635 | 36.5k | const FIXP_CHB *RESTRICT pLnormInv1d = |
636 | 36.5k | bsd[numBands - BSD_IDX_OFFSET].LnormInv1d; |
637 | 36.5k | const SCHAR *RESTRICT pLnormInv1d_sf = |
638 | 36.5k | bsd[numBands - BSD_IDX_OFFSET].LnormInv1d_sf; |
639 | | |
640 | 36.5k | x[POLY_ORDER] = b[POLY_ORDER]; |
641 | | |
642 | 146k | for (i = POLY_ORDER - 1, m = 0; i >= 0; i--) { |
643 | 109k | FIXP_DBL sum = b[i] >> SUM_SAFETY; |
644 | 109k | int sum_sf = x_sf[i]; /* sum's sf but disregarding SUM_SAFETY (added at the |
645 | | iteration's end) */ |
646 | | |
647 | 328k | for (k = i + 1; k <= POLY_ORDER; ++k, ++m) { |
648 | 219k | int e; |
649 | 219k | FIXP_DBL mult = fMultNorm(FX_CHB2FX_DBL(pLnormInv1d[m]), x[k], &e); |
650 | 219k | int mult_sf = pLnormInv1d_sf[m] + x_sf[k] + e; |
651 | | |
652 | | /* check if the new summand mult has a different sf than sum currently has |
653 | | */ |
654 | 219k | int diff = mult_sf - sum_sf; |
655 | | |
656 | 219k | if (diff > 0) { |
657 | | /* yes, and it requires the sum v to be adjusted (scaled down) */ |
658 | 107k | sum >>= diff; |
659 | 107k | sum_sf = mult_sf; |
660 | 111k | } else if (diff < 0) { |
661 | | /* yes, but here mult needs to be scaled down */ |
662 | 91.0k | mult >>= -diff; |
663 | 91.0k | } |
664 | | |
665 | | /* mult has now the same sf than what it is about to be added to. */ |
666 | | /* scale mult down additionally so that building the sum is overflow-safe. |
667 | | */ |
668 | 219k | sum -= (mult >> SUM_SAFETY); |
669 | 219k | } |
670 | | |
671 | 109k | x_sf[i] = sum_sf + SUM_SAFETY; |
672 | 109k | x[i] = sum; |
673 | 109k | } |
674 | 36.5k | } |
675 | | |
676 | | /** |
677 | | * \brief Solves a system of linear equations (L*x=b) with the Cholesky |
678 | | * algorithm. |
679 | | * |
680 | | * \param[in] numBands SBR crossover band index |
681 | | * \param[in,out] b input: vector b, output: solution vector p. |
682 | | * \param[in,out] b_sf input: exponent of b; output: exponent of solution |
683 | | * p. |
684 | | */ |
685 | | static void choleskySolve(const int numBands, FIXP_DBL *RESTRICT b, |
686 | 36.5k | int *RESTRICT b_sf) { |
687 | 36.5k | int i, e; |
688 | | |
689 | 36.5k | const FIXP_CHB *RESTRICT pBmul0 = bsd[numBands - BSD_IDX_OFFSET].Bmul0; |
690 | 36.5k | const SCHAR *RESTRICT pBmul0_sf = bsd[numBands - BSD_IDX_OFFSET].Bmul0_sf; |
691 | 36.5k | const FIXP_CHB *RESTRICT pBmul1 = bsd[numBands - BSD_IDX_OFFSET].Bmul1; |
692 | 36.5k | const SCHAR *RESTRICT pBmul1_sf = bsd[numBands - BSD_IDX_OFFSET].Bmul1_sf; |
693 | | |
694 | | /* normalize b */ |
695 | 36.5k | FIXP_DBL bnormed[POLY_ORDER + 1]; |
696 | 182k | for (i = 0; i <= POLY_ORDER; ++i) { |
697 | 146k | bnormed[i] = fMultNorm(b[i], FX_CHB2FX_DBL(pBmul0[i]), &e); |
698 | 146k | b_sf[i] += pBmul0_sf[i] + e; |
699 | 146k | } |
700 | | |
701 | 36.5k | backsubst_fw(numBands, bnormed, b, b_sf); |
702 | | |
703 | | /* normalize b again */ |
704 | 182k | for (i = 0; i <= POLY_ORDER; ++i) { |
705 | 146k | bnormed[i] = fMultNorm(b[i], FX_CHB2FX_DBL(pBmul1[i]), &e); |
706 | 146k | b_sf[i] += pBmul1_sf[i] + e; |
707 | 146k | } |
708 | | |
709 | 36.5k | backsubst_bw(numBands, bnormed, b, b_sf); |
710 | 36.5k | } |
711 | | |
712 | | /** |
713 | | * \brief Find polynomial approximation of vector y with implicit abscisas |
714 | | * x=0,1,2,3..n-1 |
715 | | * |
716 | | * The problem (V^T * V * p = V^T * y) is solved with Cholesky. |
717 | | * V is the Vandermode Matrix constructed with x = 0...n-1; |
718 | | * A = V^T * V; b = V^T * y; |
719 | | * |
720 | | * \param[in] numBands SBR crossover band index (BSD_IDX_OFFSET <= numBands <= |
721 | | * MAXLOWBANDS) |
722 | | * \param[in] y input vector (mantissa) |
723 | | * \param[in] y_sf exponents of y[] |
724 | | * \param[out] p output polynomial coefficients (mantissa) |
725 | | * \param[out] p_sf exponents of p[] |
726 | | */ |
727 | | static void polyfit(const int numBands, const FIXP_DBL *const y, const int y_sf, |
728 | 36.5k | FIXP_DBL *RESTRICT p, int *RESTRICT p_sf) { |
729 | 36.5k | int i, k; |
730 | 36.5k | LONG v[POLY_ORDER + 1]; |
731 | 36.5k | int sum_saftey = getLog2[numBands - 1]; |
732 | | |
733 | 36.5k | FDK_ASSERT((numBands >= BSD_IDX_OFFSET) && (numBands <= MAXLOWBANDS)); |
734 | | |
735 | | /* construct vector b[] temporarily stored in array p[] */ |
736 | 36.5k | FDKmemclear(p, (POLY_ORDER + 1) * sizeof(FIXP_DBL)); |
737 | | |
738 | | /* p[] are the sums over n values and each p[i] has its own sf */ |
739 | 182k | for (i = 0; i <= POLY_ORDER; ++i) p_sf[i] = 1 - DFRACT_BITS; |
740 | | |
741 | 630k | for (k = 0; k < numBands; k++) { |
742 | 593k | v[0] = (LONG)1; |
743 | 2.37M | for (i = 1; i <= POLY_ORDER; i++) { |
744 | 1.78M | v[i] = k * v[i - 1]; |
745 | 1.78M | } |
746 | | |
747 | 2.96M | for (i = 0; i <= POLY_ORDER; i++) { |
748 | 2.37M | if (v[POLY_ORDER - i] != 0 && y[k] != FIXP_DBL(0)) { |
749 | 1.65M | int e; |
750 | 1.65M | FIXP_DBL mult = fMultNorm((FIXP_DBL)v[POLY_ORDER - i], y[k], &e); |
751 | 1.65M | int sf = DFRACT_BITS - 1 + y_sf + e; |
752 | | |
753 | | /* check if the new summand has a different sf than the sum p[i] |
754 | | * currently has */ |
755 | 1.65M | int diff = sf - p_sf[i]; |
756 | | |
757 | 1.65M | if (diff > 0) { |
758 | | /* yes, and it requires the sum p[i] to be adjusted (scaled down) */ |
759 | 467k | p[i] >>= fMin(DFRACT_BITS - 1, diff); |
760 | 467k | p_sf[i] = sf; |
761 | 1.18M | } else if (diff < 0) { |
762 | | /* yes, but here mult needs to be scaled down */ |
763 | 717k | mult >>= -diff; |
764 | 717k | } |
765 | | |
766 | | /* mult has now the same sf than what it is about to be added to. |
767 | | scale mult down additionally so that building the sum is |
768 | | overflow-safe. */ |
769 | 1.65M | p[i] += mult >> sum_saftey; |
770 | 1.65M | } |
771 | 2.37M | } |
772 | 593k | } |
773 | | |
774 | 36.5k | p_sf[0] += sum_saftey; |
775 | 36.5k | p_sf[1] += sum_saftey; |
776 | 36.5k | p_sf[2] += sum_saftey; |
777 | 36.5k | p_sf[3] += sum_saftey; |
778 | | |
779 | 36.5k | choleskySolve(numBands, p, p_sf); |
780 | 36.5k | } |
781 | | |
782 | | /** |
783 | | * \brief Calculates the output of a POLY_ORDER-degree polynomial function |
784 | | * with Horner scheme: |
785 | | * |
786 | | * y(x) = p3 + p2*x + p1*x^2 + p0*x^3 |
787 | | * = p3 + x*(p2 + x*(p1 + x*p0)) |
788 | | * |
789 | | * The for loop iterates through the mult/add parts in y(x) as above, |
790 | | * during which regular upscaling ensures a stable exponent of the |
791 | | * result. |
792 | | * |
793 | | * \param[in] p coefficients as in y(x) |
794 | | * \param[in] p_sf exponents of p[] |
795 | | * \param[in] x_int non-fractional integer representation of x as in y(x) |
796 | | * \param[out] out_sf exponent of return value |
797 | | * |
798 | | * \return result y(x) |
799 | | */ |
800 | | static FIXP_DBL polyval(const FIXP_DBL *const p, const int *const p_sf, |
801 | 593k | const int x_int, int *out_sf) { |
802 | 593k | FDK_ASSERT(x_int <= 31); /* otherwise getLog2[] needs more elements */ |
803 | | |
804 | 593k | int k, x_sf; |
805 | 593k | int result_sf; /* working space to compute return value *out_sf */ |
806 | 593k | FIXP_DBL x; /* fractional value of x_int */ |
807 | 593k | FIXP_DBL result; /* return value */ |
808 | | |
809 | | /* if x == 0, then y(x) is just p3 */ |
810 | 593k | if (x_int != 0) { |
811 | 557k | x_sf = getLog2[x_int]; |
812 | 557k | x = (FIXP_DBL)x_int << (DFRACT_BITS - 1 - x_sf); |
813 | 557k | } else { |
814 | 36.5k | *out_sf = p_sf[3]; |
815 | 36.5k | return p[3]; |
816 | 36.5k | } |
817 | | |
818 | 557k | result = p[0]; |
819 | 557k | result_sf = p_sf[0]; |
820 | | |
821 | 2.22M | for (k = 1; k <= POLY_ORDER; ++k) { |
822 | 1.67M | FIXP_DBL mult = fMult(x, result); |
823 | 1.67M | int mult_sf = x_sf + result_sf; |
824 | | |
825 | 1.67M | int room = CountLeadingBits(mult); |
826 | 1.67M | mult <<= room; |
827 | 1.67M | mult_sf -= room; |
828 | | |
829 | 1.67M | FIXP_DBL pp = p[k]; |
830 | 1.67M | int pp_sf = p_sf[k]; |
831 | | |
832 | | /* equalize the shift factors of pp and mult so that we can sum them up */ |
833 | 1.67M | int diff = pp_sf - mult_sf; |
834 | | |
835 | 1.67M | if (diff > 0) { |
836 | 931k | diff = fMin(diff, DFRACT_BITS - 1); |
837 | 931k | mult >>= diff; |
838 | 931k | } else if (diff < 0) { |
839 | 463k | diff = fMax(diff, 1 - DFRACT_BITS); |
840 | 463k | pp >>= -diff; |
841 | 463k | } |
842 | | |
843 | | /* downshift by 1 to ensure safe summation */ |
844 | 1.67M | mult >>= 1; |
845 | 1.67M | mult_sf++; |
846 | 1.67M | pp >>= 1; |
847 | 1.67M | pp_sf++; |
848 | | |
849 | 1.67M | result_sf = fMax(pp_sf, mult_sf); |
850 | | |
851 | 1.67M | result = mult + pp; |
852 | | /* rarely, mult and pp happen to be almost equal except their sign, |
853 | | and then upon summation, result becomes so small, that it is within |
854 | | the inaccuracy range of a few bits, and then the relative error |
855 | | produced by this function may become HUGE */ |
856 | 1.67M | } |
857 | | |
858 | 557k | *out_sf = result_sf; |
859 | 557k | return result; |
860 | 593k | } |
861 | | |
862 | | void sbrDecoder_calculateGainVec(FIXP_DBL **sourceBufferReal, |
863 | | FIXP_DBL **sourceBufferImag, |
864 | | int sourceBuf_e_overlap, |
865 | | int sourceBuf_e_current, int overlap, |
866 | | FIXP_DBL *RESTRICT GainVec, int *GainVec_exp, |
867 | | int numBands, const int startSample, |
868 | 36.7k | const int stopSample) { |
869 | 36.7k | FIXP_DBL p[POLY_ORDER + 1]; |
870 | 36.7k | FIXP_DBL meanNrg; |
871 | 36.7k | FIXP_DBL LowEnv[MAXLOWBANDS]; |
872 | 36.7k | FIXP_DBL invNumBands = GetInvInt(numBands); |
873 | 36.7k | FIXP_DBL invNumSlots = GetInvInt(stopSample - startSample); |
874 | 36.7k | int i, loBand, exp, scale_nrg, scale_nrg_ov; |
875 | 36.7k | int sum_scale = 5, sum_scale_ov = 3; |
876 | | |
877 | 36.7k | if (overlap > 8) { |
878 | 15.5k | FDK_ASSERT(overlap <= 16); |
879 | 15.5k | sum_scale_ov += 1; |
880 | 15.5k | sum_scale += 1; |
881 | 15.5k | } |
882 | | |
883 | | /* exponents of energy values */ |
884 | 36.7k | sourceBuf_e_overlap = sourceBuf_e_overlap * 2 + sum_scale_ov; |
885 | 36.7k | sourceBuf_e_current = sourceBuf_e_current * 2 + sum_scale; |
886 | 36.7k | exp = fMax(sourceBuf_e_overlap, sourceBuf_e_current); |
887 | 36.7k | scale_nrg = sourceBuf_e_current - exp; |
888 | 36.7k | scale_nrg_ov = sourceBuf_e_overlap - exp; |
889 | | |
890 | 36.7k | meanNrg = (FIXP_DBL)0; |
891 | | /* Calculate the spectral envelope in dB over the current copy-up frame. */ |
892 | 631k | for (loBand = 0; loBand < numBands; loBand++) { |
893 | 594k | FIXP_DBL nrg_ov, nrg; |
894 | 594k | INT reserve = 0, exp_new; |
895 | 594k | FIXP_DBL maxVal = FL2FX_DBL(0.0f); |
896 | | |
897 | 25.9M | for (i = startSample; i < stopSample; i++) { |
898 | 25.3M | maxVal |= |
899 | 25.3M | (FIXP_DBL)((LONG)(sourceBufferReal[i][loBand]) ^ |
900 | 25.3M | ((LONG)sourceBufferReal[i][loBand] >> (DFRACT_BITS - 1))); |
901 | 25.3M | maxVal |= |
902 | 25.3M | (FIXP_DBL)((LONG)(sourceBufferImag[i][loBand]) ^ |
903 | 25.3M | ((LONG)sourceBufferImag[i][loBand] >> (DFRACT_BITS - 1))); |
904 | 25.3M | } |
905 | | |
906 | 594k | if (maxVal != FL2FX_DBL(0.0f)) { |
907 | 429k | reserve = CntLeadingZeros(maxVal) - 2; |
908 | 429k | } |
909 | | |
910 | 594k | nrg_ov = nrg = (FIXP_DBL)0; |
911 | 594k | if (scale_nrg_ov > -31) { |
912 | 4.47M | for (i = startSample; i < overlap; i++) { |
913 | 3.88M | nrg_ov += |
914 | 3.88M | (fPow2Div2(scaleValue(sourceBufferReal[i][loBand], reserve)) + |
915 | 3.88M | fPow2Div2(scaleValue(sourceBufferImag[i][loBand], reserve))) >> |
916 | 3.88M | sum_scale_ov; |
917 | 3.88M | } |
918 | 593k | } else { |
919 | 893 | scale_nrg_ov = 0; |
920 | 893 | } |
921 | 594k | if (scale_nrg > -31) { |
922 | 21.8M | for (i = overlap; i < stopSample; i++) { |
923 | 21.3M | nrg += (fPow2Div2(scaleValue(sourceBufferReal[i][loBand], reserve)) + |
924 | 21.3M | fPow2Div2(scaleValue(sourceBufferImag[i][loBand], reserve))) >> |
925 | 21.3M | sum_scale; |
926 | 21.3M | } |
927 | 591k | } else { |
928 | 2.69k | scale_nrg = 0; |
929 | 2.69k | } |
930 | | |
931 | 594k | nrg = (scaleValue(nrg_ov, scale_nrg_ov) >> 1) + |
932 | 594k | (scaleValue(nrg, scale_nrg) >> 1); |
933 | 594k | nrg = fMult(nrg, invNumSlots); |
934 | | |
935 | 594k | exp_new = |
936 | 594k | exp - (2 * reserve) + |
937 | 594k | 2; /* +1 for addition directly above, +1 for fPow2Div2 in loops above */ |
938 | | |
939 | | /* LowEnv = 10*log10(nrg) = log2(nrg) * 10/log2(10) */ |
940 | | /* exponent of logarithmic energy is 8 */ |
941 | 594k | if (nrg > (FIXP_DBL)0) { |
942 | 418k | int exp_log2; |
943 | 418k | nrg = CalcLog2(nrg, exp_new, &exp_log2); |
944 | 418k | nrg = scaleValue(nrg, exp_log2 - 6); |
945 | 418k | nrg = fMult(FL2FXCONST_SGL(LOG10FAC), nrg); |
946 | 418k | } else { |
947 | 176k | nrg = (FIXP_DBL)0; |
948 | 176k | } |
949 | 594k | LowEnv[loBand] = nrg; |
950 | 594k | meanNrg += fMult(nrg, invNumBands); |
951 | 594k | } |
952 | 36.7k | exp = 6 + 2; /* exponent of LowEnv: +2 is exponent of LOG10FAC */ |
953 | | |
954 | | /* subtract mean before polynomial approximation to reduce dynamic of p[] */ |
955 | 631k | for (loBand = 0; loBand < numBands; loBand++) { |
956 | 594k | LowEnv[loBand] = meanNrg - LowEnv[loBand]; |
957 | 594k | } |
958 | | |
959 | | /* For numBands < BSD_IDX_OFFSET (== POLY_ORDER+2) we dont get an |
960 | | overdetermined equation system. The calculated polynomial will exactly fit |
961 | | the input data and evaluating the polynomial will lead to the same vector |
962 | | than the original input vector: lowEnvSlope[] == lowEnv[] |
963 | | */ |
964 | 36.7k | if (numBands > POLY_ORDER + 1) { |
965 | | /* Find polynomial approximation of LowEnv */ |
966 | 36.5k | int p_sf[POLY_ORDER + 1]; |
967 | | |
968 | 36.5k | polyfit(numBands, LowEnv, exp, p, p_sf); |
969 | | |
970 | 630k | for (i = 0; i < numBands; i++) { |
971 | 593k | int sf; |
972 | | |
973 | | /* lowBandEnvSlope[i] = tmp; */ |
974 | 593k | FIXP_DBL tmp = polyval(p, p_sf, i, &sf); |
975 | | |
976 | | /* GainVec = 10^((mean(y)-y)/20) = 2^( (mean(y)-y) * log2(10)/20 ) */ |
977 | 593k | tmp = fMult(tmp, FL2FXCONST_SGL(LOG10FAC_INV)); |
978 | 593k | GainVec[i] = f2Pow(tmp, sf - 2, |
979 | 593k | &GainVec_exp[i]); /* -2 is exponent of LOG10FAC_INV */ |
980 | 593k | } |
981 | 36.5k | } else { /* numBands <= POLY_ORDER+1 */ |
982 | 825 | for (i = 0; i < numBands; i++) { |
983 | 660 | int sf = exp; /* exponent of LowEnv[] */ |
984 | | |
985 | | /* lowBandEnvSlope[i] = LowEnv[i]; */ |
986 | 660 | FIXP_DBL tmp = LowEnv[i]; |
987 | | |
988 | | /* GainVec = 10^((mean(y)-y)/20) = 2^( (mean(y)-y) * log2(10)/20 ) */ |
989 | 660 | tmp = fMult(tmp, FL2FXCONST_SGL(LOG10FAC_INV)); |
990 | 660 | GainVec[i] = f2Pow(tmp, sf - 2, |
991 | 660 | &GainVec_exp[i]); /* -2 is exponent of LOG10FAC_INV */ |
992 | 660 | } |
993 | 165 | } |
994 | 36.7k | } |