Line data Source code
1 : /******************************************************************************************************
2 :
3 : (C) 2022-2025 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB,
4 : Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,
5 : Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,
6 : Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other
7 : contributors to this repository. All Rights Reserved.
8 :
9 : This software is protected by copyright law and by international treaties.
10 : The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB,
11 : Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,
12 : Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,
13 : Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other
14 : contributors to this repository retain full ownership rights in their respective contributions in
15 : the software. This notice grants no license of any kind, including but not limited to patent
16 : license, nor is any license granted by implication, estoppel or otherwise.
17 :
18 : Contributors are required to enter into the IVAS codec Public Collaboration agreement before making
19 : contributions.
20 :
21 : This software is provided "AS IS", without any express or implied warranties. The software is in the
22 : development stage. It is intended exclusively for experts who have experience with such software and
23 : solely for the purpose of inspection. All implied warranties of non-infringement, merchantability
24 : and fitness for a particular purpose are hereby disclaimed and excluded.
25 :
26 : Any dispute, controversy or claim arising under or in relation to providing this software shall be
27 : submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in
28 : accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and
29 : the United Nations Convention on Contracts on the International Sales of Goods.
30 :
31 : *******************************************************************************************************/
32 :
33 : /*====================================================================================
34 : EVS Codec 3GPP TS26.443 Nov 04, 2021. Version 12.14.0 / 13.10.0 / 14.6.0 / 15.4.0 / 16.3.0
35 : ====================================================================================*/
36 :
37 : #include <stdint.h>
38 : #include "options.h"
39 : #ifdef DEBUGGING
40 : #include "debug.h"
41 : #endif
42 : #include <math.h>
43 : #include "cnst.h"
44 : #include "rom_com.h"
45 : #include "prot.h"
46 : #include "stl.h"
47 : #include "wmc_auto.h"
48 :
49 : /*-------------------------------------------------------------------*
50 : * Local definitions
51 : *-------------------------------------------------------------------*/
52 :
53 : /*! r: Approximate integer division for negative input */
54 1305876 : int16_t shrtCDivSignedApprox(
55 : const int16_t num, /* i : numerator */
56 : const int16_t den /* i : denominator */
57 : )
58 : {
59 : Word16 pool_part;
60 1305876 : pool_part = extract_h( L_mult( negate( abs_s( num ) ), lim_neg_inv_tbl_fx[den] ) );
61 : /* neg_in always, positive out always, so that positive truncation(rounding) towards zero is used */
62 1305876 : if ( num < 0 )
63 : {
64 26034 : pool_part = negate( pool_part ); /* make negative, one op */
65 : }
66 :
67 1305876 : return pool_part;
68 : }
69 :
70 : /*! r: Approximate integer division for positive input using lookup table */
71 2736979 : int32_t intLimCDivPos(
72 : const int32_t NUM, /* i : numerator */
73 : const int16_t DEN /* i : denominator */
74 : )
75 : {
76 : int32_t result;
77 :
78 2736979 : result = (int32_t) ( ( (uint64_t) ( NUM << 1 ) * intLimCDivInvDQ31[DEN] ) >> 32 );
79 :
80 2736979 : return result;
81 : }
82 :
83 : /*! r: Approximate integer division for signed input using lookup table */
84 1182091 : static int32_t intLimCDivSigned(
85 : const int32_t NUM, /* i : numerator */
86 : const int16_t DEN /* i : denominator */
87 : )
88 : {
89 : int32_t result;
90 : int32_t tmp_num;
91 :
92 1182091 : tmp_num = max( NUM, -NUM );
93 1182091 : result = (int32_t) sign( (float) NUM ) * (int32_t) ( ( (uint64_t) ( tmp_num << 1 ) * intLimCDivInvDQ31[DEN] ) >> 32 );
94 :
95 1182091 : return result;
96 : }
97 :
98 251434 : static void nearProjQ15(
99 : const int16_t x, /* i : input coefficient */
100 : int16_t *result /* o : projection */
101 : )
102 : {
103 251434 : int16_t NearProjC[4] = { 14967, -25518, 3415, 32351 };
104 : int32_t P12, P3;
105 :
106 251434 : P12 = ( (int32_t) x * (int32_t) NearProjC[0] >> 16 ) + (int32_t) NearProjC[1];
107 251434 : P3 = ( (int32_t) x * P12 >> 14 ) + (int32_t) NearProjC[2];
108 251434 : *result = (int16_t) ( ( (int16_t) x * P3 >> 15 ) + (int32_t) NearProjC[3] );
109 :
110 251434 : return;
111 : }
112 :
113 :
114 127182 : void obtainEnergyQuantizerDensity(
115 : const int16_t L_in, /* i : left vector energy */
116 : const int16_t R_in, /* i : right vector energy */
117 : int16_t *Density /* o : quantizer density */
118 : )
119 : {
120 : int32_t Rnrg;
121 : int16_t den;
122 : int16_t L, R;
123 :
124 127182 : L = L_in;
125 127182 : R = R_in;
126 127182 : den = ( 2 * L - 1 );
127 127182 : if ( den <= 67 )
128 : {
129 123753 : Rnrg = (int32_t) ( ( (uint64_t) ( (int32_t) R << 1 ) * intLimCDivInvDQ31[den] ) >> 32 );
130 : }
131 : else
132 : {
133 3429 : Rnrg = ( R / den );
134 : }
135 127182 : Rnrg += 28;
136 :
137 127182 : if ( ( R - 96 < 56 ) && ( R - 96 < Rnrg ) )
138 : {
139 740 : Rnrg = R - 96;
140 : }
141 126442 : else if ( 56 < Rnrg )
142 : {
143 776 : Rnrg = 56;
144 : }
145 :
146 127182 : Rnrg = max( 0, Rnrg ); /* _f table set to 1 for low Rnrg */
147 127182 : *Density = (int16_t) obtainEnergyQuantizerDensity_f[Rnrg];
148 :
149 127182 : return;
150 : }
151 :
152 :
153 1305876 : static void dsDirac2Dirac(
154 : const int16_t dsDiracIndex, /* i : input index */
155 : int16_t *diracs /* o : number of diracs */
156 : )
157 : {
158 : int16_t dsHighIndex;
159 :
160 1305876 : dsHighIndex = dsDiracIndex - DS_INDEX_LINEAR_END - 1;
161 1305876 : *diracs = dsDiracIndex;
162 1305876 : if ( dsDiracIndex > DS_INDEX_LINEAR_END )
163 : {
164 83384 : *diracs = dsHighDiracsTab[dsHighIndex];
165 : }
166 :
167 1305876 : return;
168 : }
169 :
170 :
171 1432898 : static void dsDiracPerQuanta(
172 : const int16_t td, /* i : Length of vector segment */
173 : const int16_t t_quanta, /* i : Assigned number of quanta */
174 : const int16_t dsm, /* i : Conservative rounding flag */
175 : const uint8_t *const *frQuanta, /* i : Quanta lookup table */
176 : int16_t *DsIdx /* o : Lookup table index */
177 : )
178 : {
179 : const uint8_t *sv;
180 : int16_t nsv;
181 : int16_t t_quanta_o;
182 : int16_t dsIndex;
183 : int16_t i;
184 :
185 1432898 : sv = frQuanta[td];
186 1432898 : nsv = sv[0];
187 :
188 1432898 : t_quanta_o = t_quanta - QUANTAQ3OFFSET;
189 1432898 : if ( t_quanta_o >= sv[nsv] )
190 : {
191 75852 : *DsIdx = nsv;
192 75852 : return;
193 : }
194 1357046 : if ( t_quanta_o <= sv[1] )
195 : {
196 182646 : *DsIdx = 1;
197 182646 : return;
198 : }
199 1174400 : dsIndex = 1 << frQuanta[0][td];
200 1174400 : if ( t_quanta_o > sv[nsv >> 1] )
201 : {
202 384553 : dsIndex = nsv - dsIndex; /*single op*/
203 : }
204 4189252 : for ( i = frQuanta[0][td] - 1; i >= 0; i-- )
205 : {
206 3014852 : dsIndex += shl( sub( shl( lshr( sub( sv[dsIndex], t_quanta_o ), 15 ), 1 ), 1 ), i );
207 : }
208 1174400 : dsIndex += ( t_quanta_o > sv[dsIndex] );
209 1174400 : dsIndex -= ( dsIndex > 1 );
210 :
211 1174400 : if ( dsm == PVQ_CONS )
212 : {
213 462062 : *DsIdx = dsIndex;
214 462062 : return;
215 : }
216 712338 : *DsIdx = dsIndex + ( ( sv[dsIndex + 1] + sv[dsIndex] ) < ( t_quanta_o << 1 ) );
217 :
218 712338 : return;
219 : }
220 :
221 :
222 132124 : void QuantaPerDsDirac(
223 : const int16_t td, /* i : Length of vector segment */
224 : const int16_t dsDiracIndex, /* i : Quanta table index */
225 : const uint8_t *const *dimFrQuanta, /* i : Quanta lookup table */
226 : int16_t *Quanta /* i : Quanta */
227 : )
228 : {
229 132124 : if ( !dsDiracIndex )
230 : {
231 0 : *Quanta = 0;
232 0 : return;
233 : }
234 132124 : *Quanta = ( (int16_t) dimFrQuanta[td][dsDiracIndex] ) + QUANTAQ3OFFSET;
235 :
236 132124 : return;
237 : }
238 :
239 :
240 1305876 : void conservativeL1Norm(
241 : const int16_t L, /* i : Length of vector segment */
242 : const int16_t Qvec, /* i : Assigned number of quanta */
243 : const int16_t Fcons, /* i : Conservative rounding flag */
244 : const int16_t Qavail, /* i : Input quanta remaining */
245 : const int16_t Qreserv, /* i : Input quanta in reservoir */
246 : const int16_t Dspec, /* i : assigned diracs from bitalloc */
247 : int16_t *Dvec, /* o : actual number of diracs */
248 : int16_t *Qspare, /* o : Output quanta remaining */
249 : int16_t *Qreservplus, /* o : Output quanta in reservoir */
250 : int16_t *Dspecplus /* o : Output number of diracs */
251 : )
252 : {
253 : int16_t Minit, Mprime;
254 : int16_t Qtestminus;
255 : const uint8_t *frQuantaL;
256 :
257 1305876 : frQuantaL = hBitsN[L];
258 :
259 1305876 : *Qreservplus = Qreserv + Qvec - QUANTAQ3OFFSET;
260 :
261 1305876 : dsDiracPerQuanta( L, Qvec, Fcons, hBitsN, &Minit );
262 1305876 : Mprime = Minit;
263 : do
264 : {
265 1438464 : Qtestminus = frQuantaL[Mprime];
266 1438464 : *Qspare = Qavail - Qtestminus;
267 1438464 : Mprime = Mprime - 1;
268 1438464 : } while ( Mprime >= 0 && *Qspare < QUANTAQ3OFFSET );
269 :
270 1305876 : if ( Mprime < 0 )
271 : {
272 120626 : *Qspare = Qavail + QUANTAQ3OFFSET;
273 : }
274 :
275 1305876 : dsDirac2Dirac( Mprime + 1, Dvec );
276 :
277 1305876 : *Dspecplus = Dspec + *Dvec;
278 :
279 1305876 : *Qreservplus -= frQuantaL[Minit];
280 :
281 1305876 : *Qspare -= QUANTAQ3OFFSET;
282 :
283 1305876 : return;
284 : }
285 :
286 1178694 : void bandBitsAdjustment(
287 : const int16_t Brc, /* i : Current number of read quanta in range coder */
288 : const uint32_t INTrc, /* i : Range coder state */
289 : const int16_t Bavail, /* i : Available number of quanta */
290 : const int16_t Nbands, /* i : Number of bands */
291 : const int16_t D, /* i : Remaining number of bands to encode */
292 : const int16_t L, /* i : Size of current band */
293 : const int16_t Bband, /* i : Quanta allocation for current band */
294 : const int16_t Breserv, /* i : Quanta reservoir */
295 : int16_t *Bband_adj, /* o : Actual used number of quanta */
296 : int16_t *Brem, /* o : Quanta remaining */
297 : int16_t *Breservplus /* o : Quanta pool size */
298 : )
299 : {
300 : int16_t Btemp;
301 : int16_t Bff;
302 : int16_t Dtmp, tmp_short1, tmp_short2;
303 :
304 1178694 : rangeCoderFinalizationFBits( Brc, INTrc, &Bff );
305 :
306 1178694 : if ( D < Nbands )
307 : {
308 1055069 : Dtmp = min( D, 3 ); /* avoid branch cost */
309 1055069 : Btemp = (int16_t) intLimCDivSigned( Breserv - Bff, Dtmp ); /* result always fits in Word16 */
310 :
311 1055069 : *Breservplus = Bband + Breserv;
312 : }
313 : else
314 : {
315 123625 : Btemp = 0;
316 123625 : *Breservplus = Bband + Bff;
317 : }
318 :
319 1178694 : *Bband_adj = Bband;
320 1178694 : tmp_short1 = L * 80;
321 1178694 : if ( tmp_short1 < Bband )
322 : {
323 52 : *Bband_adj = tmp_short1;
324 : }
325 :
326 1178694 : tmp_short1 = Bavail - Bff;
327 1178694 : tmp_short2 = *Bband_adj + Btemp;
328 :
329 1178694 : *Bband_adj = tmp_short2;
330 1178694 : if ( tmp_short1 < tmp_short2 )
331 : {
332 151937 : *Bband_adj = tmp_short1;
333 : }
334 1178694 : if ( *Bband_adj < 0 )
335 : {
336 263 : *Bband_adj = 0;
337 : }
338 1178694 : *Brem = Bavail - Bff;
339 :
340 1178694 : return;
341 : }
342 :
343 :
344 : /*-------------------------------------------------------------------*
345 : * local_norm_s()
346 : *
347 : *
348 : *-------------------------------------------------------------------*/
349 :
350 : /*! r: n shifts needed to normalize */
351 251434 : static int16_t local_norm_s(
352 : int16_t short_var /* i/o: signed 16 bit variable / normalized value */
353 : )
354 : {
355 : int16_t short_res;
356 :
357 251434 : if ( short_var == -1 )
358 : {
359 0 : return ( 16 - 1 );
360 : }
361 : else
362 : {
363 251434 : if ( short_var == 0 )
364 : {
365 0 : return 0;
366 : }
367 :
368 : else
369 : {
370 251434 : if ( short_var < 0 )
371 : {
372 0 : short_var = ~short_var;
373 : }
374 :
375 343663 : for ( short_res = 0; short_var < 0x4000; short_res++ )
376 : {
377 92229 : short_var <<= 1;
378 : }
379 : }
380 : }
381 :
382 251434 : return ( short_res );
383 : }
384 :
385 : /*! r: ratio Q11 */
386 125717 : static int16_t Ratio_base2Q11(
387 : const int16_t opp, /* i : opposite Q15 */
388 : const int16_t near /* i : near Q15 */
389 : )
390 : {
391 : Word16 mc, nc, ms, ns, d, z;
392 : Word16 result;
393 : Word32 acc;
394 :
395 125717 : ns = local_norm_s( opp ); /* exponent*/
396 125717 : nc = local_norm_s( near ); /* exponent */
397 :
398 125717 : ms = opp << ns; /* mantissa */
399 125717 : mc = near << nc; /* mantissa */
400 :
401 125717 : acc = L_mac( 538500224L, mc, -2776 ); /* a0*mc + a1, acc(Q27), a0(Q11), a1(Q27)*/
402 125717 : z = mac_r( acc, ms, -2776 ); /* z in Q11, a0 in Q11 */
403 125717 : d = sub( ms, mc ); /* d in Q15 */
404 125717 : z = mult_r( z, d ); /* z in Q11 */
405 :
406 125717 : result = add( z, shl( sub( nc, ns ), 11 ) );
407 :
408 125717 : return result;
409 : }
410 :
411 :
412 125717 : static void Ratio_rQ3(
413 : int16_t opp, /* i : opposite */
414 : int16_t near, /* i : near */
415 : int16_t *result /* o : ratio */
416 : )
417 : {
418 125717 : int16_t tmp = 128;
419 :
420 125717 : tmp += Ratio_base2Q11( opp, near );
421 125717 : *result = ( tmp >> 8 );
422 :
423 125717 : return;
424 : }
425 :
426 :
427 127182 : void densityAngle2RmsProjDec(
428 : const int16_t D, /* i : density */
429 : const int16_t indexphi, /* i : decoded index from AR dec */
430 : int16_t *oppQ15, /* o : opposite */
431 : int16_t *nearQ15, /* o : near */
432 : int16_t *oppRatioQ3 /* o : ratio */
433 : )
434 : {
435 : int16_t phiQ14q;
436 : int32_t L_2xPhiQ14q;
437 : int16_t oppTail, nearTail;
438 :
439 127182 : phiQ14q = (int16_t) intLimCDivPos( labs( indexphi ) << 13, D >> 1 );
440 127182 : if ( ( 0xFFFe & D ) == 0 ) /* even -> positive, odd -> 0 */
441 : {
442 712 : phiQ14q = 1 << 13; /* one op */
443 : }
444 :
445 : #define C_TOP ( ( 1 << 14 ) - ( 1 << 6 ) )
446 : #define C_BOT ( 1 << 6 )
447 : #define C_SHRT_MAXABS ( 1L << 15 )
448 : #define C_SHRT_MAX_POS ( ( 1 << 15 ) - 1 )
449 : #define C_Q14_MAX ( 1 << 14 )
450 :
451 127182 : oppTail = phiQ14q >= C_TOP;
452 127182 : nearTail = phiQ14q < C_BOT;
453 127182 : if ( !( oppTail || nearTail ) )
454 : {
455 125717 : L_2xPhiQ14q = 2 * (int32_t) phiQ14q;
456 125717 : nearProjQ15( (int16_t) ( C_SHRT_MAXABS - L_2xPhiQ14q ), oppQ15 );
457 125717 : nearProjQ15( phiQ14q << 1, nearQ15 );
458 125717 : Ratio_rQ3( *oppQ15, *nearQ15, oppRatioQ3 );
459 : }
460 : else
461 : {
462 1465 : *oppRatioQ3 = ( 1 - 2 * nearTail ) * C_Q14_MAX;
463 1465 : *oppQ15 = oppTail * C_SHRT_MAX_POS;
464 1465 : *nearQ15 = nearTail * C_SHRT_MAX_POS;
465 : }
466 :
467 : #undef C_TOP
468 : #undef C_BOT
469 : #undef C_SHRT_MAXABS
470 : #undef C_SHRT_MAX_POS
471 : #undef C_Q14_MAX
472 :
473 127182 : return;
474 : }
475 :
476 :
477 32073 : void densityAngle2RmsProjEnc(
478 : const int16_t D, /* i : density */
479 : const int16_t phiQ14uq, /* i : angle */
480 : int16_t *indexphi, /* o : index */
481 : int16_t *oppQ15, /* o : opposite */
482 : int16_t *nearQ15, /* o : near */
483 : int16_t *oppRatioQ3 /* o : ratio */
484 : )
485 : {
486 : #define C_MAX13 ( 1L << 13 )
487 32073 : *indexphi = ( D * phiQ14uq + C_MAX13 ) >> 14;
488 : #undef C_MAX13
489 :
490 32073 : if ( ( 0xFFFE & D ) == 0 )
491 : {
492 181 : *indexphi = -1; /* one op */
493 : }
494 32073 : densityAngle2RmsProjDec( D, *indexphi, oppQ15, nearQ15, oppRatioQ3 );
495 :
496 32073 : return;
497 : }
498 :
499 127182 : void NearOppSplitAdjustment(
500 : const int16_t qband, /* i : quanta for current band */
501 : const int16_t qzero, /* i : range coder finalization quanta */
502 : const int16_t Qac, /* i : range coder current quanta */
503 : const uint32_t INTac, /* i : range coder state */
504 : const int16_t qglobal, /* i : quanta input */
505 : const int16_t FlagCons, /* i : conservative rounding flag */
506 : const int16_t Np, /* i : number of parts */
507 : const int16_t Nhead, /* i : first part */
508 : const int16_t Ntail, /* i : remaining parts */
509 : const int16_t Nnear, /* i : length of near component */
510 : const int16_t Nopp, /* i : length of opposite component */
511 : int16_t oppRQ3, /* i : ratio */
512 : int16_t *qnear, /* o : quantized near */
513 : int16_t *qopp, /* o : quantized opposite */
514 : int16_t *qglobalupd /* o : quanta remaining */
515 : )
516 : {
517 : int16_t qac, qskew, qboth, QIa, QIb;
518 : int32_t L_QIb, qnum;
519 : int16_t qmin, qavg, Midx;
520 :
521 127182 : rangeCoderFinalizationFBits( Qac, INTac, &qac );
522 127182 : qboth = qband - ( qac - qzero );
523 127182 : qskew = 0; /* skew calc code */
524 127182 : if ( Nhead > 1 )
525 : {
526 127022 : qavg = (int16_t) intLimCDivSigned( qboth, Np );
527 127022 : dsDiracPerQuanta( Ntail, qavg, FlagCons, hBitsN, &Midx );
528 127022 : QuantaPerDsDirac( Nhead, Midx, hBitsN, &qmin );
529 :
530 127022 : qskew = qavg - qmin;
531 127022 : qskew = max( 0, qskew );
532 : } /* end of skew calc code*/
533 :
534 127182 : QIa = (int16_t) intLimCDivPos( Nopp, Nnear ) + 1; /* always positive Word16 out */
535 127182 : qnum = qband + qzero - qac - qskew - Nopp * (int32_t) oppRQ3;
536 127182 : L_QIb = 0;
537 127182 : if ( qnum > 0 )
538 : {
539 125154 : L_QIb = intLimCDivPos( qnum, QIa );
540 : }
541 127182 : QIb = (int16_t) min( 32767L, L_QIb ); /* saturate, L_QIb >= 0 */
542 127182 : *qnear = QIb;
543 127182 : if ( QIb > qboth )
544 : {
545 1743 : *qnear = qboth; /*single op*/
546 : }
547 127182 : *qopp = qboth - *qnear;
548 127182 : *qglobalupd = qglobal - ( qac - qzero );
549 :
550 127182 : return;
551 : }
552 :
553 : /*--------------------------------------------------------------------------*
554 : * apply_gain()
555 : *
556 : * Apply gain
557 : *--------------------------------------------------------------------------*/
558 :
559 122326 : void apply_gain(
560 : const int16_t *ord, /* i : Indices for energy order */
561 : const int16_t *band_start, /* i : Sub band start indices */
562 : const int16_t *band_end, /* i : Sub band end indices */
563 : const int16_t num_sfm, /* i : Number of bands */
564 : const float *gains, /* i : Band gain vector */
565 : float *xq /* i/o: float synthesis / gain adjusted synth */
566 : )
567 : {
568 : int16_t band, i;
569 : float g;
570 :
571 1903013 : for ( band = 0; band < num_sfm; band++ )
572 : {
573 1780687 : g = gains[ord[band]];
574 :
575 32738127 : for ( i = band_start[band]; i < band_end[band]; i++ )
576 : {
577 30957440 : xq[i] *= g;
578 : }
579 : }
580 :
581 122326 : return;
582 : }
583 :
584 :
585 : /*--------------------------------------------------------------------------*
586 : * fine_gain_quant()
587 : *
588 : * Fine gain quantization
589 : *--------------------------------------------------------------------------*/
590 :
591 29902 : void fine_gain_quant(
592 : BSTR_ENC_HANDLE hBstr, /* i/o: encoder bitstream handle */
593 : const int16_t *ord, /* i : Indices for energy order */
594 : const int16_t num_sfm, /* i : Number of bands */
595 : const int16_t *gain_bits, /* i : Gain adjustment bits per sub band */
596 : float *fg_pred, /* i/o: Predicted gains / Corrected gains */
597 : const float *gopt /* i : Optimal gains */
598 : )
599 : {
600 : int16_t band;
601 : int16_t gbits;
602 : int16_t idx;
603 : float gain_db, gain_dbq;
604 : float err;
605 :
606 479087 : for ( band = 0; band < num_sfm; band++ )
607 : {
608 449185 : gbits = gain_bits[ord[band]];
609 449185 : if ( fg_pred[band] != 0 && gbits > 0 )
610 : {
611 36121 : err = gopt[band] / fg_pred[band];
612 36121 : gain_db = 20.0f * (float) log10( err );
613 36121 : idx = (int16_t) squant( gain_db, &gain_dbq, finegain[gbits - 1], gain_cb_size[gbits - 1] );
614 36121 : push_indice( hBstr, IND_PVQ_FINE_GAIN, idx, gbits );
615 :
616 : /* Update prediced gain with quantized correction */
617 36121 : fg_pred[band] *= (float) pow( 10, gain_dbq * 0.05f );
618 : }
619 : }
620 :
621 29902 : return;
622 : }
623 :
624 : /*-------------------------------------------------------------------*
625 : * srt_vec_ind()
626 : *
627 : * sort vector and save sorting indices
628 : *-------------------------------------------------------------------*/
629 :
630 1269219 : void srt_vec_ind(
631 : const int16_t *linear, /* i : linear input */
632 : int16_t *srt, /* o : sorted output */
633 : int16_t *I, /* o : index for sorted output */
634 : const int16_t length /* i : length of vector */
635 : )
636 : {
637 : float linear_f[MAX_SRT_LEN];
638 : float srt_f[MAX_SRT_LEN];
639 :
640 1269219 : mvs2r( linear, linear_f, length );
641 1269219 : srt_vec_ind_f( linear_f, srt_f, I, length );
642 1269219 : mvr2s( srt_f, srt, length );
643 :
644 1269219 : return;
645 : }
646 :
647 : /*-------------------------------------------------------------------*
648 : * srt_vec_ind_f()
649 : *
650 : * sort vector and save sorting indices, using float input values
651 : *-------------------------------------------------------------------*/
652 :
653 1269219 : void srt_vec_ind_f(
654 : const float *linear, /* i : linear input */
655 : float *srt, /* o : sorted output */
656 : int16_t *I, /* o : index for sorted output */
657 : const int16_t length /* i : length of vector */
658 : )
659 : {
660 : int16_t pos, npos;
661 : int16_t idxMem;
662 : float valMem;
663 :
664 : /*initialize */
665 3225246 : for ( pos = 0; pos < length; pos++ )
666 : {
667 1956027 : I[pos] = pos;
668 : }
669 :
670 1269219 : mvr2r( linear, srt, length );
671 :
672 : /* now iterate */
673 1956027 : for ( pos = 0; pos < ( length - 1 ); pos++ )
674 : {
675 3030092 : for ( npos = ( pos + 1 ); npos < length; npos++ )
676 : {
677 2343284 : if ( srt[npos] < srt[pos] )
678 : {
679 1073135 : idxMem = I[pos];
680 1073135 : I[pos] = I[npos];
681 1073135 : I[npos] = idxMem;
682 :
683 1073135 : valMem = srt[pos];
684 1073135 : srt[pos] = srt[npos];
685 1073135 : srt[npos] = valMem;
686 : }
687 : }
688 : }
689 :
690 1269219 : return;
691 : }
692 :
693 : #define WMC_TOOL_SKIP
694 :
695 : /*-------------------------------------------------------------------*
696 : * UMult_32_32()
697 : *
698 : *
699 : *-------------------------------------------------------------------*/
700 :
701 : /*! r: Multiplication result */
702 14427072 : uint32_t UMult_32_32(
703 : const uint32_t UL_var1, /* i : factor 1 */
704 : const uint32_t UL_var2 /* i : factor 2 */
705 : )
706 : {
707 : uint64_t tmp;
708 :
709 14427072 : tmp = (uint64_t) UL_var1 * (uint64_t) UL_var2;
710 :
711 14427072 : return (uint32_t) ( tmp >> 32 );
712 : }
713 :
714 : /*-------------------------------------------------------------------*
715 : * UL_div()
716 : *
717 : * Calculate UL_num/UL_den. UL_num assumed to be Q31, UL_den assumed
718 : * to be Q32, then result is in Q32.
719 : *-------------------------------------------------------------------*/
720 :
721 : /*! r: Division result */
722 1311552 : static uint32_t UL_div(
723 : const uint32_t UL_num, /* i : numerator */
724 : const uint32_t UL_den /* i : denominator */
725 : )
726 : {
727 : uint32_t UL_e, UL_Q;
728 : uint32_t UL_msb;
729 : int16_t i;
730 :
731 1311552 : UL_e = 0xffffffff - UL_den;
732 1311552 : UL_Q = UL_num;
733 :
734 7869312 : for ( i = 0; i < 5; i++ )
735 : {
736 6557760 : UL_msb = UMult_32_32( UL_Q, UL_e );
737 6557760 : UL_Q = UL_Q + UL_msb;
738 6557760 : UL_e = UMult_32_32( UL_e, UL_e );
739 : }
740 :
741 1311552 : return UL_Q;
742 : }
743 :
744 :
745 : /*-------------------------------------------------------------------*
746 : * UL_inverse()
747 : *
748 : * Calculate inverse of UL_val. Output in Q_exp.
749 : *-------------------------------------------------------------------*/
750 :
751 : /*! r: inverse */
752 1311552 : uint32_t UL_inverse(
753 : const uint32_t UL_val, /* i : input value Q_exp */
754 : int16_t *exp /* i/o: input exp / result exp */
755 : )
756 : {
757 : uint32_t UL_tmp;
758 :
759 1311552 : *exp = norm_ul( UL_val ); /* aligned to BASOP */
760 1311552 : UL_tmp = UL_val << ( *exp ); /* Q32*/
761 :
762 1311552 : *exp = 32 + 31 - *exp;
763 :
764 1311552 : return UL_div( 0x80000000, UL_tmp );
765 : }
766 :
767 : /*-----------------------------------------------------------------------------
768 : * ratio()
769 : *
770 : * Divide the numerator by the denominator.
771 : *----------------------------------------------------------------------------*/
772 :
773 : /*! r: ratio */
774 126470 : Word16 ratio(
775 : const Word32 numer, /* i : numerator */
776 : const Word32 denom, /* i : denominator */
777 : Word16 *expo /* i/o: input exp / result exp */
778 : )
779 : {
780 : Word16 expNumer, expDenom;
781 : Word16 manNumer, manDenom;
782 : Word16 quotient;
783 :
784 126470 : expDenom = norm_l( denom ); /* exponent */
785 126470 : manDenom = extract_h( L_shl( denom, expDenom ) ); /* mantissa */
786 126470 : expNumer = norm_l( numer ); /* exponent */
787 126470 : manNumer = extract_h( L_shl( numer, expNumer ) ); /* mantissa */
788 126470 : manNumer = shr( manNumer, 1 ); /* Ensure the numerator < the denominator */
789 126470 : quotient = div_s( manNumer, manDenom ); /* in Q14 */
790 :
791 126470 : *expo = sub( expNumer, expDenom );
792 :
793 126470 : return quotient; /* Q14 */
794 : }
795 :
796 : /*-----------------------------------------------------------------------------
797 : * atan2_fx():
798 : *
799 : * Approximates arctan piecewise with various 4th to 5th order least square fit
800 : * polynomials for input in 5 segments:
801 : * - 0.0 to 1.0
802 : * - 1.0 to 2.0
803 : * - 2.0 to 4.0
804 : * - 4.0 to 8.0
805 : * - 8.0 to infinity
806 : *---------------------------------------------------------------------------*/
807 :
808 : /*! r: Angle between 0 and EVS_PI/2 radian (Q14) */
809 126470 : Word16 atan2_fx(
810 : const Word32 y, /* i : near side (Argument must be positive) (Q15) */
811 : const Word32 x /* i : opposite side (Q15) */
812 : )
813 : {
814 : Word32 acc, arg;
815 : Word16 man, expo, reciprocal;
816 : Word16 angle, w, z;
817 :
818 126470 : IF( L_sub( x, 0 ) == 0 )
819 : {
820 0 : return 25736; /* EVS_PI/2 in Q14 */
821 : }
822 126470 : man = ratio( y, x, &expo ); /* man in Q14 */
823 126470 : expo = sub( expo, ( 15 - 14 ) ); /* Now, man is considered in Q15 */
824 126470 : arg = L_shr( (Word32) man, expo );
825 :
826 126470 : IF( L_shr( arg, 3 + 15 ) != 0 )
827 : /*===============================*
828 : * 8.0 <= x < infinity *
829 : *===============================*/
830 : {
831 : /* atan(x) = EVS_PI/2 - 1/x + 1/(3x^3) - 1/(5x^5) + ...
832 : * ~ EVS_PI/2 - 1/x, for x >= 8.
833 : */
834 0 : expo = norm_l( arg );
835 0 : man = extract_h( L_shl( arg, expo ) );
836 0 : reciprocal = div_s( 0x3fff, man );
837 0 : expo = sub( 15 + 1, expo );
838 0 : reciprocal = shr( reciprocal, expo ); /* Q14*/
839 0 : angle = sub( 25736, reciprocal ); /* Q14 (EVS_PI/2 - 1/x) */
840 :
841 : /* For 8.0 <= x < 10.0, 1/(5x^5) is not completely negligible.
842 : * For more accurate result, add very small correction term.
843 : */
844 0 : IF( L_sub( L_shr( arg, 15 ), 10L ) < 0 )
845 : {
846 0 : angle = add( angle, 8 ); /* Add tiny correction term. */
847 : }
848 : }
849 126470 : ELSE IF( L_shr( arg, 2 + 15 ) != 0 )
850 : /*==========================*
851 : * 4.0 <= x < 8.0 *
852 : *==========================*/
853 : {
854 : /* interval: [3.999, 8.001]
855 : * atan(x) ~ (((a0*x + a1)*x + a2)*x + a3)*x + a4
856 : * = (((a0*8*y + a1)*8*y + a2)*8*y + a3)*8*y + a4 Substitute 8*y -> x
857 : * = (((a0*8^3*y + a1*8^2)*y + a2*8)*y + a3)*8*y + a4
858 : * = ((( c0*y + c1)*y + c2)*y + c3)*8*y + c4,
859 : * where y = x/8
860 : * and a0 = -1.28820869667651e-04, a1 = 3.88263533346295e-03,
861 : * a2 = -4.64216306484597e-02, a3 = 2.75986060068931e-01,
862 : * a4 = 7.49208077809799e-01.
863 : */
864 0 : w = extract_l( L_shr( arg, 3 ) ); /* Q15 y = x/8 */
865 0 : acc = 533625337L; /* Q31 c1 = a1*8^2 */
866 0 : move32();
867 0 : z = mac_r( acc, w, -2161 ); /* Q15 c0 = a0*8^3 */
868 0 : acc = -797517542L; /* Q31 c2 = a2*8 */
869 0 : move32();
870 0 : z = mac_r( acc, w, z ); /* Q15 */
871 0 : acc = 592675551L; /* Q31 c3 = a3 */
872 0 : move32();
873 0 : z = mac_r( acc, w, z ); /* z (in:Q15, out:Q12) */
874 0 : acc = 201114012L; /* Q28 c4 = a4 */
875 0 : move32();
876 0 : acc = L_mac( acc, w, z ); /* Q28 */
877 0 : angle = extract_l( L_shr( acc, ( 28 - 14 ) ) ); /* Q14 result of atan(x), where 4 <= x < 8 */
878 : }
879 126470 : ELSE IF( L_shr( arg, 1 + 15 ) != 0 )
880 : /*==========================*
881 : * 2.0 <= x < 4.0 *
882 : *==========================*/
883 : {
884 : /* interval: [1.999, 4.001]
885 : * atan(x) ~ (((a0*x + a1)*x + a2)*x + a3)*x + a4
886 : * = (((a0*4*y + a1)*4*y + a2)*4*y + a3)*4*y + a4 Substitute 4*y -> x
887 : * = (((a0*16*y + a1*4)*y + a2)*4*y + a3)*4*y + a4
888 : * = (((a0*32*y + a1*8)*y + a2*2)*2*y + a3)*4*y + a4
889 : * = ((( c0*y + c1)*y + c2)*2*y + c3)*4*y + c4,
890 : * where y = x/4
891 : * and a0 = -0.00262378195660943, a1 = 0.04089687039888652,
892 : * a2 = -0.25631148958325911, a3 = 0.81685854627399479,
893 : * a4 = 0.21358070563097167
894 : * */
895 40 : w = extract_l( L_shr( arg, 2 ) ); /* Q15 y = x/4 */
896 40 : acc = 702602883L; /* Q31 c1 = a1*8 */
897 40 : move32();
898 40 : z = mac_r( acc, w, -2751 ); /* Q15 c0 = a0*32 */
899 40 : acc = -1100849465L; /* Q31 c2 = a2*2 */
900 40 : move32();
901 40 : z = mac_r( acc, w, z ); /* z (in:Q15, out:Q14) */
902 40 : acc = 877095185L; /* Q30 c3 = a3 */
903 40 : move32();
904 40 : z = mac_r( acc, w, z ); /* z (in:Q14, out:Q12) */
905 40 : acc = 57332634L; /* Q28 c4 = a4 */
906 40 : move32();
907 40 : acc = L_mac( acc, w, z ); /* Q28 */
908 40 : angle = extract_l( L_shr( acc, ( 28 - 14 ) ) ); /* Q14 result of atan(x) where 2 <= x < 4 */
909 : }
910 126430 : ELSE IF( L_shr( arg, 15 ) != 0 )
911 : /*==========================*
912 : * 1.0 <= x < 2.0 *
913 : *==========================*/
914 : {
915 : /* interval: [0.999, 2.001]
916 : * atan(x) ~ (((a0*x + 1)*x + a2)*x + a3)*x + a4
917 : * = (((a0*2*y + a1)*2*y + a2)*2*y + a3)*2*y + a4 Substitute 2*y -> x
918 : * = (((a0*4*y + a1*2)*y + a2)*2*y + a3)*2*y + a4
919 : * = (((a0*4*y + a1*2)*y + a2)*y + a3/2)*4*y + a4
920 : * = ((( c0*y + c1)*y + c2)*y + c3)*4*y + c4,
921 : * where y = x/2
922 : * and a0 = -0.0160706457245251, a1 = 0.1527106504065224,
923 : * a2 = -0.6123208404800871, a3 = 1.3307896976322915,
924 : * a4 = -0.0697089375247448
925 : */
926 124717 : w = extract_l( L_shr( arg, 1 ) ); /* Q15 y= x/2 */
927 124717 : acc = 655887249L; /* Q31 c1 = a1*2 */
928 124717 : move32();
929 124717 : z = mac_r( acc, w, -2106 ); /* Q15 c0 = a0*4 */
930 124717 : acc = -1314948992L; /* Q31 c2 = a2 */
931 124717 : move32();
932 124717 : z = mac_r( acc, w, z );
933 124717 : acc = 1428924557L; /* Q31 c3 = a3/2 */
934 124717 : move32();
935 124717 : z = mac_r( acc, w, z ); /* z (in:Q15, out:Q13) */
936 124717 : acc = -37424701L; /* Q29 c4 = a4 */
937 124717 : move32();
938 124717 : acc = L_mac( acc, w, z ); /* Q29 */
939 124717 : angle = extract_l( L_shr( acc, ( 29 - 14 ) ) ); /* Q14 result of atan(x) where 1 <= x < 2 */
940 : }
941 : ELSE
942 : /*==========================*
943 : * 0.0 <= x < 1.0 *
944 : *==========================*/
945 : {
946 : /* interval: [-0.001, 1.001]
947 : * atan(x) ~ ((((a0*x + a1)*x + a2)*x + a3)*x + a4)*x + a5
948 : * = ((((a0*2*x + a1*2)*x/2 + a2)*x + a3)*x + a4)*x + a5
949 : * = (((( c0*x + c1)*x/2 + c2)*x + c3)*x + c4)*x + c5
950 : * where
951 : * a0 = -5.41182677118661e-02, a1 = 2.76690449232515e-01,
952 : * a2 = -4.63358392562492e-01, a3 = 2.87188466598566e-02,
953 : * a4 = 9.97438122814383e-01, a5 = 5.36158556179092e-05.
954 : */
955 1713 : w = extract_l( arg ); /* Q15 */
956 1713 : acc = 1188376431L; /* Q31 c1 = a1*2 */
957 1713 : move32();
958 1713 : z = mac_r( acc, w, -3547 ); /* Q15 c0 = a0*2 */
959 1713 : acc = -995054571L; /* Q31 c2 = a2 */
960 1713 : move32();
961 1713 : z = extract_h( L_mac0( acc, w, z ) ); /* Q15 non-fractional mode multiply */
962 1713 : acc = 61673254L; /* Q31 c3 = a3 */
963 1713 : move32();
964 1713 : z = mac_r( acc, w, z );
965 1713 : acc = 2141982059L; /* Q31 c4 = a4 */
966 1713 : move32();
967 1713 : z = mac_r( acc, w, z );
968 1713 : acc = 115139L; /* Q31 c5 = a5 */
969 1713 : move32();
970 1713 : acc = L_mac( acc, w, z ); /* Q31 */
971 1713 : angle = extract_l( L_shr( acc, 31 - 14 ) ); /* Q14 result of atan(x), where 0 <= x < 1 */
972 : }
973 :
974 126470 : return angle; /* Q14 between 0 and EVS_PI/2 radian. */
975 : }
976 :
977 : #undef WMC_TOOL_SKIP
|
