Coverage Report

Created: 2026-07-24 06:56

next uncovered line (L), next uncovered region (R), next uncovered branch (B)
/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
13.4M
#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
665k
#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
332k
#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
798k
#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
33.2k
                         FIXP_DBL *RESTRICT x, int *RESTRICT x_sf) {
569
33.2k
  int i, k;
570
33.2k
  int m; /* the trip counter that indexes incrementally through Lnorm1d[] */
571
572
33.2k
  const FIXP_CHB *RESTRICT pLnorm1d = bsd[numBands - BSD_IDX_OFFSET].Lnorm1d;
573
33.2k
  const SCHAR *RESTRICT pLnorm1d_sf = bsd[numBands - BSD_IDX_OFFSET].Lnorm1d_sf;
574
33.2k
  const FIXP_CHB *RESTRICT pLnormii = bsd[numBands - BSD_IDX_OFFSET].Lnormii;
575
33.2k
  const SCHAR *RESTRICT pLnormii_sf = bsd[numBands - BSD_IDX_OFFSET].Lnormii_sf;
576
577
33.2k
  x[0] = b[0];
578
579
133k
  for (i = 1, m = 0; i <= POLY_ORDER; ++i) {
580
99.8k
    FIXP_DBL sum = b[i] >> SUM_SAFETY;
581
99.8k
    int sum_sf = x_sf[i];
582
199k
    for (k = i - 1; k > 0; --k, ++m) {
583
99.8k
      int e;
584
99.8k
      FIXP_DBL mult = fMultNorm(FX_CHB2FX_DBL(pLnorm1d[m]), x[k], &e);
585
99.8k
      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
99.8k
      int diff = mult_sf - sum_sf;
590
591
99.8k
      if (diff > 0) {
592
        /* yes, and it requires the sum to be adjusted (scaled down) */
593
43.8k
        sum >>= diff;
594
43.8k
        sum_sf = mult_sf;
595
55.9k
      } else if (diff < 0) {
596
        /* yes, but here mult needs to be scaled down */
597
51.2k
        mult >>= -diff;
598
51.2k
      }
599
99.8k
      sum -= (mult >> SUM_SAFETY);
600
99.8k
    }
601
602
    /* - x[0] */
603
99.8k
    if (x_sf[0] > sum_sf) {
604
25.2k
      sum >>= (x_sf[0] - sum_sf);
605
25.2k
      sum_sf = x_sf[0];
606
25.2k
    }
607
99.8k
    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
99.8k
    int e;
611
99.8k
    x[i] = fMultNorm(sum, FX_CHB2FX_DBL(pLnormii[i - 1]), &e);
612
99.8k
    x_sf[i] = sum_sf + pLnormii_sf[i - 1] + e + SUM_SAFETY;
613
99.8k
  }
614
33.2k
}
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
33.2k
                         FIXP_DBL *RESTRICT x, int *RESTRICT x_sf) {
632
33.2k
  int i, k;
633
33.2k
  int m; /* the trip counter that indexes incrementally through LnormInv1d[] */
634
635
33.2k
  const FIXP_CHB *RESTRICT pLnormInv1d =
636
33.2k
      bsd[numBands - BSD_IDX_OFFSET].LnormInv1d;
637
33.2k
  const SCHAR *RESTRICT pLnormInv1d_sf =
638
33.2k
      bsd[numBands - BSD_IDX_OFFSET].LnormInv1d_sf;
639
640
33.2k
  x[POLY_ORDER] = b[POLY_ORDER];
641
642
133k
  for (i = POLY_ORDER - 1, m = 0; i >= 0; i--) {
643
99.8k
    FIXP_DBL sum = b[i] >> SUM_SAFETY;
644
99.8k
    int sum_sf = x_sf[i]; /* sum's sf but disregarding SUM_SAFETY (added at the
645
                             iteration's end) */
646
647
299k
    for (k = i + 1; k <= POLY_ORDER; ++k, ++m) {
648
199k
      int e;
649
199k
      FIXP_DBL mult = fMultNorm(FX_CHB2FX_DBL(pLnormInv1d[m]), x[k], &e);
650
199k
      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
199k
      int diff = mult_sf - sum_sf;
655
656
199k
      if (diff > 0) {
657
        /* yes, and it requires the sum v to be adjusted (scaled down) */
658
97.0k
        sum >>= diff;
659
97.0k
        sum_sf = mult_sf;
660
102k
      } else if (diff < 0) {
661
        /* yes, but here mult needs to be scaled down */
662
80.2k
        mult >>= -diff;
663
80.2k
      }
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
199k
      sum -= (mult >> SUM_SAFETY);
669
199k
    }
670
671
99.8k
    x_sf[i] = sum_sf + SUM_SAFETY;
672
99.8k
    x[i] = sum;
673
99.8k
  }
674
33.2k
}
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
33.2k
                          int *RESTRICT b_sf) {
687
33.2k
  int i, e;
688
689
33.2k
  const FIXP_CHB *RESTRICT pBmul0 = bsd[numBands - BSD_IDX_OFFSET].Bmul0;
690
33.2k
  const SCHAR *RESTRICT pBmul0_sf = bsd[numBands - BSD_IDX_OFFSET].Bmul0_sf;
691
33.2k
  const FIXP_CHB *RESTRICT pBmul1 = bsd[numBands - BSD_IDX_OFFSET].Bmul1;
692
33.2k
  const SCHAR *RESTRICT pBmul1_sf = bsd[numBands - BSD_IDX_OFFSET].Bmul1_sf;
693
694
  /* normalize b */
695
33.2k
  FIXP_DBL bnormed[POLY_ORDER + 1];
696
166k
  for (i = 0; i <= POLY_ORDER; ++i) {
697
133k
    bnormed[i] = fMultNorm(b[i], FX_CHB2FX_DBL(pBmul0[i]), &e);
698
133k
    b_sf[i] += pBmul0_sf[i] + e;
699
133k
  }
700
701
33.2k
  backsubst_fw(numBands, bnormed, b, b_sf);
702
703
  /* normalize b again */
704
166k
  for (i = 0; i <= POLY_ORDER; ++i) {
705
133k
    bnormed[i] = fMultNorm(b[i], FX_CHB2FX_DBL(pBmul1[i]), &e);
706
133k
    b_sf[i] += pBmul1_sf[i] + e;
707
133k
  }
708
709
33.2k
  backsubst_bw(numBands, bnormed, b, b_sf);
710
33.2k
}
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
33.2k
                    FIXP_DBL *RESTRICT p, int *RESTRICT p_sf) {
729
33.2k
  int i, k;
730
33.2k
  LONG v[POLY_ORDER + 1];
731
33.2k
  int sum_saftey = getLog2[numBands - 1];
732
733
33.2k
  FDK_ASSERT((numBands >= BSD_IDX_OFFSET) && (numBands <= MAXLOWBANDS));
734
735
  /* construct vector b[] temporarily stored in array p[] */
736
33.2k
  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
166k
  for (i = 0; i <= POLY_ORDER; ++i) p_sf[i] = 1 - DFRACT_BITS;
740
741
665k
  for (k = 0; k < numBands; k++) {
742
631k
    v[0] = (LONG)1;
743
2.52M
    for (i = 1; i <= POLY_ORDER; i++) {
744
1.89M
      v[i] = k * v[i - 1];
745
1.89M
    }
746
747
3.15M
    for (i = 0; i <= POLY_ORDER; i++) {
748
2.52M
      if (v[POLY_ORDER - i] != 0 && y[k] != FIXP_DBL(0)) {
749
1.75M
        int e;
750
1.75M
        FIXP_DBL mult = fMultNorm((FIXP_DBL)v[POLY_ORDER - i], y[k], &e);
751
1.75M
        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.75M
        int diff = sf - p_sf[i];
756
757
1.75M
        if (diff > 0) {
758
          /* yes, and it requires the sum p[i] to be adjusted (scaled down) */
759
453k
          p[i] >>= fMin(DFRACT_BITS - 1, diff);
760
453k
          p_sf[i] = sf;
761
1.30M
        } else if (diff < 0) {
762
          /* yes, but here mult needs to be scaled down */
763
802k
          mult >>= -diff;
764
802k
        }
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.75M
        p[i] += mult >> sum_saftey;
770
1.75M
      }
771
2.52M
    }
772
631k
  }
773
774
33.2k
  p_sf[0] += sum_saftey;
775
33.2k
  p_sf[1] += sum_saftey;
776
33.2k
  p_sf[2] += sum_saftey;
777
33.2k
  p_sf[3] += sum_saftey;
778
779
33.2k
  choleskySolve(numBands, p, p_sf);
780
33.2k
}
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
631k
                        const int x_int, int *out_sf) {
802
631k
  FDK_ASSERT(x_int <= 31); /* otherwise getLog2[] needs more elements */
803
804
631k
  int k, x_sf;
805
631k
  int result_sf;   /* working space to compute return value *out_sf */
806
631k
  FIXP_DBL x;      /* fractional value of x_int */
807
631k
  FIXP_DBL result; /* return value */
808
809
  /* if x == 0, then y(x) is just p3 */
810
631k
  if (x_int != 0) {
811
598k
    x_sf = getLog2[x_int];
812
598k
    x = (FIXP_DBL)x_int << (DFRACT_BITS - 1 - x_sf);
813
598k
  } else {
814
33.2k
    *out_sf = p_sf[3];
815
33.2k
    return p[3];
816
33.2k
  }
817
818
598k
  result = p[0];
819
598k
  result_sf = p_sf[0];
820
821
2.39M
  for (k = 1; k <= POLY_ORDER; ++k) {
822
1.79M
    FIXP_DBL mult = fMult(x, result);
823
1.79M
    int mult_sf = x_sf + result_sf;
824
825
1.79M
    int room = CountLeadingBits(mult);
826
1.79M
    mult <<= room;
827
1.79M
    mult_sf -= room;
828
829
1.79M
    FIXP_DBL pp = p[k];
830
1.79M
    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.79M
    int diff = pp_sf - mult_sf;
834
835
1.79M
    if (diff > 0) {
836
1.02M
      diff = fMin(diff, DFRACT_BITS - 1);
837
1.02M
      mult >>= diff;
838
1.02M
    } else if (diff < 0) {
839
491k
      diff = fMax(diff, 1 - DFRACT_BITS);
840
491k
      pp >>= -diff;
841
491k
    }
842
843
    /* downshift by 1 to ensure safe summation */
844
1.79M
    mult >>= 1;
845
1.79M
    mult_sf++;
846
1.79M
    pp >>= 1;
847
1.79M
    pp_sf++;
848
849
1.79M
    result_sf = fMax(pp_sf, mult_sf);
850
851
1.79M
    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.79M
  }
857
858
598k
  *out_sf = result_sf;
859
598k
  return result;
860
631k
}
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
33.4k
                                 const int stopSample) {
869
33.4k
  FIXP_DBL p[POLY_ORDER + 1];
870
33.4k
  FIXP_DBL meanNrg;
871
33.4k
  FIXP_DBL LowEnv[MAXLOWBANDS];
872
33.4k
  FIXP_DBL invNumBands = GetInvInt(numBands);
873
33.4k
  FIXP_DBL invNumSlots = GetInvInt(stopSample - startSample);
874
33.4k
  int i, loBand, exp, scale_nrg, scale_nrg_ov;
875
33.4k
  int sum_scale = 5, sum_scale_ov = 3;
876
877
33.4k
  if (overlap > 8) {
878
7.51k
    FDK_ASSERT(overlap <= 16);
879
7.51k
    sum_scale_ov += 1;
880
7.51k
    sum_scale += 1;
881
7.51k
  }
882
883
  /* exponents of energy values */
884
33.4k
  sourceBuf_e_overlap = sourceBuf_e_overlap * 2 + sum_scale_ov;
885
33.4k
  sourceBuf_e_current = sourceBuf_e_current * 2 + sum_scale;
886
33.4k
  exp = fMax(sourceBuf_e_overlap, sourceBuf_e_current);
887
33.4k
  scale_nrg = sourceBuf_e_current - exp;
888
33.4k
  scale_nrg_ov = sourceBuf_e_overlap - exp;
889
890
33.4k
  meanNrg = (FIXP_DBL)0;
891
  /* Calculate the spectral envelope in dB over the current copy-up frame. */
892
666k
  for (loBand = 0; loBand < numBands; loBand++) {
893
632k
    FIXP_DBL nrg_ov, nrg;
894
632k
    INT reserve = 0, exp_new;
895
632k
    FIXP_DBL maxVal = FL2FX_DBL(0.0f);
896
897
24.3M
    for (i = startSample; i < stopSample; i++) {
898
23.7M
      maxVal |=
899
23.7M
          (FIXP_DBL)((LONG)(sourceBufferReal[i][loBand]) ^
900
23.7M
                     ((LONG)sourceBufferReal[i][loBand] >> (DFRACT_BITS - 1)));
901
23.7M
      maxVal |=
902
23.7M
          (FIXP_DBL)((LONG)(sourceBufferImag[i][loBand]) ^
903
23.7M
                     ((LONG)sourceBufferImag[i][loBand] >> (DFRACT_BITS - 1)));
904
23.7M
    }
905
906
632k
    if (maxVal != FL2FX_DBL(0.0f)) {
907
451k
      reserve = CntLeadingZeros(maxVal) - 2;
908
451k
    }
909
910
632k
    nrg_ov = nrg = (FIXP_DBL)0;
911
632k
    if (scale_nrg_ov > -31) {
912
4.18M
      for (i = startSample; i < overlap; i++) {
913
3.55M
        nrg_ov +=
914
3.55M
            (fPow2Div2(scaleValue(sourceBufferReal[i][loBand], reserve)) +
915
3.55M
             fPow2Div2(scaleValue(sourceBufferImag[i][loBand], reserve))) >>
916
3.55M
            sum_scale_ov;
917
3.55M
      }
918
631k
    } else {
919
1.28k
      scale_nrg_ov = 0;
920
1.28k
    }
921
632k
    if (scale_nrg > -31) {
922
20.7M
      for (i = overlap; i < stopSample; i++) {
923
20.0M
        nrg += (fPow2Div2(scaleValue(sourceBufferReal[i][loBand], reserve)) +
924
20.0M
                fPow2Div2(scaleValue(sourceBufferImag[i][loBand], reserve))) >>
925
20.0M
               sum_scale;
926
20.0M
      }
927
629k
    } else {
928
2.92k
      scale_nrg = 0;
929
2.92k
    }
930
931
632k
    nrg = (scaleValue(nrg_ov, scale_nrg_ov) >> 1) +
932
632k
          (scaleValue(nrg, scale_nrg) >> 1);
933
632k
    nrg = fMult(nrg, invNumSlots);
934
935
632k
    exp_new =
936
632k
        exp - (2 * reserve) +
937
632k
        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
632k
    if (nrg > (FIXP_DBL)0) {
942
432k
      int exp_log2;
943
432k
      nrg = CalcLog2(nrg, exp_new, &exp_log2);
944
432k
      nrg = scaleValue(nrg, exp_log2 - 6);
945
432k
      nrg = fMult(FL2FXCONST_SGL(LOG10FAC), nrg);
946
432k
    } else {
947
200k
      nrg = (FIXP_DBL)0;
948
200k
    }
949
632k
    LowEnv[loBand] = nrg;
950
632k
    meanNrg += fMult(nrg, invNumBands);
951
632k
  }
952
33.4k
  exp = 6 + 2; /* exponent of LowEnv: +2 is exponent of LOG10FAC */
953
954
  /* subtract mean before polynomial approximation to reduce dynamic of p[] */
955
666k
  for (loBand = 0; loBand < numBands; loBand++) {
956
632k
    LowEnv[loBand] = meanNrg - LowEnv[loBand];
957
632k
  }
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
33.4k
  if (numBands > POLY_ORDER + 1) {
965
    /* Find polynomial approximation of LowEnv */
966
33.2k
    int p_sf[POLY_ORDER + 1];
967
968
33.2k
    polyfit(numBands, LowEnv, exp, p, p_sf);
969
970
665k
    for (i = 0; i < numBands; i++) {
971
631k
      int sf;
972
973
      /* lowBandEnvSlope[i] = tmp; */
974
631k
      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
631k
      tmp = fMult(tmp, FL2FXCONST_SGL(LOG10FAC_INV));
978
631k
      GainVec[i] = f2Pow(tmp, sf - 2,
979
631k
                         &GainVec_exp[i]); /* -2 is exponent of LOG10FAC_INV */
980
631k
    }
981
33.2k
  } else { /* numBands <= POLY_ORDER+1 */
982
865
    for (i = 0; i < numBands; i++) {
983
692
      int sf = exp; /* exponent of LowEnv[] */
984
985
      /* lowBandEnvSlope[i] = LowEnv[i]; */
986
692
      FIXP_DBL tmp = LowEnv[i];
987
988
      /* GainVec = 10^((mean(y)-y)/20) = 2^( (mean(y)-y) * log2(10)/20 ) */
989
692
      tmp = fMult(tmp, FL2FXCONST_SGL(LOG10FAC_INV));
990
692
      GainVec[i] = f2Pow(tmp, sf - 2,
991
692
                         &GainVec_exp[i]); /* -2 is exponent of LOG10FAC_INV */
992
692
    }
993
173
  }
994
33.4k
}