Line data Source code
1 : /*====================================================================================
2 : EVS Codec 3GPP TS26.452 Aug 12, 2021. Version 16.3.0
3 : ====================================================================================*/
4 :
5 : #include <assert.h>
6 : #include <stdint.h>
7 : #include "options.h"
8 : #include "basop_util.h"
9 : #include "options.h"
10 : #include "rom_basop_util.h"
11 : #include "rom_com.h"
12 : #include "prot_fx.h"
13 : #include "prot_fx_enc.h"
14 : #include "ivas_prot_fx.h"
15 :
16 : #define FFT_SCALING_512 1073741824 // Q22
17 : #define FFT_SCALING_640 1342177280 // Q22
18 :
19 : #define DELTA_SHIFT 2
20 : #define DELTA_SHIFT_LD64 67108864l /*DELTA_SHIFT/64.0 Q31*/
21 : #define CNG_HS 4 /* 4 bit headroom for dot product */
22 : #define CNG_S 6 /* 1 sign bit, 6 bit integer part, 9 bit frational part for input and output data */
23 :
24 :
25 : /*-------------------------------------------------------------------
26 : * Local function prototypes
27 : *-------------------------------------------------------------------*/
28 : static void getmidbands( const Word16 *part, const Word16 npart, Word16 *midband, Word16 *psize, Word16 *psize_norm, Word16 *psize_norm_exp, Word16 *psize_inv );
29 :
30 :
31 : /*-------------------------------------------------------------------
32 : * createFdCngCom()
33 : *
34 : * Create an instance of type FD_CNG_COM
35 : *-------------------------------------------------------------------*/
36 7135 : ivas_error createFdCngCom_fx(
37 : HANDLE_FD_CNG_COM *hFdCngCom /* i/o: FD_CNG structure containing all buffers and variables */
38 : )
39 : {
40 : HANDLE_FD_CNG_COM hs;
41 :
42 : /* Allocate memory */
43 7135 : hs = (HANDLE_FD_CNG_COM) malloc( sizeof( FD_CNG_COM ) );
44 :
45 7135 : if ( hs == NULL )
46 : {
47 0 : return IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for FD CNG COM" );
48 : }
49 :
50 7135 : *hFdCngCom = hs;
51 :
52 7135 : return IVAS_ERR_OK;
53 : }
54 : /*-------------------------------------------------------------------
55 : * initFdCngCom()
56 : *
57 : *
58 : *-------------------------------------------------------------------*/
59 :
60 7132 : void ivas_initFdCngCom_fx( HANDLE_FD_CNG_COM hFdCngCom, Word16 scale )
61 : {
62 : /* Calculate CLDFB scaling factor */
63 : /* shl(i_mult2(scale, scale), 3) does not fit in 16 bit */
64 : /*hFdCngCom->scalingFactor = div_s(1, shl(i_mult2(scale, scale), 3));*/
65 7132 : assert( 2048 /*1.0/(1<<4) Q15*/ < mult( scale, scale ) );
66 : /* Exponent invScalingFactor: -16 = -(2*7 (scale) + 2 (8.0) */
67 7132 : hFdCngCom->invScalingFactor = shl( mult( scale, scale ), 1 );
68 : /* Exponent scalingFactor: -15 = -(2*7 (scale) + 2 (8.0) - 1 (1.0)) */
69 7132 : hFdCngCom->scalingFactor = div_s( 0x4000, hFdCngCom->invScalingFactor );
70 :
71 : /* Initialize the overlap-add */
72 7132 : set16_fx( hFdCngCom->timeDomainBuffer, 0, L_FRAME16k );
73 7132 : hFdCngCom->olapBufferAna = NULL;
74 7132 : set16_fx( hFdCngCom->olapBufferAna_fx, 0, FFTLEN );
75 :
76 7132 : move16();
77 7132 : set16_fx( hFdCngCom->olapBufferSynth, 0, FFTLEN );
78 7132 : hFdCngCom->olapBufferSynth2 = NULL;
79 7132 : move16();
80 :
81 : /* Initialize the comfort noise generation */
82 7132 : set32_fx( hFdCngCom->fftBuffer, 0, FFTLEN );
83 7132 : set32_fx( hFdCngCom->cngNoiseLevel, 0, FFTCLDFBLEN );
84 7132 : set16_fx( &hFdCngCom->cngNoiseLevelExp, 0, 1 );
85 :
86 : /* Initialize quantizer */
87 7132 : set32_fx( hFdCngCom->sidNoiseEst, 0, NPART );
88 7132 : set16_fx( &hFdCngCom->sidNoiseEstExp, 0, 1 );
89 7132 : set16_fx( hFdCngCom->A_cng, 0, M + 1 );
90 7132 : hFdCngCom->A_cng[0] = 4096 /*1.f Q12*/; /* 3Q12 */
91 7132 : move16();
92 :
93 : /* Set some counters and flags */
94 7132 : hFdCngCom->inactive_frame_counter = 0; /* Either SID or zero frames */
95 7132 : move16();
96 7132 : hFdCngCom->active_frame_counter = 0;
97 7132 : move16();
98 7132 : hFdCngCom->frame_type_previous = ACTIVE_FRAME;
99 7132 : move16();
100 7132 : hFdCngCom->flag_noisy_speech = 0;
101 7132 : move16();
102 7132 : hFdCngCom->likelihood_noisy_speech = 0;
103 7132 : move16();
104 7132 : hFdCngCom->numCoreBands = 0;
105 7132 : move16();
106 7132 : hFdCngCom->stopBand = 0;
107 7132 : move16();
108 7132 : hFdCngCom->startBand = 0;
109 7132 : move16();
110 7132 : hFdCngCom->stopFFTbin = 0;
111 7132 : move16();
112 7132 : hFdCngCom->frameSize = 0;
113 7132 : move16();
114 7132 : hFdCngCom->fftlen = 0;
115 7132 : move16();
116 7132 : hFdCngCom->seed = 0;
117 7132 : move16();
118 7132 : hFdCngCom->seed2 = 1;
119 7132 : move16();
120 7132 : hFdCngCom->seed3 = 2;
121 7132 : move16();
122 7132 : hFdCngCom->CngBitrate = -1;
123 7132 : move16();
124 7132 : hFdCngCom->olapBufferSynth_exp = 0;
125 7132 : move16();
126 :
127 : /* Initialize noise estimation algorithm */
128 7132 : set32_fx( hFdCngCom->periodog, 0, PERIODOGLEN );
129 7132 : mhvals( MSNUMSUBFR * MSSUBFRLEN, &( hFdCngCom->msM_win ) );
130 7132 : mhvals( MSSUBFRLEN, &( hFdCngCom->msM_subwin ) );
131 7132 : set32_fx( hFdCngCom->msPeriodogSum, 0, 2 );
132 7132 : hFdCngCom->msPeriodogSum_exp[0] = 0;
133 7132 : move16();
134 7132 : hFdCngCom->msPeriodogSum_exp[1] = 0;
135 7132 : move16();
136 7132 : set32_fx( hFdCngCom->msPsdSum, 0, 2 );
137 7132 : set16_fx( hFdCngCom->msSlope, 0, 2 );
138 7132 : set32_fx( hFdCngCom->msQeqInvAv, 0, 2 );
139 7132 : hFdCngCom->msQeqInvAv_exp[0] = 0;
140 7132 : move16();
141 7132 : hFdCngCom->msQeqInvAv_exp[1] = 0;
142 7132 : move16();
143 7132 : hFdCngCom->msFrCnt_init_counter = 0;
144 7132 : move16();
145 7132 : hFdCngCom->msFrCnt_init_thresh = 1;
146 7132 : move16();
147 7132 : hFdCngCom->init_old = 0;
148 7132 : move16();
149 7132 : hFdCngCom->offsetflag = 0;
150 7132 : move16();
151 7132 : hFdCngCom->msFrCnt = MSSUBFRLEN;
152 7132 : move16();
153 7132 : hFdCngCom->msMinBufferPtr = 0;
154 7132 : move16();
155 7132 : hFdCngCom->msAlphaCor[0] = 644245120l /*0.3f Q31*/;
156 7132 : move16();
157 7132 : hFdCngCom->msAlphaCor[1] = 644245120l /*0.3f Q31*/;
158 7132 : move16();
159 7132 : set16_fx( hFdCngCom->psize, 0, NPART );
160 : /* Initialize exponents */
161 7132 : hFdCngCom->exp_cldfb_periodog = 0;
162 7132 : move16();
163 :
164 7132 : hFdCngCom->coherence_fx = 16384; /* 0.5 in Q15 */
165 7132 : move16();
166 :
167 7132 : set32_fx( hFdCngCom->olapBufferSynth_fx, 0, FFTLEN );
168 7132 : set32_fx( hFdCngCom->olapBufferSynth2_fx, 0, FFTLEN );
169 7132 : set32_fx( hFdCngCom->exc_cng_32fx, 0, L_FRAME16k );
170 7132 : set16_fx( hFdCngCom->exc_cng, 0, L_FRAME16k );
171 :
172 7132 : return;
173 : }
174 3 : void initFdCngCom( HANDLE_FD_CNG_COM hFdCngCom, Word16 scale )
175 : {
176 : /* Calculate CLDFB scaling factor */
177 : /* shl(i_mult2(scale, scale), 3) does not fit in 16 bit */
178 : /*hFdCngCom->scalingFactor = div_s(1, shl(i_mult2(scale, scale), 3));*/
179 3 : assert( 2048 /*1.0/(1<<4) Q15*/ < mult( scale, scale ) );
180 : /* Exponent invScalingFactor: -16 = -(2*7 (scale) + 2 (8.0) */
181 3 : hFdCngCom->invScalingFactor = shl( mult( scale, scale ), 1 );
182 3 : move16();
183 : /* Exponent scalingFactor: -15 = -(2*7 (scale) + 2 (8.0) - 1 (1.0)) */
184 3 : hFdCngCom->scalingFactor = div_s( 0x4000, hFdCngCom->invScalingFactor );
185 3 : move16();
186 :
187 : /* Initialize the overlap-add */
188 3 : set16_fx( hFdCngCom->timeDomainBuffer, 0, L_FRAME16k );
189 3 : hFdCngCom->olapBufferAna = NULL;
190 3 : set16_fx( hFdCngCom->olapBufferAna_fx, 0, FFTLEN );
191 :
192 3 : set16_fx( hFdCngCom->olapBufferSynth, 0, FFTLEN );
193 3 : hFdCngCom->olapBufferSynth2 = NULL;
194 :
195 : /* Initialize the comfort noise generation */
196 3 : set32_fx( hFdCngCom->fftBuffer, 0, FFTLEN );
197 3 : set32_fx( hFdCngCom->cngNoiseLevel, 0, FFTCLDFBLEN );
198 3 : set16_fx( &hFdCngCom->cngNoiseLevelExp, 0, 1 );
199 :
200 : /* Initialize quantizer */
201 3 : set32_fx( hFdCngCom->sidNoiseEst, 0, NPART );
202 3 : set16_fx( &hFdCngCom->sidNoiseEstExp, 0, 1 );
203 3 : set16_fx( hFdCngCom->A_cng, 0, M + 1 );
204 3 : hFdCngCom->A_cng[0] = 4096 /*1.f Q12*/; /* 3Q12 */
205 3 : move16();
206 :
207 : /* Set some counters and flags */
208 3 : hFdCngCom->inactive_frame_counter = 0; /* Either SID or zero frames */
209 3 : move16();
210 3 : hFdCngCom->active_frame_counter = 0;
211 3 : move16();
212 3 : hFdCngCom->frame_type_previous = ACTIVE_FRAME;
213 3 : move16();
214 3 : hFdCngCom->flag_noisy_speech = 0;
215 3 : move16();
216 3 : hFdCngCom->likelihood_noisy_speech = 0;
217 3 : move16();
218 3 : hFdCngCom->numCoreBands = 0;
219 3 : move16();
220 3 : hFdCngCom->stopBand = 0;
221 3 : move16();
222 3 : hFdCngCom->startBand = 0;
223 3 : move16();
224 3 : hFdCngCom->stopFFTbin = 0;
225 3 : move16();
226 3 : hFdCngCom->frameSize = 0;
227 3 : move16();
228 3 : hFdCngCom->fftlen = 0;
229 3 : move16();
230 3 : hFdCngCom->seed = 0;
231 3 : move16();
232 : // hFdCngCom->seed2 = 1;
233 : // move16();
234 : // hFdCngCom->seed3 = 2;
235 : // move16();
236 3 : hFdCngCom->CngBitrate = -1;
237 3 : move16();
238 :
239 : /* Initialize noise estimation algorithm */
240 3 : set32_fx( hFdCngCom->periodog, 0, PERIODOGLEN );
241 3 : mhvals( MSNUMSUBFR * MSSUBFRLEN, &( hFdCngCom->msM_win ) );
242 3 : mhvals( MSSUBFRLEN, &( hFdCngCom->msM_subwin ) );
243 3 : set32_fx( hFdCngCom->msPeriodogSum, 0, 2 );
244 3 : hFdCngCom->msPeriodogSum_exp[0] = 0;
245 3 : move16();
246 3 : hFdCngCom->msPeriodogSum_exp[1] = 0;
247 3 : move16();
248 3 : set32_fx( hFdCngCom->msPsdSum, 0, 2 );
249 3 : set16_fx( hFdCngCom->msSlope, 0, 2 );
250 3 : set32_fx( hFdCngCom->msQeqInvAv, 0, 2 );
251 3 : hFdCngCom->msQeqInvAv_exp[0] = 0;
252 3 : move16();
253 3 : hFdCngCom->msQeqInvAv_exp[1] = 0;
254 3 : move16();
255 3 : hFdCngCom->msFrCnt_init_counter = 0;
256 3 : move16();
257 3 : hFdCngCom->msFrCnt_init_thresh = 1;
258 3 : move16();
259 3 : hFdCngCom->init_old = 0;
260 3 : move16();
261 3 : hFdCngCom->offsetflag = 0;
262 3 : move16();
263 3 : hFdCngCom->msFrCnt = MSSUBFRLEN;
264 3 : move16();
265 3 : hFdCngCom->msMinBufferPtr = 0;
266 3 : move16();
267 3 : hFdCngCom->msAlphaCor[0] = 644245120l /*0.3f Q31*/;
268 3 : move16();
269 3 : hFdCngCom->msAlphaCor[1] = 644245120l /*0.3f Q31*/;
270 3 : move16();
271 3 : set16_fx( hFdCngCom->psize, 0, NPART );
272 : /* Initialize exponents */
273 3 : hFdCngCom->exp_cldfb_periodog = 0;
274 3 : move16();
275 :
276 3 : return;
277 : }
278 :
279 : /*-------------------------------------------------------------------
280 : * deleteFdCngCom()
281 : *
282 : * Delete an instance of type FD_CNG_COM
283 : *-------------------------------------------------------------------*/
284 7135 : void deleteFdCngCom_fx( HANDLE_FD_CNG_COM *hFdCngCom ) /* i/o: Contains the variables related to the CLDFB-based CNG process */
285 : {
286 : HANDLE_FD_CNG_COM hsCom;
287 7135 : hsCom = *hFdCngCom;
288 7135 : move16();
289 7135 : IF( hsCom != NULL )
290 : {
291 7135 : free( hsCom );
292 7135 : *hFdCngCom = NULL;
293 7135 : move16();
294 : }
295 7135 : }
296 :
297 : /*-------------------------------------------------------------------
298 : * initPartitions()
299 : *
300 : * Initialize the spectral partitioning
301 : *-------------------------------------------------------------------*/
302 1125334 : void initPartitions(
303 : const Word16 *part_in,
304 : const Word16 npart_in,
305 : const Word16 startBand,
306 : const Word16 stopBand,
307 : Word16 *part_out,
308 : Word16 *npart_out,
309 : Word16 *midband,
310 : Word16 *psize,
311 : Word16 *psize_norm,
312 : Word16 *psize_norm_exp,
313 : Word16 *psize_inv,
314 : const Word16 stopBandFR )
315 : {
316 : Word16 i, j, len_out, tmp16;
317 :
318 :
319 1125334 : IF( part_in != NULL )
320 : {
321 1125334 : len_out = 0;
322 1125334 : move16();
323 1125334 : IF( GT_16( stopBandFR, startBand ) )
324 : {
325 562667 : len_out = sub( stopBandFR, startBand );
326 21944013 : FOR( i = 0; i < len_out; i++ )
327 : {
328 21381346 : part_out[i] = i;
329 21381346 : move16();
330 : }
331 : }
332 36570080 : FOR( j = 0; j < npart_in; j++ )
333 : {
334 35444746 : IF( GE_16( part_in[j], stopBand ) )
335 : {
336 0 : BREAK;
337 : }
338 35444746 : tmp16 = sub( part_in[j], startBand );
339 35444746 : test();
340 35444746 : if ( GE_16( part_in[j], stopBandFR ) && tmp16 >= 0 )
341 : {
342 26442074 : part_out[len_out++] = tmp16;
343 26442074 : move16();
344 : }
345 : }
346 : }
347 : ELSE
348 : {
349 0 : len_out = sub( stopBand, startBand );
350 0 : FOR( i = 0; i < len_out; i++ )
351 : {
352 0 : part_out[i] = i;
353 0 : move16();
354 : }
355 : }
356 :
357 1125334 : *npart_out = len_out;
358 1125334 : move16();
359 1125334 : getmidbands( part_out, len_out, midband, psize, psize_norm, psize_norm_exp, psize_inv );
360 1125334 : }
361 :
362 :
363 : /*-------------------------------------------------------------------
364 : * compress_range()
365 : *
366 : * Apply some dynamic range compression based on the log
367 : *-------------------------------------------------------------------*/
368 44318 : void compress_range(
369 : Word32 *in, // Q(31 - in_exp)
370 : Word16 in_exp,
371 : Word16 *out,
372 : const Word16 len )
373 : {
374 : Word16 i;
375 : Word32 in_s;
376 : Word32 one_s;
377 : Word32 in_exp32;
378 : Word32 L_tmp;
379 :
380 :
381 : /* out = log2( 1 + in ) */
382 44318 : IF( in_exp >= 0 )
383 : {
384 34685 : one_s = L_shr( 1073741824l /*0.5 Q31*/, in_exp ); // Q(31 - in_exp)
385 34685 : in_exp32 = L_shl( L_deposit_l( add( in_exp, 1 ) ), WORD32_BITS - 1 - LD_DATA_SCALE );
386 2163922 : FOR( i = 0; i < len; i++ )
387 : {
388 2129237 : in_s = L_add( L_shr( in[i], 1 ), one_s ); // Q(31 - in_exp)
389 2129237 : L_tmp = L_add( BASOP_Util_Log2( in_s ), in_exp32 );
390 2129237 : if ( in_s == 0 )
391 : {
392 56965 : out[i] = 0;
393 56965 : move16();
394 : }
395 2129237 : if ( in_s != 0 )
396 : {
397 2072272 : out[i] = extract_h( L_tmp );
398 2072272 : move16();
399 : }
400 2129237 : if ( out[i] == 0 )
401 : {
402 104827 : out[i] = 1;
403 104827 : move16();
404 : }
405 : }
406 : }
407 : ELSE
408 : {
409 9633 : in_exp = sub( in_exp, 1 );
410 9633 : in_exp32 = L_shl( L_deposit_l( 1 ), WORD32_BITS - 1 - LD_DATA_SCALE );
411 603086 : FOR( i = 0; i < len; i++ )
412 : {
413 593453 : L_tmp = L_add( BASOP_Util_Log2( L_add( L_shl( in[i], in_exp ), 1073741824l /*0.5 Q31*/ ) ), in_exp32 );
414 593453 : if ( in[i] == 0 )
415 : {
416 6960 : out[i] = 0;
417 6960 : move16();
418 : }
419 593453 : if ( in[i] != 0 )
420 : {
421 586493 : out[i] = extract_h( L_tmp );
422 586493 : move16();
423 : }
424 593453 : if ( out[i] == 0 )
425 : {
426 167614 : out[i] = 1;
427 167614 : move16();
428 : }
429 : }
430 : }
431 44318 : }
432 :
433 : /*-------------------------------------------------------------------
434 : * expand_range()
435 : *
436 : * Apply some dynamic range expansion to undo the compression
437 : *-------------------------------------------------------------------*/
438 15 : void expand_range(
439 : Word16 *in, // Q15
440 : Word32 *out, // Q(31 - out_exp)
441 : Word16 *out_exp,
442 : const Word16 len )
443 : {
444 : Word16 i;
445 : Word16 s;
446 : Word16 tmp;
447 : Word32 one_s, tmp32;
448 : Word16 maxVal;
449 : Word16 maxVal2;
450 :
451 :
452 15 : maxVal = 0;
453 15 : move16();
454 937 : FOR( i = 0; i < len; i++ )
455 : {
456 922 : maxVal = s_max( maxVal, in[i] ); // Q15
457 : }
458 :
459 15 : maxVal2 = maxVal;
460 15 : move16();
461 15 : s = 0;
462 15 : move16();
463 61 : WHILE( maxVal >= 0 )
464 : {
465 46 : maxVal = sub( maxVal, 512 /*0.015625 Q15*/ ); // Q15
466 46 : s = add( s, 1 );
467 : }
468 15 : tmp = sub( maxVal2, maxVal );
469 :
470 15 : one_s = L_shr( 1073741824l /*0.5 Q31*/, sub( s, 1 ) ); // Q(31 - sub( s, 1 ))
471 15 : tmp32 = L_shr( 726941l /*0.0003385080526823181 Q31*/, s );
472 : /* out = (2^(in) - 1) */
473 937 : FOR( i = 0; i < len; i++ )
474 : {
475 922 : out[i] = L_sub( BASOP_Util_InvLog2( L_deposit_h( sub( in[i], tmp ) ) ), one_s );
476 922 : move32();
477 922 : if ( out[i] == 0 )
478 : {
479 50 : out[i] = tmp32;
480 50 : move32();
481 : }
482 : }
483 15 : *out_exp = s;
484 15 : move16();
485 15 : }
486 :
487 44303 : void expand_range_fx(
488 : Word32 *in, // Q25
489 : Word32 *out, // exp:out_exp
490 : Word16 *out_exp,
491 : const Word16 len )
492 : {
493 : Word16 i, tmp_e;
494 : Word32 maxVal, tmp32;
495 : Word16 maxOutExp;
496 :
497 44303 : const Word32 low_lim = 726940; /* 0.0003385080526823181f in Q31 */
498 44303 : move32();
499 :
500 : /* Find max possible output exponent. */
501 44303 : maxVal = 0;
502 44303 : move32();
503 2766071 : FOR( i = 0; i < len; i++ )
504 : {
505 2721768 : maxVal = L_max( maxVal, in[i] );
506 : }
507 44303 : tmp32 = BASOP_util_Pow2( maxVal, Q31 - Q25, &tmp_e );
508 44303 : maxOutExp = tmp_e;
509 44303 : move16();
510 :
511 : /* out = (2^(in) - 1) */
512 2766071 : FOR( i = 0; i < len; i++ )
513 : {
514 2721768 : tmp32 = BASOP_util_Pow2( in[i], Q31 - Q25, &tmp_e ); // 2^x
515 2721768 : tmp32 = L_sub( tmp32, L_shl( 1, sub( 31, tmp_e ) ) ); // 2^x - 1
516 2721768 : tmp32 = L_shr( tmp32, sub( maxOutExp, tmp_e ) ); // make exp same as maxExpOut
517 :
518 2721768 : out[i] = tmp32;
519 2721768 : move32();
520 :
521 2721768 : Word32 tmp_low_lim = L_shr( low_lim, maxOutExp );
522 2721768 : if ( LT_32( out[i], tmp_low_lim ) )
523 : {
524 14424 : out[i] = tmp_low_lim;
525 14424 : move32();
526 : }
527 : }
528 :
529 44303 : *out_exp = maxOutExp;
530 44303 : move16();
531 :
532 44303 : return;
533 : }
534 :
535 : /*-------------------------------------------------------------------
536 : * expand_range_var_e()
537 : *
538 : * Apply some dynamic range expansion to undo the compression. Similar
539 : * to expand_range but has variable input exponent.
540 : *
541 : *-------------------------------------------------------------------*/
542 0 : void expand_range_var_exp(
543 : Word16 *in,
544 : Word16 in_exp,
545 : Word32 *out,
546 : Word16 *out_exp,
547 : const Word16 len )
548 : {
549 : Word16 i, tmp_e;
550 : Word32 tmp32;
551 : Word16 maxVal, maxOutExp;
552 :
553 0 : const Word32 low_lim = 726940; /* 0.0003385080526823181f in Q31 */
554 0 : move32();
555 :
556 : /* Find max possible output exponent. */
557 0 : maxVal = 0;
558 0 : move16();
559 0 : FOR( i = 0; i < len; i++ )
560 : {
561 0 : maxVal = s_max( maxVal, in[i] );
562 : }
563 0 : tmp32 = BASOP_util_Pow2( L_deposit_h( maxVal ), in_exp, &tmp_e );
564 0 : maxOutExp = tmp_e;
565 0 : move16();
566 :
567 : /* out = (2^(in) - 1) */
568 0 : FOR( i = 0; i < len; i++ )
569 : {
570 0 : tmp32 = BASOP_util_Pow2( L_deposit_h( in[i] ), in_exp, &tmp_e ); // 2^x
571 0 : tmp32 = L_sub( tmp32, L_shl( 1, sub( 31, tmp_e ) ) ); // 2^x - 1
572 0 : tmp32 = L_shr( tmp32, sub( maxOutExp, tmp_e ) ); // make exp same as maxExpOut
573 :
574 0 : out[i] = tmp32;
575 0 : move32();
576 :
577 0 : Word32 tmp_low_lim = L_shr( low_lim, maxOutExp );
578 0 : if ( LT_32( out[i], tmp_low_lim ) )
579 : {
580 0 : out[i] = tmp_low_lim;
581 0 : move32();
582 : }
583 : }
584 :
585 0 : *out_exp = maxOutExp;
586 0 : move16();
587 :
588 0 : return;
589 : }
590 :
591 : /*-------------------------------------------------------------------
592 : * minimum_statistics()
593 : *
594 : * Noise estimation using Minimum Statistics (MS)
595 : *-------------------------------------------------------------------*/
596 :
597 15 : void minimum_statistics(
598 : Word16 len, /* i : Total number of partitions (CLDFB or FFT) */
599 : Word16 lenFFT, /* i : Number of FFT partitions */
600 : Word16 *psize, /* i : Partition sizes, fractional */
601 : Word16 *msPeriodog, /* i : Periodogram (energies) */
602 : Word16 *msNoiseFloor, /* i/o: Noise floors (energies) */
603 : Word16 *msNoiseEst, /* i/o: Noise estimates (energies) */
604 : Word32 *msAlpha, /* i/o: Forgetting factors */
605 : Word16 *msPsd, /* i/o: Power Spectral Density (smoothed periodogram => energies) */
606 : Word16 *msPsdFirstMoment, /* i/o: PSD statistics of 1st order (energy means) */
607 : Word32 *msPsdSecondMoment, /* i/o: PSD statistics of 2nd order (energy variances) */
608 : Word32 *msMinBuf, /* i/o: Buffer of minima (energies) */
609 : Word32 *msBminWin, /* o : Bias correction factors */
610 : Word32 *msBminSubWin, /* o : Bias correction factors */
611 : Word32 *msCurrentMin, /* i/o: Local minima (energies) */
612 : Word32 *msCurrentMinOut, /* i/o: Local minima (energies) */
613 : Word32 *msCurrentMinSubWindow, /* i/o: Local minima (energies) */
614 : Word16 *msLocalMinFlag, /* i : Binary flag */
615 : Word16 *msNewMinFlag, /* i : Binary flag */
616 : Word16 *msPeriodogBuf, /* i/o: Buffer of periodograms (energies) */
617 : Word16 *msPeriodogBufPtr, /* i/o: Counter */
618 : HANDLE_FD_CNG_COM hFdCngCom /* i/o: FD_CNG structure containing all buffers and variables */
619 : )
620 : {
621 : Word16 i, j, k, s, s1, s2, s3;
622 : Word16 len2;
623 : Word16 current_len;
624 : Word16 start, stop, cnt;
625 : Word16 totsize;
626 : Word16 inv_totsize;
627 :
628 : Word32 tmp, tmp0, tmp1;
629 : Word32 scalar, scalar2, scalar3;
630 : Word32 snr;
631 : Word32 msAlphaHatMin2;
632 : Word32 BminCorr;
633 : Word32 QeqInvAv;
634 : Word32 *ptr;
635 : Word32 *msPsdSum;
636 : Word32 *msPeriodogSum;
637 :
638 : Word16 tmp16, tmp16_1;
639 : Word16 beta;
640 : Word16 slope;
641 : Word16 QeqInv;
642 : Word16 scalar16;
643 : Word16 scalar216;
644 : Word16 scalar316;
645 : Word16 msM_win;
646 : Word16 msM_subwin;
647 : Word16 msAlphaCorAlpha;
648 : Word16 msAlphaCorAlpha2;
649 : Word16 msPeriodogSum16;
650 : Word16 msNoiseFloor16;
651 : #ifndef ISSUE_1836_replace_overflow_libcom
652 : #ifdef BASOP_NOGLOB_DECLARE_LOCAL
653 : Flag Overflow = 0;
654 : move32();
655 : #endif
656 : #endif
657 :
658 15 : len2 = i_mult( MSNUMSUBFR, len );
659 :
660 15 : msM_win = hFdCngCom->msM_win;
661 15 : move16();
662 15 : msM_subwin = hFdCngCom->msM_subwin;
663 15 : move16();
664 :
665 15 : msAlphaCorAlpha = MSALPHACORALPHA;
666 15 : move16();
667 15 : msAlphaCorAlpha2 = MSALPHACORALPHA2;
668 15 : move16();
669 :
670 15 : msPsdSum = hFdCngCom->msPsdSum;
671 15 : msPeriodogSum = hFdCngCom->msPeriodogSum;
672 :
673 : /* No minimum statistics at initialization */
674 15 : IF( LT_16( hFdCngCom->msFrCnt_init_counter, hFdCngCom->msFrCnt_init_thresh ) )
675 : {
676 6 : Copy( msPeriodog, msPsd, len ); /* 6Q9 */
677 6 : Copy( msPeriodog, msNoiseFloor, len ); /* 6Q9 */
678 6 : Copy( msPeriodog, msNoiseEst, len ); /* 6Q9 */
679 6 : Copy( msPeriodog, msPsdFirstMoment, len ); /* 6Q9 */
680 :
681 6 : set32_fx( msPsdSecondMoment, 0l /*0.0 Q31*/, len );
682 6 : msPeriodogSum[0] = dotp_s_fx( msPeriodog, psize, lenFFT, CNG_HS );
683 6 : move32();
684 6 : msPsdSum[0] = msPeriodogSum[0]; /* 16Q15 */
685 6 : move32();
686 :
687 6 : IF( LT_16( lenFFT, len ) )
688 : {
689 0 : msPeriodogSum[1] = dotp_s_fx( msPeriodog + lenFFT, psize + lenFFT, sub( len, lenFFT ), CNG_HS );
690 0 : move32();
691 0 : msPsdSum[1] = msPeriodogSum[1]; /* 16Q15 */
692 0 : move32();
693 : }
694 :
695 : /* Increment frame counter at initialization */
696 : /* Some frames are sometimes zero at initialization => ignore them */
697 6 : IF( LT_16( msPeriodog[0], hFdCngCom->init_old ) )
698 : {
699 2 : set32_fx( msCurrentMinOut, 2147483647l /*1.0 Q31*/, len ); /* 16Q15 */
700 2 : set32_fx( msCurrentMin, 2147483647l /*1.0 Q31*/, len ); /* 16Q15 */
701 2 : set32_fx( msMinBuf, 2147483647l /*1.0 Q31*/, len2 ); /* 16Q15 */
702 2 : set32_fx( msCurrentMinSubWindow, 2147483647l /*1.0 Q31*/, len ); /* 16Q15 */
703 :
704 2 : hFdCngCom->msFrCnt_init_counter = add( hFdCngCom->msFrCnt_init_counter, 1 );
705 2 : move16();
706 : }
707 6 : hFdCngCom->init_old = msPeriodog[0]; /* 6Q9 */
708 6 : move16();
709 : }
710 : ELSE
711 : {
712 : /* Consider the FFT and CLDFB bands separately
713 : - first iteration for FFT bins,
714 : - second one for CLDFB bands in SWB mode */
715 9 : cnt = 0;
716 9 : move16();
717 9 : start = 0;
718 9 : move16();
719 9 : stop = lenFFT;
720 9 : move16();
721 9 : totsize = sub( hFdCngCom->stopFFTbin, hFdCngCom->startBand );
722 18 : WHILE( GT_16( stop, start ) )
723 : {
724 9 : current_len = sub( stop, start );
725 :
726 : /* Compute smoothed correction factor for PSD smoothing */
727 :
728 : /* msPeriodogSum[cnt] with format 16Q15 */
729 9 : msPeriodogSum[cnt] = dotp_s_fx( msPeriodog + start, psize + start, current_len, CNG_HS );
730 9 : move32();
731 :
732 9 : IF( msPeriodogSum[cnt] == 0 )
733 : {
734 0 : hFdCngCom->msAlphaCor[cnt] = Mpy_32_16_1( hFdCngCom->msAlphaCor[cnt], msAlphaCorAlpha );
735 0 : move32();
736 : }
737 : ELSE
738 : {
739 : /* calculate scalar with normalized msPeriodogSum[cnt], exponent -2*s1 */
740 9 : s1 = norm_l( msPeriodogSum[cnt] );
741 9 : msPeriodogSum16 = round_fx_sat( L_shl_sat( msPeriodogSum[cnt], s1 ) );
742 9 : scalar = L_mult( msPeriodogSum16, msPeriodogSum16 );
743 :
744 : /* calculate difference, both elements in 16Q15 format, use absolute value
745 : to avoid -1.0 x -1.0 multiplications later */
746 9 : scalar2 = L_abs( L_sub( msPsdSum[cnt], msPeriodogSum[cnt] ) );
747 :
748 9 : s2 = WORD32_BITS - 1;
749 9 : move16();
750 9 : if ( scalar2 != 0 )
751 : {
752 : /* use absolute value to avoid -1.0 x -1.0 multiplications */
753 9 : s2 = norm_l( scalar2 );
754 : }
755 : #ifdef ISSUE_1836_replace_overflow_libcom
756 9 : scalar216 = round_fx_sat( L_shl( scalar2, s2 ) );
757 : #else
758 : scalar216 = round_fx_o( L_shl_o( scalar2, s2, &Overflow ), &Overflow );
759 : #endif
760 9 : scalar2 = L_mult( scalar216, scalar216 );
761 :
762 : /* find common exponent */
763 9 : tmp16_1 = sub( s1, s2 );
764 9 : tmp16 = s_min( shl( abs_s( tmp16_1 ), 1 ), WORD32_BITS - 1 );
765 9 : if ( tmp16_1 < 0 )
766 : {
767 7 : scalar2 = L_shr( scalar2, tmp16 );
768 : }
769 9 : if ( tmp16_1 > 0 )
770 : {
771 1 : scalar = L_shr( scalar, tmp16 );
772 : }
773 :
774 :
775 : /* add scalar and scalar2, avoid overflows */
776 9 : scalar = L_shr( scalar, 1 );
777 9 : scalar2 = L_shr( scalar2, 1 );
778 9 : scalar3 = L_add( scalar, scalar2 );
779 :
780 : /* calculate division */
781 9 : scalar16 = BASOP_Util_Divide3232_uu_1616_Scale( scalar, scalar3, &s3 );
782 9 : s3 = s_max( s3, -( WORD16_BITS - 1 ) );
783 9 : scalar16 = shl( scalar16, s3 );
784 9 : scalar16 = s_max( scalar16, MSALPHACORMAX );
785 :
786 9 : hFdCngCom->msAlphaCor[cnt] = L_add( Mpy_32_16_1( hFdCngCom->msAlphaCor[cnt], msAlphaCorAlpha ),
787 : L_mult( scalar16, msAlphaCorAlpha2 ) );
788 9 : move32();
789 : }
790 :
791 : /* Compute SNR */
792 :
793 : /* msPeriodogSum[cnt] with format 16Q15 */
794 9 : snr = dotp_s_fx( msNoiseFloor + start, psize + start, current_len, CNG_HS );
795 :
796 9 : IF( GT_32( L_shr( Mpy_32_16_1( msPsdSum[cnt], 18431 /*0.56246299817 Q15*/ ), 13 ), snr ) )
797 : {
798 0 : tmp0 = BASOP_Util_Log2( msPsdSum[cnt] );
799 0 : tmp1 = BASOP_Util_Log2( snr );
800 0 : tmp1 = L_sub( tmp0, tmp1 );
801 0 : tmp1 = Mpy_32_16_1( tmp1, MSSNREXP );
802 0 : msAlphaHatMin2 = BASOP_Util_InvLog2( tmp1 );
803 : }
804 : ELSE
805 : {
806 9 : msAlphaHatMin2 = MSALPHAHATMIN;
807 9 : move32();
808 : }
809 9 : scalar = Mpy_32_16_1( hFdCngCom->msAlphaCor[cnt], MSALPHAMAX );
810 :
811 560 : FOR( j = start; j < stop; j++ )
812 : {
813 : /* Compute optimal smoothing parameter for PSD estimation */
814 551 : test();
815 551 : IF( ( scalar == 0 ) || ( msNoiseFloor[j] == 0 ) )
816 : {
817 0 : msAlpha[j] = msAlphaHatMin2;
818 0 : move32();
819 : }
820 : ELSE
821 : {
822 : /* calculate scalar2 with normalized msNoiseFloor[j], exponent -2*s1 */
823 551 : s1 = WORD16_BITS - 1;
824 551 : move16();
825 551 : if ( msNoiseFloor[j] != 0 )
826 : {
827 551 : s1 = norm_s( msNoiseFloor[j] );
828 : }
829 551 : msNoiseFloor16 = shl( msNoiseFloor[j], s1 );
830 551 : scalar2 = L_mult( msNoiseFloor16, msNoiseFloor16 );
831 :
832 : /* calculate difference, both elements in 6Q9 format, use absolute value
833 : to avoid -1.0 x -1.0 multiplications later */
834 551 : scalar316 = abs_s( sub( msPsd[j], msNoiseFloor[j] ) );
835 :
836 551 : s2 = WORD16_BITS - 1;
837 551 : move16();
838 551 : if ( scalar316 != 0 )
839 : {
840 : /* use absolute value to avoid -1.0 x -1.0 multiplications */
841 397 : s2 = norm_s( scalar316 );
842 : }
843 551 : scalar316 = shl( scalar316, s2 );
844 551 : scalar3 = L_mult( scalar316, scalar316 );
845 :
846 : /* find common exponent */
847 551 : tmp16_1 = sub( s1, s2 );
848 551 : tmp16 = s_min( shl( abs_s( tmp16_1 ), 1 ), WORD32_BITS - 1 );
849 551 : if ( tmp16_1 < 0 )
850 : {
851 460 : scalar3 = L_shr( scalar3, tmp16 );
852 : }
853 551 : if ( tmp16_1 > 0 )
854 : {
855 3 : scalar2 = L_shr( scalar2, tmp16 );
856 : }
857 :
858 : /* add scalar2 and scalar3, avoid overflows */
859 551 : scalar2 = L_shr( scalar2, 1 );
860 551 : scalar3 = L_shr( scalar3, 1 );
861 551 : scalar3 = L_add( scalar2, scalar3 );
862 :
863 : /* calculate division */
864 551 : tmp16 = BASOP_Util_Divide3232_uu_1616_Scale( scalar2, scalar3, &s3 );
865 551 : scalar2 = Mpy_32_16_1( scalar, tmp16 );
866 551 : s3 = s_max( s3, -( WORD32_BITS - 1 ) );
867 551 : scalar2 = L_shl( scalar2, s3 );
868 551 : msAlpha[j] = L_max( scalar2, msAlphaHatMin2 );
869 551 : move32();
870 : }
871 :
872 : /* Compute the PSD (smoothed periodogram) in each band */
873 551 : msPsd[j] = round_fx( L_add( Mpy_32_16_1( msAlpha[j], msPsd[j] ),
874 551 : Mpy_32_16_1( L_sub( 2147483647l /*1.0 Q31*/, msAlpha[j] ), msPeriodog[j] ) ) );
875 551 : move16();
876 : }
877 9 : msPsdSum[cnt] = dotp_s_fx( msPsd + start, psize + start, current_len, CNG_HS );
878 9 : move32();
879 :
880 9 : QeqInvAv = 0l /*0.0 Q31*/;
881 9 : move32();
882 :
883 : /* scalar: 7Q24 format */
884 9 : tmp = 1191182336l /*(float)(MSNUMSUBFR*MSSUBFRLEN)-1.0 Q24*/;
885 9 : move32();
886 9 : scalar = L_sub( tmp, L_mult( round_fx( tmp ), msM_win ) );
887 :
888 : /* scalar2: 4Q27 format */
889 9 : tmp = 1476395008l /*(float)MSSUBFRLEN-1.0 Q27*/;
890 9 : move32();
891 9 : scalar2 = L_sub( tmp, L_mult( round_fx( tmp ), msM_subwin ) );
892 :
893 560 : FOR( j = start; j < stop; j++ )
894 : {
895 : /* Compute variance of PSD */
896 551 : tmp = L_min( msAlpha[j], MSBETAMAX_SQRT );
897 :
898 551 : s1 = WORD32_BITS - 1;
899 551 : move16();
900 551 : if ( tmp != 0 )
901 : {
902 551 : s1 = norm_l( tmp );
903 : }
904 551 : s2 = shl( s1, 1 );
905 :
906 551 : s2 = s_min( s2, WORD32_BITS - 1 );
907 :
908 : /* beta: scaled by s2 */
909 : #ifdef ISSUE_1836_replace_overflow_libcom
910 551 : tmp16 = round_fx_sat( L_shl( tmp, s1 ) );
911 : #else
912 : tmp16 = round_fx_o( L_shl_o( tmp, s1, &Overflow ), &Overflow );
913 : #endif
914 551 : beta = mult_r( tmp16, tmp16 );
915 :
916 : /* scalar3: scaled by s3 */
917 551 : scalar3 = L_sub( L_deposit_l( msPsd[j] ), L_deposit_l( msPsdFirstMoment[j] ) );
918 551 : s3 = norm_l( scalar3 );
919 551 : scalar3 = L_shl( scalar3, s3 );
920 :
921 : /* msPsdFirstMoment: 6Q9 */
922 551 : tmp = L_msu( L_mult( beta, msPsdFirstMoment[j] ), beta, msPsd[j] );
923 551 : msPsdFirstMoment[j] = add( round_fx( L_shr( tmp, s2 ) ), msPsd[j] );
924 551 : move16();
925 :
926 : /* msPsdSecondMoment: 12Q19 */
927 551 : tmp0 = L_shr( Mpy_32_16_r( msPsdSecondMoment[j], beta ), s2 );
928 551 : tmp1 = Mpy_32_32( scalar3, scalar3 );
929 551 : tmp1 = L_shr( L_sub( tmp1, L_shr( Mpy_32_16_r( tmp1, beta ), s2 ) ), sub( shl( s3, 1 ), 32 ) );
930 551 : msPsdSecondMoment[j] = L_add( tmp0, tmp1 );
931 551 : move32();
932 :
933 : /* Compute inverse of amount of degrees of freedom */
934 551 : QeqInv = MSQEQINVMAX;
935 551 : move16();
936 :
937 551 : IF( msNoiseFloor[j] != 0 /*0.0 Q15*/ )
938 : {
939 551 : tmp = L_mult( msNoiseFloor[j], msNoiseFloor[j] );
940 551 : tmp16 = BASOP_Util_Divide3232_uu_1616_Scale( msPsdSecondMoment[j], tmp, &s );
941 : /* consider factor of 2 */
942 551 : s = s_min( s_max( sub( s, 1 ), -( WORD16_BITS - 1 ) ), ( WORD16_BITS - 1 ) );
943 551 : if ( s < 0 )
944 : {
945 494 : QeqInv = shl( tmp16, s );
946 : }
947 551 : QeqInv = s_min( QeqInv, MSQEQINVMAX );
948 : }
949 551 : QeqInvAv = L_add( QeqInvAv, L_mult( QeqInv, psize[j] ) );
950 :
951 : /* Compute bias correction Bmin */
952 551 : tmp0 = Mpy_32_16_1( scalar, QeqInv );
953 551 : tmp1 = L_sub( 1073741824l /*0.5 Q31*/, L_mult( msM_win, QeqInv ) );
954 551 : tmp16 = BASOP_Util_Divide3232_uu_1616_Scale( tmp0, tmp1, &s );
955 551 : msBminWin[j] = L_add( 134217728l /*1.0 Q27*/, L_shl( L_deposit_h( tmp16 ), add( s, 7 - 4 ) ) );
956 551 : move32();
957 :
958 551 : tmp0 = Mpy_32_16_1( scalar2, QeqInv );
959 551 : tmp1 = L_sub( 1073741824l /*0.5 Q31*/, L_mult( msM_subwin, QeqInv ) );
960 551 : tmp16 = BASOP_Util_Divide3232_uu_1616_Scale( tmp0, tmp1, &s );
961 551 : msBminSubWin[j] = L_add( 134217728l /*1.0 Q27*/, L_shl( L_deposit_h( tmp16 ), s ) );
962 551 : move32();
963 : }
964 :
965 9 : inv_totsize = BASOP_Util_Divide3232_uu_1616_Scale( 1, totsize, &s );
966 9 : QeqInvAv = Mpy_32_16_1( QeqInvAv, inv_totsize );
967 9 : QeqInvAv = L_shl( QeqInvAv, s );
968 9 : hFdCngCom->msQeqInvAv[cnt] = QeqInvAv;
969 9 : move32();
970 :
971 : /* New minimum? */
972 :
973 : /* exponent QeqInvAv: CNG_S, exponent MSAV: (4>>1) */
974 9 : s = CNG_S + 2 * MSAV_EXP;
975 9 : move16();
976 9 : BminCorr = Mpy_32_16_1( Sqrt32( QeqInvAv, &s ), MSAV );
977 9 : BminCorr = L_shl( BminCorr, sub( s, 1 ) );
978 :
979 : /* exponent BminCorr: 1 */
980 9 : BminCorr = L_add( BminCorr, 1073741824l /*0.5 Q31*/ );
981 :
982 560 : FOR( j = start; j < stop; j++ )
983 : {
984 : /* exponent scalar: CNG_S+1 */
985 551 : scalar = Mpy_32_16_1( BminCorr, msPsd[j] );
986 :
987 : /* exponent scalar2: CNG_S+1+4 */
988 551 : scalar2 = Mpy_32_32( scalar, msBminWin[j] );
989 :
990 551 : msNewMinFlag[j] = 0;
991 551 : move16();
992 551 : IF( LT_32( scalar2, msCurrentMin[j] ) /*0.0 Q31*/ )
993 : {
994 323 : msNewMinFlag[j] = 1;
995 323 : move16();
996 : /* exponent msCurrentMin[j]: CNG_S+1+4 */
997 323 : msCurrentMin[j] = scalar2;
998 323 : move32();
999 : /* exponent msCurrentMinSubWindow[j]: CNG_S */
1000 : BASOP_SATURATE_WARNING_OFF_EVS;
1001 : #ifdef ISSUE_1836_replace_overflow_libcom
1002 323 : msCurrentMinSubWindow[j] = L_shl_sat( Mpy_32_32( scalar, msBminSubWin[j] ), 5 );
1003 323 : move32();
1004 : #else
1005 : msCurrentMinSubWindow[j] = L_shl_o( Mpy_32_32( scalar, msBminSubWin[j] ), 5, &Overflow );
1006 : move32();
1007 : #endif
1008 : BASOP_SATURATE_WARNING_ON_EVS;
1009 : }
1010 : }
1011 :
1012 : /* This is used later to identify local minima */
1013 9 : IF( GE_16( hFdCngCom->msFrCnt, MSSUBFRLEN ) )
1014 : {
1015 7 : FOR( i = 0; i < 3; i++ )
1016 : {
1017 6 : IF( LT_32( hFdCngCom->msQeqInvAv[cnt], L_shr( L_deposit_h( msQeqInvAv_thresh[i] ), CNG_S ) ) /*0.0 Q31*/ )
1018 : {
1019 1 : BREAK;
1020 : }
1021 : }
1022 : /* format 1Q30 */
1023 2 : hFdCngCom->msSlope[cnt] = msNoiseSlopeMax[i];
1024 2 : move32();
1025 : }
1026 :
1027 : /* Consider the FFT and CLDFB bands separately */
1028 9 : start = stop;
1029 9 : move16();
1030 9 : stop = len;
1031 9 : move16();
1032 9 : totsize = sub( hFdCngCom->stopBand, hFdCngCom->stopFFTbin );
1033 9 : cnt = add( cnt, 1 );
1034 : } /*while (stop > start)*/
1035 :
1036 : /* Update minimum between sub windows */
1037 9 : test();
1038 9 : IF( GT_16( hFdCngCom->msFrCnt, 1 ) && LT_16( hFdCngCom->msFrCnt, MSSUBFRLEN ) )
1039 : {
1040 310 : FOR( j = 0; j < len; j++ )
1041 : {
1042 305 : if ( msNewMinFlag[j] > 0 )
1043 : {
1044 77 : msLocalMinFlag[j] = 1;
1045 77 : move16();
1046 : }
1047 305 : if ( LT_32( msCurrentMinSubWindow[j], msCurrentMinOut[j] ) /*0.0 Q31*/ )
1048 : {
1049 : /* msCurrentMinOut[j] scaled with CNG_S */
1050 88 : msCurrentMinOut[j] = msCurrentMinSubWindow[j];
1051 88 : move32();
1052 : }
1053 : }
1054 : /* Get the current noise floor */
1055 5 : Copy_Scale_sig_32_16( msCurrentMinOut, msNoiseFloor, len, -16 );
1056 : }
1057 : ELSE /* sub window complete */
1058 : {
1059 4 : IF( GE_16( hFdCngCom->msFrCnt, MSSUBFRLEN ) )
1060 : {
1061 : /* Collect buffers */
1062 2 : Copy32( msCurrentMinSubWindow, msMinBuf + len * hFdCngCom->msMinBufferPtr, len );
1063 :
1064 : /* Compute minimum among all buffers */
1065 2 : Copy32( msMinBuf, msCurrentMinOut, len );
1066 2 : ptr = msMinBuf + len;
1067 12 : FOR( i = 1; i < MSNUMSUBFR; i++ )
1068 : {
1069 625 : FOR( j = 0; j < len; j++ )
1070 : {
1071 615 : if ( LT_32( *ptr, msCurrentMinOut[j] ) /*0.0 Q31*/ )
1072 : {
1073 0 : msCurrentMinOut[j] = *ptr;
1074 0 : move32();
1075 : }
1076 615 : ptr++;
1077 : }
1078 : }
1079 :
1080 : /* Take over local minima */
1081 2 : slope = hFdCngCom->msSlope[0];
1082 2 : move16();
1083 125 : FOR( j = 0; j < len; j++ )
1084 : {
1085 123 : if ( EQ_16( j, lenFFT ) )
1086 : {
1087 0 : slope = hFdCngCom->msSlope[1];
1088 0 : move16();
1089 : }
1090 123 : test();
1091 123 : test();
1092 123 : test();
1093 123 : IF( ( msLocalMinFlag[j] != 0 ) && ( msNewMinFlag[j] == 0 ) && ( LT_32( L_shr( msCurrentMinSubWindow[j], 1 ), Mpy_32_16_1( msCurrentMinOut[j], slope ) ) /*0.0 Q31*/ ) && ( GT_32( msCurrentMinSubWindow[j], msCurrentMinOut[j] ) /*0.0 Q31*/ ) )
1094 : {
1095 0 : msCurrentMinOut[j] = msCurrentMinSubWindow[j];
1096 0 : move32();
1097 0 : i = j;
1098 0 : move16();
1099 0 : FOR( k = 0; k < MSNUMSUBFR; k++ )
1100 : {
1101 0 : msMinBuf[i] = msCurrentMinOut[j];
1102 0 : move32();
1103 0 : i = add( i, len );
1104 : }
1105 : }
1106 : }
1107 :
1108 : /* Reset */
1109 2 : set16_fx( msLocalMinFlag, 0, len );
1110 2 : set32_fx( msCurrentMin, 2147483647l /*1.0 Q31*/, len );
1111 :
1112 : /* Get the current noise floor */
1113 2 : Copy_Scale_sig_32_16( msCurrentMinOut, msNoiseFloor, len, -16 );
1114 : }
1115 : }
1116 :
1117 :
1118 : /* Detect sudden offsets based on the FFT bins (core bandwidth) */
1119 9 : IF( GT_32( Mpy_32_16_1( msPsdSum[0], 655 /*0.02 Q15*/ ), msPeriodogSum[0] ) /*0.0 Q31*/ )
1120 : {
1121 0 : IF( hFdCngCom->offsetflag > 0 )
1122 : {
1123 0 : Copy( msPeriodog, msPsd, len );
1124 0 : FOR( j = 0; j < len; j++ )
1125 : {
1126 0 : msCurrentMinOut[j] = L_deposit_h( msPeriodog[j] );
1127 0 : move32();
1128 : }
1129 0 : set32_fx( hFdCngCom->msAlphaCor, 2147483647l /*1.0 Q31*/, cnt );
1130 0 : set32_fx( msAlpha, 0l /*0.0 Q31*/, len );
1131 0 : Copy( msPeriodog, msPsdFirstMoment, len );
1132 0 : set32_fx( msPsdSecondMoment, 0l /*0.0 Q31*/, len );
1133 :
1134 0 : msPsdSum[0] = dotp_s_fx( msPeriodog, psize, lenFFT, CNG_HS );
1135 0 : move32();
1136 0 : IF( LT_16( lenFFT, len ) )
1137 : {
1138 0 : msPsdSum[1] = dotp_s_fx( msPeriodog + lenFFT, psize + lenFFT, sub( len, lenFFT ), CNG_HS );
1139 0 : move32();
1140 : }
1141 : }
1142 0 : hFdCngCom->offsetflag = 1;
1143 0 : move16();
1144 : }
1145 : ELSE
1146 : {
1147 9 : hFdCngCom->offsetflag = 0;
1148 9 : move16();
1149 : }
1150 :
1151 :
1152 : /* Increment frame counter */
1153 9 : IF( EQ_16( hFdCngCom->msFrCnt, MSSUBFRLEN ) )
1154 : {
1155 2 : hFdCngCom->msFrCnt = 1;
1156 2 : move16();
1157 2 : hFdCngCom->msMinBufferPtr = add( hFdCngCom->msMinBufferPtr, 1 );
1158 2 : move16();
1159 2 : if ( EQ_16( hFdCngCom->msMinBufferPtr, MSNUMSUBFR ) )
1160 : {
1161 0 : hFdCngCom->msMinBufferPtr = 0;
1162 0 : move16();
1163 : }
1164 : }
1165 : ELSE
1166 : {
1167 7 : hFdCngCom->msFrCnt = add( hFdCngCom->msFrCnt, 1 );
1168 : }
1169 :
1170 : /* Smooth noise estimate during CNG phases */
1171 560 : FOR( j = 0; j < len; j++ )
1172 : {
1173 551 : msNoiseEst[j] = round_fx( L_mac( L_mult( 31130 /*0.95 Q15*/, msNoiseEst[j] ), 1638 /*0.05 Q15*/, msNoiseFloor[j] ) );
1174 551 : move16();
1175 : }
1176 : }
1177 : /* Collect buffers */
1178 15 : Copy( msPeriodog, msPeriodogBuf + len * ( *msPeriodogBufPtr ), len );
1179 :
1180 15 : *msPeriodogBufPtr = add( *msPeriodogBufPtr, 1 );
1181 15 : move16();
1182 15 : if ( EQ_16( *msPeriodogBufPtr, MSBUFLEN ) )
1183 : {
1184 2 : ( *msPeriodogBufPtr ) = 0;
1185 2 : move16();
1186 : }
1187 :
1188 : /* Upper limit the noise floors with the averaged input energy */
1189 937 : FOR( j = 0; j < len; j++ )
1190 : {
1191 922 : scalar = L_mult( msPeriodogBuf[j], 6554 /*1.0/MSBUFLEN Q15*/ );
1192 :
1193 4610 : FOR( i = j + len; i < MSBUFLEN * len; i += len )
1194 : {
1195 3688 : scalar = L_mac( scalar, msPeriodogBuf[i], 6554 /*1.0/MSBUFLEN Q15*/ );
1196 : }
1197 922 : scalar16 = round_fx( scalar );
1198 922 : if ( GT_16( msNoiseEst[j], scalar16 ) /*0.0 Q15*/ )
1199 : {
1200 361 : msNoiseEst[j] = scalar16;
1201 361 : move16();
1202 : }
1203 :
1204 922 : assert( msNoiseEst[j] >= 0 /*0.0 Q15*/ );
1205 : }
1206 15 : }
1207 :
1208 44303 : void minimum_statistics_fx(
1209 : Word16 len, /* i : Total number of partitions (CLDFB or FFT) */
1210 : Word16 lenFFT, /* i : Number of FFT partitions */
1211 : Word16 *psize, /* i : Partition sizes, fractional Q9*/
1212 : Word16 *msPeriodog, /* i : Periodogram (energies) */
1213 : Word32 *msNoiseFloor, /* i/o: Noise floors (energies) Q25*/
1214 : Word32 *msNoiseEst, /* i/o: Noise estimates (energies) Q25*/
1215 : Word32 *msAlpha, /* i/o: Forgetting factors */
1216 : Word16 *msPsd, /* i/o: Power Spectral Density (smoothed periodogram => energies) */
1217 : Word16 *msPsdFirstMoment, /* i/o: PSD statistics of 1st order (energy means) */
1218 : Word32 *msPsdSecondMoment, /* i/o: PSD statistics of 2nd order (energy variances) */
1219 : Word32 *msMinBuf, /* i/o: Buffer of minima (energies) */
1220 : Word32 *msBminWin, /* o : Bias correction factors */
1221 : Word32 *msBminSubWin, /* o : Bias correction factors */
1222 : Word32 *msCurrentMin, /* i/o: Local minima (energies) */
1223 : Word32 *msCurrentMinOut, /* i/o: Local minima (energies) */
1224 : Word32 *msCurrentMinSubWindow, /* i/o: Local minima (energies) */
1225 : Word16 *msLocalMinFlag, /* i : Binary flag */
1226 : Word16 *msNewMinFlag, /* i : Binary flag */
1227 : Word16 *msPeriodogBuf, /* i/o: Buffer of periodograms (energies) */
1228 : Word16 *msPeriodogBufPtr, /* i/o: Counter */
1229 : HANDLE_FD_CNG_COM hFdCngCom, /* i/o: FD_CNG structure containing all buffers and variables */
1230 : const Word16 enc_dec, /* i : encoder/decoder indicator */
1231 : const Word16 element_mode /* i : IVAS element mode type */
1232 : )
1233 : {
1234 : Word16 i, j, k, s, s1, s2, s3;
1235 : Word16 len2;
1236 : Word16 current_len;
1237 : Word16 start, stop, cnt;
1238 : Word16 totsize;
1239 : Word16 inv_totsize;
1240 :
1241 : Word32 tmp, tmp0, tmp1;
1242 : Word32 scalar, scalar2, scalar3;
1243 : Word32 snr;
1244 : Word32 msAlphaHatMin2;
1245 : Word32 BminCorr;
1246 : Word32 QeqInvAv;
1247 : Word32 *ptr;
1248 : Word32 *msPsdSum;
1249 : Word32 *msPeriodogSum;
1250 :
1251 : Word16 tmp16, tmp16_1;
1252 : Word16 beta;
1253 : Word16 slope;
1254 : Word16 QeqInv;
1255 : Word16 scalar16;
1256 : Word16 scalar216;
1257 : Word16 scalar316;
1258 : Word16 msM_win;
1259 : Word16 msM_subwin;
1260 : Word16 msAlphaCorAlpha;
1261 : Word16 msAlphaCorAlpha2;
1262 : Word16 msPeriodogSum16;
1263 : Word16 msNoiseFloor16;
1264 : #ifndef ISSUE_1836_replace_overflow_libcom
1265 : #ifdef BASOP_NOGLOB_DECLARE_LOCAL
1266 : Flag Overflow = 0;
1267 : move32();
1268 : #endif
1269 : #endif
1270 :
1271 :
1272 44303 : len2 = i_mult( MSNUMSUBFR, len );
1273 :
1274 44303 : msM_win = hFdCngCom->msM_win;
1275 44303 : move16();
1276 44303 : msM_subwin = hFdCngCom->msM_subwin;
1277 44303 : move16();
1278 :
1279 44303 : msAlphaCorAlpha = MSALPHACORALPHA;
1280 44303 : move16();
1281 44303 : msAlphaCorAlpha2 = MSALPHACORALPHA2;
1282 44303 : move16();
1283 :
1284 44303 : msPsdSum = hFdCngCom->msPsdSum;
1285 44303 : msPeriodogSum = hFdCngCom->msPeriodogSum;
1286 :
1287 : /* No minimum statistics at initialization */
1288 44303 : IF( LT_16( hFdCngCom->msFrCnt_init_counter, hFdCngCom->msFrCnt_init_thresh ) )
1289 : {
1290 3919 : Copy( msPeriodog, msPsd, len ); /* 6Q9 */
1291 3919 : Copy_Scale_sig_16_32_no_sat( msPeriodog, msNoiseFloor, len, Q16 ); /* Q25 */
1292 3919 : Copy_Scale_sig_16_32_no_sat( msPeriodog, msNoiseEst, len, Q16 ); /* Q25 */
1293 3919 : Copy( msPeriodog, msPsdFirstMoment, len ); /* 6Q9 */
1294 :
1295 3919 : set32_fx( msPsdSecondMoment, 0l /*0.0 Q31*/, len );
1296 3919 : msPeriodogSum[0] = dotp_s_fx( msPeriodog, psize, lenFFT, CNG_HS );
1297 3919 : move32();
1298 3919 : msPsdSum[0] = msPeriodogSum[0]; /* 16Q15 */
1299 3919 : move32();
1300 :
1301 3919 : IF( LT_16( lenFFT, len ) )
1302 : {
1303 0 : msPeriodogSum[1] = dotp_s_fx( msPeriodog + lenFFT, psize + lenFFT, sub( len, lenFFT ), CNG_HS );
1304 0 : move32();
1305 0 : msPsdSum[1] = msPeriodogSum[1]; /* 16Q15 */
1306 0 : move32();
1307 : }
1308 :
1309 : /* Increment frame counter at initialization */
1310 : /* Some frames are sometimes zero at initialization => ignore them */
1311 3919 : IF( LT_16( msPeriodog[0], hFdCngCom->init_old ) )
1312 : {
1313 710 : set32_fx( msCurrentMinOut, 2147483647l /*1.0 Q31*/, len ); /* 16Q15 */
1314 710 : set32_fx( msCurrentMin, 2147483647l /*1.0 Q31*/, len ); /* 16Q15 */
1315 710 : set32_fx( msMinBuf, 2147483647l /*1.0 Q31*/, len2 ); /* 16Q15 */
1316 710 : set32_fx( msCurrentMinSubWindow, 2147483647l /*1.0 Q31*/, len ); /* 16Q15 */
1317 :
1318 710 : hFdCngCom->msFrCnt_init_counter = add( hFdCngCom->msFrCnt_init_counter, 1 );
1319 710 : move16();
1320 : }
1321 3919 : hFdCngCom->init_old = msPeriodog[0]; /* 6Q9 */
1322 3919 : move16();
1323 : }
1324 : ELSE
1325 : {
1326 : /* Consider the FFT and CLDFB bands separately
1327 : - first iteration for FFT bins,
1328 : - second one for CLDFB bands in SWB mode */
1329 40384 : cnt = 0;
1330 40384 : move16();
1331 40384 : start = 0;
1332 40384 : move16();
1333 40384 : stop = lenFFT;
1334 40384 : move16();
1335 40384 : totsize = sub( hFdCngCom->stopFFTbin, hFdCngCom->startBand );
1336 80768 : WHILE( GT_16( stop, start ) )
1337 : {
1338 40384 : current_len = sub( stop, start );
1339 :
1340 : /* Compute smoothed correction factor for PSD smoothing */
1341 :
1342 : /* msPeriodogSum[cnt] with format 16Q15 */
1343 40384 : msPeriodogSum[cnt] = dotp_s_fx( msPeriodog + start, psize + start, current_len, CNG_HS );
1344 40384 : move32();
1345 :
1346 40384 : IF( msPeriodogSum[cnt] == 0 )
1347 : {
1348 0 : hFdCngCom->msAlphaCor[cnt] = Mpy_32_16_1( hFdCngCom->msAlphaCor[cnt], msAlphaCorAlpha );
1349 0 : move32();
1350 : }
1351 : ELSE
1352 : {
1353 : /* calculate scalar with normalized msPeriodogSum[cnt], exponent -2*s1 */
1354 40384 : s1 = norm_l( msPeriodogSum[cnt] );
1355 40384 : msPeriodogSum16 = round_fx_sat( L_shl( msPeriodogSum[cnt], s1 ) );
1356 40384 : scalar = L_mult( msPeriodogSum16, msPeriodogSum16 );
1357 :
1358 : /* calculate difference, both elements in 16Q15 format, use absolute value
1359 : to avoid -1.0 x -1.0 multiplications later */
1360 40384 : scalar2 = L_abs( L_sub( msPsdSum[cnt], msPeriodogSum[cnt] ) );
1361 :
1362 40384 : s2 = WORD32_BITS - 1;
1363 40384 : move16();
1364 40384 : if ( scalar2 != 0 )
1365 : {
1366 : /* use absolute value to avoid -1.0 x -1.0 multiplications */
1367 40128 : s2 = norm_l( scalar2 );
1368 : }
1369 : #ifdef ISSUE_1836_replace_overflow_libcom
1370 40384 : scalar216 = round_fx_sat( L_shl( scalar2, s2 ) );
1371 : #else
1372 : scalar216 = round_fx_o( L_shl_o( scalar2, s2, &Overflow ), &Overflow );
1373 : #endif
1374 40384 : scalar2 = L_mult( scalar216, scalar216 );
1375 :
1376 : /* find common exponent */
1377 40384 : tmp16_1 = sub( s1, s2 );
1378 40384 : tmp16 = s_min( shl( abs_s( tmp16_1 ), 1 ), WORD32_BITS - 1 );
1379 40384 : if ( tmp16_1 < 0 )
1380 : {
1381 39066 : scalar2 = L_shr( scalar2, tmp16 );
1382 : }
1383 40384 : if ( tmp16_1 > 0 )
1384 : {
1385 550 : scalar = L_shr( scalar, tmp16 );
1386 : }
1387 :
1388 :
1389 : /* add scalar and scalar2, avoid overflows */
1390 40384 : scalar = L_shr( scalar, 1 );
1391 40384 : scalar2 = L_shr( scalar2, 1 );
1392 40384 : scalar3 = L_add( scalar, scalar2 );
1393 :
1394 : /* calculate division */
1395 40384 : scalar16 = BASOP_Util_Divide3232_uu_1616_Scale( scalar, scalar3, &s3 );
1396 40384 : s3 = s_max( s3, -( WORD16_BITS - 1 ) );
1397 40384 : scalar16 = shl( scalar16, s3 );
1398 40384 : scalar16 = s_max( scalar16, MSALPHACORMAX );
1399 :
1400 40384 : hFdCngCom->msAlphaCor[cnt] = Madd_32_16( L_mult( scalar16, msAlphaCorAlpha2 ), hFdCngCom->msAlphaCor[cnt], msAlphaCorAlpha );
1401 40384 : move32();
1402 : }
1403 :
1404 : /* Compute SNR */
1405 :
1406 : /* msPeriodogSum[cnt] with format 16Q15 */
1407 : Word16 msNoiseFloor16tmp[NPART_SHAPING];
1408 40384 : IF( enc_dec == ENC )
1409 : {
1410 0 : Copy_Scale_sig32_16( msNoiseFloor, msNoiseFloor16tmp, NPART, 0 );
1411 : }
1412 : ELSE
1413 : {
1414 40384 : Copy_Scale_sig32_16( msNoiseFloor, msNoiseFloor16tmp, NPART_SHAPING, 0 );
1415 : }
1416 40384 : snr = dotp_s_fx( &msNoiseFloor16tmp[start], psize + start, current_len, CNG_HS );
1417 :
1418 40384 : IF( GT_32( L_shr( Mpy_32_16_1( msPsdSum[cnt], 18431 /*0.56246299817 Q15*/ ), 13 ), snr ) )
1419 : {
1420 0 : tmp0 = BASOP_Util_Log2( msPsdSum[cnt] );
1421 0 : tmp1 = BASOP_Util_Log2( snr );
1422 0 : tmp1 = L_sub( tmp0, tmp1 );
1423 0 : tmp1 = Mpy_32_16_1( tmp1, MSSNREXP );
1424 0 : msAlphaHatMin2 = BASOP_Util_InvLog2( tmp1 );
1425 : }
1426 : ELSE
1427 : {
1428 40384 : msAlphaHatMin2 = MSALPHAHATMIN;
1429 40384 : move32();
1430 : }
1431 40384 : scalar = Mpy_32_16_1( hFdCngCom->msAlphaCor[cnt], MSALPHAMAX );
1432 :
1433 2521709 : FOR( j = start; j < stop; j++ )
1434 : {
1435 : /* Compute optimal smoothing parameter for PSD estimation */
1436 2481325 : test();
1437 2481325 : IF( ( scalar == 0 ) || ( msNoiseFloor[j] == 0 ) )
1438 : {
1439 81 : msAlpha[j] = msAlphaHatMin2;
1440 81 : move32();
1441 : }
1442 : ELSE
1443 : {
1444 : /* calculate scalar2 with normalized msNoiseFloor[j], exponent -2*s1 */
1445 2481244 : s1 = WORD16_BITS - 1;
1446 2481244 : move16();
1447 2481244 : if ( msNoiseFloor16tmp[j] != 0 )
1448 : {
1449 2481244 : s1 = norm_s( msNoiseFloor16tmp[j] );
1450 : }
1451 2481244 : msNoiseFloor16 = shl( msNoiseFloor16tmp[j], s1 );
1452 2481244 : scalar2 = L_mult( msNoiseFloor16, msNoiseFloor16 );
1453 :
1454 : /* calculate difference, both elements in 6Q9 format, use absolute value
1455 : to avoid -1.0 x -1.0 multiplications later */
1456 2481244 : scalar316 = abs_s( sub( msPsd[j], msNoiseFloor16tmp[j] ) );
1457 :
1458 2481244 : s2 = WORD16_BITS - 1;
1459 2481244 : move16();
1460 2481244 : if ( scalar316 != 0 )
1461 : {
1462 : /* use absolute value to avoid -1.0 x -1.0 multiplications */
1463 2036607 : s2 = norm_s( scalar316 );
1464 : }
1465 2481244 : scalar316 = shl( scalar316, s2 );
1466 2481244 : scalar3 = L_mult( scalar316, scalar316 );
1467 :
1468 : /* find common exponent */
1469 2481244 : tmp16_1 = sub( s1, s2 );
1470 2481244 : tmp16 = s_min( shl( abs_s( tmp16_1 ), 1 ), WORD32_BITS - 1 );
1471 2481244 : if ( tmp16_1 < 0 )
1472 : {
1473 2176242 : scalar3 = L_shr( scalar3, tmp16 );
1474 : }
1475 2481244 : if ( tmp16_1 > 0 )
1476 : {
1477 184976 : scalar2 = L_shr( scalar2, tmp16 );
1478 : }
1479 :
1480 : /* add scalar2 and scalar3, avoid overflows */
1481 2481244 : scalar2 = L_shr( scalar2, 1 );
1482 2481244 : scalar3 = L_shr( scalar3, 1 );
1483 2481244 : scalar3 = L_add( scalar2, scalar3 );
1484 :
1485 : /* calculate division */
1486 2481244 : tmp16 = BASOP_Util_Divide3232_uu_1616_Scale( scalar2, scalar3, &s3 );
1487 2481244 : scalar2 = Mpy_32_16_1( scalar, tmp16 );
1488 2481244 : s3 = s_max( s3, -( WORD32_BITS - 1 ) );
1489 2481244 : scalar2 = L_shl( scalar2, s3 );
1490 2481244 : msAlpha[j] = L_max( scalar2, msAlphaHatMin2 );
1491 2481244 : move32();
1492 : }
1493 :
1494 : /* Compute the PSD (smoothed periodogram) in each band */
1495 2481325 : msPsd[j] = round_fx( Madd_32_16( Mpy_32_16_1( msAlpha[j], msPsd[j] ), L_sub( 2147483647l /*1.0 Q31*/, msAlpha[j] ), msPeriodog[j] ) );
1496 2481325 : move16();
1497 : }
1498 40384 : msPsdSum[cnt] = dotp_s_fx( msPsd + start, psize + start, current_len, CNG_HS );
1499 40384 : move32();
1500 :
1501 40384 : QeqInvAv = 0l /*0.0 Q31*/;
1502 40384 : move32();
1503 :
1504 : /* scalar: 7Q24 format */
1505 40384 : tmp = 1191182336l /*(float)(MSNUMSUBFR*MSSUBFRLEN)-1.0 Q24*/;
1506 40384 : move32();
1507 40384 : scalar = L_sub( tmp, L_mult( round_fx( tmp ), msM_win ) );
1508 :
1509 : /* scalar2: 4Q27 format */
1510 40384 : tmp = 1476395008l /*(float)MSSUBFRLEN-1.0 Q27*/;
1511 40384 : move32();
1512 40384 : scalar2 = L_sub( tmp, L_mult( round_fx( tmp ), msM_subwin ) );
1513 :
1514 2521709 : FOR( j = start; j < stop; j++ )
1515 : {
1516 : /* Compute variance of PSD */
1517 2481325 : tmp = L_min( msAlpha[j], MSBETAMAX_SQRT );
1518 :
1519 2481325 : s1 = WORD32_BITS - 1;
1520 2481325 : move16();
1521 2481325 : if ( tmp != 0 )
1522 : {
1523 2481325 : s1 = norm_l( tmp );
1524 : }
1525 2481325 : s2 = shl( s1, 1 );
1526 :
1527 2481325 : s2 = s_min( s2, WORD32_BITS - 1 );
1528 :
1529 : /* beta: scaled by s2 */
1530 : #ifdef ISSUE_1836_replace_overflow_libcom
1531 2481325 : tmp16 = round_fx_sat( L_shl( tmp, s1 ) );
1532 : #else
1533 : tmp16 = round_fx_o( L_shl_o( tmp, s1, &Overflow ), &Overflow );
1534 : #endif
1535 2481325 : beta = mult_r( tmp16, tmp16 );
1536 :
1537 : /* scalar3: scaled by s3 */
1538 2481325 : scalar3 = L_sub( L_deposit_l( msPsd[j] ), L_deposit_l( msPsdFirstMoment[j] ) );
1539 2481325 : s3 = norm_l( scalar3 );
1540 2481325 : scalar3 = L_shl( scalar3, s3 );
1541 :
1542 : /* msPsdFirstMoment: 6Q9 */
1543 2481325 : tmp = L_msu( L_mult( beta, msPsdFirstMoment[j] ), beta, msPsd[j] );
1544 2481325 : msPsdFirstMoment[j] = add( round_fx( L_shr( tmp, s2 ) ), msPsd[j] );
1545 2481325 : move16();
1546 :
1547 : /* msPsdSecondMoment: 12Q19 */
1548 2481325 : tmp0 = L_shr( Mpy_32_16_r( msPsdSecondMoment[j], beta ), s2 );
1549 2481325 : tmp1 = Mpy_32_32( scalar3, scalar3 );
1550 2481325 : tmp1 = L_shr( L_sub( tmp1, L_shr( Mpy_32_16_r( tmp1, beta ), s2 ) ), sub( shl( s3, 1 ), 32 ) );
1551 2481325 : msPsdSecondMoment[j] = L_add( tmp0, tmp1 );
1552 2481325 : move32();
1553 :
1554 : /* Compute inverse of amount of degrees of freedom */
1555 2481325 : QeqInv = MSQEQINVMAX;
1556 2481325 : move16();
1557 :
1558 2481325 : IF( msNoiseFloor[j] != 0 /*0.0 Q15*/ )
1559 : {
1560 2481244 : tmp = Mpy_32_32( msNoiseFloor[j], msNoiseFloor[j] );
1561 2481244 : tmp16 = BASOP_Util_Divide3232_uu_1616_Scale( msPsdSecondMoment[j], tmp, &s );
1562 : /* consider factor of 2 */
1563 2481244 : s = s_min( s_max( sub( s, 1 ), -( WORD16_BITS - 1 ) ), ( WORD16_BITS - 1 ) );
1564 2481244 : if ( s < 0 )
1565 : {
1566 2301446 : QeqInv = shl( tmp16, s );
1567 : }
1568 2481244 : QeqInv = s_min( QeqInv, MSQEQINVMAX );
1569 : }
1570 2481325 : QeqInvAv = L_add( QeqInvAv, L_mult( QeqInv, psize[j] ) );
1571 :
1572 : /* Compute bias correction Bmin */
1573 2481325 : tmp0 = Mpy_32_16_1( scalar, QeqInv );
1574 2481325 : tmp1 = L_msu( 1073741824l /*0.5 Q31*/, msM_win, QeqInv );
1575 2481325 : tmp16 = BASOP_Util_Divide3232_uu_1616_Scale( tmp0, tmp1, &s );
1576 2481325 : msBminWin[j] = L_add( 134217728l /*1.0 Q27*/, L_shl( L_deposit_h( tmp16 ), add( s, 7 - 4 ) ) );
1577 2481325 : move32();
1578 :
1579 2481325 : tmp0 = Mpy_32_16_1( scalar2, QeqInv );
1580 2481325 : tmp1 = L_msu( 1073741824l /*0.5 Q31*/, msM_subwin, QeqInv );
1581 2481325 : tmp16 = BASOP_Util_Divide3232_uu_1616_Scale( tmp0, tmp1, &s );
1582 2481325 : msBminSubWin[j] = L_add( 134217728l /*1.0 Q27*/, L_shl( L_deposit_h( tmp16 ), s ) );
1583 2481325 : move32();
1584 : }
1585 :
1586 40384 : inv_totsize = BASOP_Util_Divide3232_uu_1616_Scale( 1, totsize, &s );
1587 40384 : QeqInvAv = Mpy_32_16_1( QeqInvAv, inv_totsize );
1588 40384 : QeqInvAv = L_shl( QeqInvAv, s );
1589 40384 : hFdCngCom->msQeqInvAv[cnt] = QeqInvAv;
1590 40384 : move32();
1591 :
1592 : /* New minimum? */
1593 :
1594 : /* exponent QeqInvAv: CNG_S, exponent MSAV: (4>>1) */
1595 40384 : s = CNG_S + 2 * MSAV_EXP;
1596 40384 : move16();
1597 40384 : BminCorr = Mpy_32_16_1( Sqrt32( QeqInvAv, &s ), MSAV );
1598 40384 : BminCorr = L_shl( BminCorr, sub( s, 1 ) );
1599 :
1600 : /* exponent BminCorr: 1 */
1601 40384 : BminCorr = L_add( BminCorr, 1073741824l /*0.5 Q31*/ );
1602 :
1603 2521709 : FOR( j = start; j < stop; j++ )
1604 : {
1605 : /* exponent scalar: CNG_S+1 */
1606 2481325 : scalar = Mpy_32_16_1( BminCorr, msPsd[j] );
1607 :
1608 : /* exponent scalar2: CNG_S+1+4 */
1609 2481325 : scalar2 = Mpy_32_32( scalar, msBminWin[j] );
1610 :
1611 2481325 : msNewMinFlag[j] = 0;
1612 2481325 : move16();
1613 2481325 : IF( LT_32( scalar2, msCurrentMin[j] ) /*0.0 Q31*/ )
1614 : {
1615 873470 : msNewMinFlag[j] = 1;
1616 873470 : move16();
1617 : /* exponent msCurrentMin[j]: CNG_S+1+4 */
1618 873470 : msCurrentMin[j] = scalar2;
1619 873470 : move32();
1620 : /* exponent msCurrentMinSubWindow[j]: CNG_S */
1621 : BASOP_SATURATE_WARNING_OFF_EVS;
1622 : #ifdef ISSUE_1836_replace_overflow_libcom
1623 873470 : msCurrentMinSubWindow[j] = L_shl_sat( Mpy_32_32( scalar, msBminSubWin[j] ), 5 );
1624 : #else
1625 : msCurrentMinSubWindow[j] = L_shl_o( Mpy_32_32( scalar, msBminSubWin[j] ), 5, &Overflow );
1626 : #endif
1627 873470 : move32();
1628 : BASOP_SATURATE_WARNING_ON_EVS;
1629 : }
1630 : }
1631 :
1632 : /* This is used later to identify local minima */
1633 40384 : IF( GE_16( hFdCngCom->msFrCnt, MSSUBFRLEN ) )
1634 : {
1635 9691 : FOR( i = 0; i < 3; i++ )
1636 : {
1637 8325 : IF( LT_32( hFdCngCom->msQeqInvAv[cnt], L_shr( L_deposit_h( msQeqInvAv_thresh[i] ), CNG_S ) ) /*0.0 Q31*/ )
1638 : {
1639 2360 : BREAK;
1640 : }
1641 : }
1642 : /* format 1Q30 */
1643 3726 : hFdCngCom->msSlope[cnt] = msNoiseSlopeMax[i];
1644 3726 : move32();
1645 : }
1646 :
1647 : /* Consider the FFT and CLDFB bands separately */
1648 40384 : start = stop;
1649 40384 : move16();
1650 40384 : stop = len;
1651 40384 : move16();
1652 40384 : totsize = sub( hFdCngCom->stopBand, hFdCngCom->stopFFTbin );
1653 40384 : cnt = add( cnt, 1 );
1654 : } /*while (stop > start)*/
1655 :
1656 : /* Update minimum between sub windows */
1657 40384 : test();
1658 40384 : IF( GT_16( hFdCngCom->msFrCnt, 1 ) && LT_16( hFdCngCom->msFrCnt, MSSUBFRLEN ) )
1659 : {
1660 2061564 : FOR( j = 0; j < len; j++ )
1661 : {
1662 2028549 : if ( msNewMinFlag[j] > 0 )
1663 : {
1664 576430 : msLocalMinFlag[j] = 1;
1665 576430 : move16();
1666 : }
1667 2028549 : if ( LT_32( msCurrentMinSubWindow[j], msCurrentMinOut[j] ) /*0.0 Q31*/ )
1668 : {
1669 : /* msCurrentMinOut[j] scaled with CNG_S */
1670 185799 : msCurrentMinOut[j] = msCurrentMinSubWindow[j];
1671 185799 : move32();
1672 : }
1673 : }
1674 : /* Get the current noise floor */
1675 33015 : Copy32( msCurrentMinOut, msNoiseFloor, len );
1676 : }
1677 : ELSE /* sub window complete */
1678 : {
1679 7369 : IF( GE_16( hFdCngCom->msFrCnt, MSSUBFRLEN ) )
1680 : {
1681 : /* Collect buffers */
1682 3726 : Copy32( msCurrentMinSubWindow, msMinBuf + len * hFdCngCom->msMinBufferPtr, len );
1683 :
1684 : /* Compute minimum among all buffers */
1685 3726 : Copy32( msMinBuf, msCurrentMinOut, len );
1686 3726 : ptr = msMinBuf + len;
1687 22356 : FOR( i = 1; i < MSNUMSUBFR; i++ )
1688 : {
1689 1163280 : FOR( j = 0; j < len; j++ )
1690 : {
1691 1144650 : if ( LT_32( *ptr, msCurrentMinOut[j] ) /*0.0 Q31*/ )
1692 : {
1693 220175 : msCurrentMinOut[j] = *ptr;
1694 220175 : move32();
1695 : }
1696 1144650 : ptr++;
1697 : }
1698 : }
1699 :
1700 : /* Take over local minima */
1701 3726 : slope = hFdCngCom->msSlope[0];
1702 3726 : move16();
1703 232656 : FOR( j = 0; j < len; j++ )
1704 : {
1705 228930 : if ( EQ_16( j, lenFFT ) )
1706 : {
1707 0 : slope = hFdCngCom->msSlope[1];
1708 0 : move16();
1709 : }
1710 228930 : test();
1711 228930 : test();
1712 228930 : test();
1713 228930 : IF( ( msLocalMinFlag[j] != 0 ) && ( msNewMinFlag[j] == 0 ) && ( LT_32( L_shr( msCurrentMinSubWindow[j], 1 ), Mpy_32_16_1( msCurrentMinOut[j], slope ) ) /*0.0 Q31*/ ) && ( GT_32( msCurrentMinSubWindow[j], msCurrentMinOut[j] ) /*0.0 Q31*/ ) )
1714 : {
1715 27545 : msCurrentMinOut[j] = msCurrentMinSubWindow[j];
1716 27545 : move32();
1717 27545 : i = j;
1718 27545 : move16();
1719 192815 : FOR( k = 0; k < MSNUMSUBFR; k++ )
1720 : {
1721 165270 : msMinBuf[i] = msCurrentMinOut[j];
1722 165270 : move32();
1723 165270 : i = add( i, len );
1724 : }
1725 : }
1726 : }
1727 :
1728 : /* Reset */
1729 3726 : set16_fx( msLocalMinFlag, 0, len );
1730 3726 : set32_fx( msCurrentMin, 2147483647l /*1.0 Q31*/, len );
1731 :
1732 : /* Get the current noise floor */
1733 3726 : Copy32( msCurrentMinOut, msNoiseFloor, len );
1734 : }
1735 : }
1736 :
1737 :
1738 : /* Detect sudden offsets based on the FFT bins (core bandwidth) */
1739 40384 : IF( GT_32( Mpy_32_16_1( msPsdSum[0], 655 /*0.02 Q15*/ ), msPeriodogSum[0] ) /*0.0 Q31*/ )
1740 : {
1741 131 : IF( hFdCngCom->offsetflag > 0 )
1742 : {
1743 59 : Copy( msPeriodog, msPsd, len );
1744 3676 : FOR( j = 0; j < len; j++ )
1745 : {
1746 3617 : msCurrentMinOut[j] = L_deposit_h( msPeriodog[j] );
1747 3617 : move32();
1748 : }
1749 59 : set32_fx( hFdCngCom->msAlphaCor, 2147483647l /*1.0 Q31*/, cnt );
1750 59 : set32_fx( msAlpha, 0l /*0.0 Q31*/, len );
1751 59 : Copy( msPeriodog, msPsdFirstMoment, len );
1752 59 : set32_fx( msPsdSecondMoment, 0l /*0.0 Q31*/, len );
1753 :
1754 59 : msPsdSum[0] = dotp_s_fx( msPeriodog, psize, lenFFT, CNG_HS );
1755 59 : move32();
1756 59 : IF( LT_16( lenFFT, len ) )
1757 : {
1758 0 : msPsdSum[1] = dotp_s_fx( msPeriodog + lenFFT, psize + lenFFT, sub( len, lenFFT ), CNG_HS );
1759 0 : move32();
1760 : }
1761 : }
1762 131 : hFdCngCom->offsetflag = 1;
1763 131 : move16();
1764 : }
1765 : ELSE
1766 : {
1767 40253 : hFdCngCom->offsetflag = 0;
1768 40253 : move16();
1769 : }
1770 :
1771 :
1772 : /* Increment frame counter */
1773 40384 : IF( EQ_16( hFdCngCom->msFrCnt, MSSUBFRLEN ) )
1774 : {
1775 3726 : hFdCngCom->msFrCnt = 1;
1776 3726 : move16();
1777 3726 : hFdCngCom->msMinBufferPtr = add( hFdCngCom->msMinBufferPtr, 1 );
1778 3726 : move16();
1779 3726 : if ( EQ_16( hFdCngCom->msMinBufferPtr, MSNUMSUBFR ) )
1780 : {
1781 419 : hFdCngCom->msMinBufferPtr = 0;
1782 419 : move16();
1783 : }
1784 : }
1785 : ELSE
1786 : {
1787 36658 : hFdCngCom->msFrCnt = add( hFdCngCom->msFrCnt, 1 );
1788 36658 : move16();
1789 : }
1790 : Word64 tmp64;
1791 : /* Smooth noise estimate during CNG phases */
1792 2521709 : FOR( j = 0; j < len; j++ )
1793 : {
1794 2481325 : tmp64 = W_add( W_mult_32_16( msNoiseEst[j], 31130 /*0.95 Q15*/ ), W_mult_32_16( msNoiseFloor[j], 1638 /*0.05 Q15*/ ) );
1795 2481325 : msNoiseEst[j] = W_extract_h( W_shl( tmp64, 16 ) ); // Q25
1796 2481325 : move32();
1797 : }
1798 : }
1799 44303 : IF( EQ_16( enc_dec, DEC ) && EQ_16( element_mode, IVAS_CPE_TD ) )
1800 : {
1801 : // v_multc(msNoiseEst, 1.4125f, msNoiseEst, NPART_SHAPING);
1802 0 : v_multc_att32( msNoiseEst, 23142, msNoiseEst, NPART_SHAPING );
1803 : }
1804 : /* Collect buffers */
1805 44303 : Copy( msPeriodog, msPeriodogBuf + len * ( *msPeriodogBufPtr ), len );
1806 :
1807 44303 : *msPeriodogBufPtr = add( *msPeriodogBufPtr, 1 );
1808 44303 : move16();
1809 44303 : if ( EQ_16( *msPeriodogBufPtr, MSBUFLEN ) )
1810 : {
1811 8487 : ( *msPeriodogBufPtr ) = 0;
1812 8487 : move16();
1813 : }
1814 :
1815 : /* Upper limit the noise floors with the averaged input energy */
1816 2766071 : FOR( j = 0; j < len; j++ )
1817 : {
1818 2721768 : scalar = L_mult( msPeriodogBuf[j], 6554 /*1.0/MSBUFLEN Q15*/ );
1819 :
1820 13608840 : FOR( i = j + len; i < MSBUFLEN * len; i += len )
1821 : {
1822 10887072 : scalar = L_mac( scalar, msPeriodogBuf[i], 6554 /*1.0/MSBUFLEN Q15*/ );
1823 : }
1824 2721768 : if ( GT_32( msNoiseEst[j], scalar ) /*0.0 Q15*/ )
1825 : {
1826 706186 : msNoiseEst[j] = scalar;
1827 706186 : move32();
1828 : }
1829 :
1830 2721768 : assert( msNoiseEst[j] >= 0 /*0.0 Q15*/ );
1831 : }
1832 44303 : }
1833 : /*-------------------------------------------------------------------
1834 : * apply_scale()
1835 : *
1836 : * Apply bitrate-dependent scale
1837 : *-------------------------------------------------------------------*/
1838 0 : void apply_scale(
1839 : Word32 *scale, /* o : scalefactor */
1840 : const Word16 bwmode, /* i : audio bandwidth */
1841 : const Word32 bitrate, /* i : Bit rate */
1842 : const SCALE_SETUP *scaleTable, /* i : Scale table Q7 */
1843 : const Word16 scaleTableSize /* i : Size of scale table */
1844 : )
1845 : {
1846 : Word16 i;
1847 : // PMT("Verifiy if the basop ued are ok for stereo too")
1848 0 : FOR( i = 0; i < scaleTableSize; i++ )
1849 : {
1850 0 : cast16();
1851 0 : IF( s_and( (Word16) EQ_16( bwmode, (Word16) scaleTable[i].bwmode ),
1852 : s_and( L_sub( bitrate, scaleTable[i].bitrateFrom ) >= 0,
1853 : L_sub( bitrate, scaleTable[i].bitrateTo ) < 0 ) ) )
1854 : {
1855 0 : BREAK;
1856 : }
1857 : }
1858 :
1859 : {
1860 0 : *scale = L_add( *scale, L_deposit_h( scaleTable[i].scale ) );
1861 0 : move32();
1862 : }
1863 0 : }
1864 :
1865 : /*-------------------------------------------------------------------
1866 : * apply_scale_ind()
1867 : *
1868 : * Apply bitrate-dependent scale
1869 : * returns index of scaleTable
1870 : *-------------------------------------------------------------------*/
1871 0 : Word16 apply_scale_ind(
1872 : Word32 *scale, /* o : scalefactor */
1873 : const Word16 bwmode, /* i : audio bandwidth */
1874 : const Word32 bitrate, /* i : Bit rate */
1875 : const SCALE_SETUP *scaleTable, /* i : Scale table */
1876 : const Word16 scaleTableSize /* i : Size of scale table */
1877 : )
1878 : {
1879 : Word16 i;
1880 : // PMT("Verifiy if the basop ued are ok for stereo too")
1881 0 : FOR( i = 0; i < scaleTableSize; i++ )
1882 : {
1883 0 : cast16();
1884 0 : IF( s_and( (Word16) EQ_16( bwmode, (Word16) scaleTable[i].bwmode ),
1885 : s_and( L_sub( bitrate, scaleTable[i].bitrateFrom ) >= 0,
1886 : L_sub( bitrate, scaleTable[i].bitrateTo ) < 0 ) ) )
1887 : {
1888 0 : BREAK;
1889 : }
1890 : }
1891 :
1892 : {
1893 0 : *scale = L_add( *scale, L_deposit_h( scaleTable[i].scale ) );
1894 0 : move32();
1895 : }
1896 0 : return i;
1897 : }
1898 :
1899 169445 : void apply_scale_ivas_fx(
1900 : Word32 *scale, /* o : scalefactor */
1901 : const Word16 bwmode, /* i : audio bandwidth */
1902 : const Word32 bitrate, /* i : Bit rate */
1903 : const SCALE_SETUP *scaleTable, /* i : Scale table Q7 */
1904 : const Word16 scaleTableSize, /* i : Size of scale table */
1905 : Word16 *index )
1906 : {
1907 : Word16 i;
1908 : // PMT("Verifiy if the basop ued are ok for stereo too")
1909 1789847 : FOR( i = 0; i < scaleTableSize; i++ )
1910 : {
1911 1789847 : cast16();
1912 1789847 : IF( s_and( (Word16) EQ_16( bwmode, (Word16) scaleTable[i].bwmode ),
1913 : s_and( L_sub( bitrate, scaleTable[i].bitrateFrom ) >= 0,
1914 : L_sub( bitrate, scaleTable[i].bitrateTo ) < 0 ) ) )
1915 : {
1916 169445 : BREAK;
1917 : }
1918 : }
1919 169445 : assert( i < scaleTableSize );
1920 :
1921 169445 : *scale = L_add( *scale, L_deposit_h( scaleTable[i].scale ) );
1922 169445 : move32();
1923 169445 : *index = i;
1924 169445 : move16();
1925 169445 : }
1926 : /*-------------------------------------------------------------------
1927 : * bandcombinepow()
1928 : *
1929 : * Compute the power for each partition
1930 : *-------------------------------------------------------------------*/
1931 53576 : void bandcombinepow(
1932 : const Word32 *bandpow, /* i : Power for each band */
1933 : const Word16 exp_bandpow, /* i : exponent of bandpow */
1934 : const Word16 nband, /* i : Number of bands */
1935 : Word16 *part, /* i : Partition upper boundaries (band indices starting from 0) */
1936 : const Word16 npart, /* i : Number of partitions */
1937 : const Word16 *psize_inv, /* i : Inverse partition sizes */
1938 : Word32 *partpow, /* o : Power for each partition */
1939 : Word16 *exp_partpow )
1940 : {
1941 :
1942 : Word16 i, p;
1943 : Word32 temp;
1944 : Word16 smin, len, prev_part;
1945 : Word16 facTabExp[NPART_SHAPING];
1946 :
1947 :
1948 53576 : IF( EQ_16( nband, npart ) )
1949 : {
1950 0 : Copy32( bandpow, partpow, nband ); // Q(31 - exp_bandpow)
1951 0 : smin = 0;
1952 0 : move16();
1953 : }
1954 : ELSE
1955 : {
1956 : /* Compute the power in each partition */
1957 53576 : prev_part = -1;
1958 53576 : move16();
1959 3189996 : FOR( p = 0; p < npart; p++ )
1960 : {
1961 3136420 : len = sub( part[p], prev_part );
1962 3136420 : facTabExp[p] = getScaleFactor32( &bandpow[prev_part + 1], len );
1963 3136420 : move16();
1964 3136420 : prev_part = part[p];
1965 3136420 : move16();
1966 : }
1967 :
1968 53576 : smin = WORD32_BITS - 1;
1969 53576 : move16();
1970 3189996 : FOR( p = 0; p < npart; p++ )
1971 : {
1972 3136420 : smin = s_min( smin, facTabExp[p] );
1973 : }
1974 :
1975 53576 : i = 0;
1976 53576 : move16();
1977 3189996 : FOR( p = 0; p < npart; p++ )
1978 : {
1979 : /* Arithmetic averaging of power for all bins in partition */
1980 3136420 : temp = 0;
1981 3136420 : move32();
1982 18653844 : FOR( ; i <= part[p]; i++ )
1983 : {
1984 15517424 : temp = L_add( temp, Mpy_32_16_1( L_shl( bandpow[i], facTabExp[p] ), psize_inv[p] ) );
1985 : }
1986 3136420 : partpow[p] = L_shr( temp, sub( facTabExp[p], smin ) );
1987 3136420 : move32();
1988 : }
1989 : }
1990 :
1991 53576 : *exp_partpow = sub( exp_bandpow, smin );
1992 53576 : move16();
1993 53576 : }
1994 :
1995 : /*-------------------------------------------------------------------
1996 : * scalebands()
1997 : *
1998 : * Scale partitions (with smoothing)
1999 : *-------------------------------------------------------------------*/
2000 3082 : void scalebands(
2001 : const Word32 *partpow, /* i : Power for each partition Qx*/
2002 : Word16 *part, /* i : Partition upper boundaries (band indices starting from 0) */
2003 : const Word16 npart, /* i : Number of partitions */
2004 : Word16 *midband, /* i : Central band of each partition */
2005 : const Word16 nFFTpart, /* i : Number of FFT partitions */
2006 : const Word16 nband, /* i : Number of bands */
2007 : Word32 *bandpow, /* o : Power for each band Qx*/
2008 : const Word16 flag_fft_en )
2009 : {
2010 : Word16 i, j, s, s1, nint, delta, delta_cmp, delta_s;
2011 : Word16 startBand, startPart, stopPart, stopPartM1;
2012 : Word32 tmp, val, partpowLD64, partpowLD64M1;
2013 :
2014 :
2015 3082 : j = 0;
2016 3082 : move16();
2017 3082 : delta = 0;
2018 3082 : move16();
2019 3082 : partpowLD64M1 = 0L; /* to avoid compilation warnings */
2020 3082 : move32();
2021 :
2022 : /* Interpolate the bin/band-wise levels from the partition levels */
2023 3082 : IF( EQ_16( nband, npart ) )
2024 : {
2025 0 : Copy32( partpow, bandpow, npart );
2026 : }
2027 : ELSE
2028 : {
2029 3082 : startBand = 0;
2030 3082 : move16();
2031 3082 : startPart = 0;
2032 3082 : move16();
2033 3082 : stopPart = nFFTpart;
2034 3082 : move16();
2035 :
2036 8624 : WHILE( startBand < nband )
2037 : {
2038 5542 : stopPartM1 = sub( stopPart, 1 );
2039 5542 : test();
2040 5542 : IF( ( flag_fft_en != 0 ) || ( GE_16( startPart, nFFTpart ) ) )
2041 : {
2042 : /* first half partition */
2043 5542 : j = startPart;
2044 5542 : move16();
2045 :
2046 21063 : FOR( i = startBand; i <= midband[j]; i++ )
2047 : {
2048 15521 : bandpow[i] = partpow[j];
2049 15521 : move32();
2050 : }
2051 5542 : j = add( j, 1 );
2052 :
2053 : /* inner partitions */
2054 5542 : IF( LT_16( j, stopPart ) )
2055 : {
2056 5394 : partpowLD64M1 = BASOP_Util_Log2( partpow[j - 1] );
2057 : }
2058 :
2059 : /* Debug values to check this variable is set. */
2060 5542 : delta = 0x4000; // 1.Q14
2061 5542 : move16();
2062 5542 : delta_cmp = 0x4000; // 1.Q14
2063 5542 : move16();
2064 5542 : s1 = 1;
2065 5542 : move16();
2066 5542 : s = 1;
2067 5542 : move16();
2068 :
2069 72263 : FOR( ; j < stopPart; j++ )
2070 : {
2071 66721 : nint = sub( midband[j], midband[j - 1] );
2072 :
2073 : /* log-linear interpolation */
2074 66721 : partpowLD64 = BASOP_Util_Log2( partpow[j] );
2075 66721 : tmp = L_sub( partpowLD64, partpowLD64M1 );
2076 66721 : tmp = Mpy_32_16_1( tmp, getNormReciprocalWord16( nint ) );
2077 :
2078 : /* scale logarithmic value */
2079 66721 : tmp = L_sub( tmp, DELTA_SHIFT_LD64 );
2080 66721 : delta_s = DELTA_SHIFT;
2081 66721 : move16();
2082 :
2083 66748 : WHILE( tmp > 0 )
2084 : {
2085 27 : tmp = L_sub( tmp, 33554432l /*0.015625 Q31*/ ); // Q31
2086 27 : delta_s = add( delta_s, 1 );
2087 : }
2088 66721 : delta_cmp = shl( 1, s_max( -15, sub( WORD16_BITS - 1, delta_s ) ) );
2089 :
2090 66721 : tmp = BASOP_Util_InvLog2( tmp );
2091 66721 : s = norm_l( tmp );
2092 66721 : s1 = sub( delta_s, s );
2093 :
2094 66721 : delta = round_fx_sat( L_shl_sat( tmp, s ) ); // Q(14+s)
2095 : /* Choose scale such that the interpolation start and end point both are representable and add 1 additional bit hr. */
2096 66721 : delta_s = sub( s_min( norm_l( partpow[j - 1] ), norm_l( partpow[j] ) ), 1 );
2097 66721 : val = L_shl( partpow[j - 1], delta_s );
2098 845595 : FOR( ; i < midband[j]; i++ )
2099 : {
2100 778874 : val = L_shl( Mpy_32_16_1( val, delta ), s1 );
2101 778874 : bandpow[i] = L_shr( val, delta_s );
2102 778874 : move32();
2103 : }
2104 66721 : bandpow[i++] = partpow[j];
2105 66721 : move32();
2106 66721 : partpowLD64M1 = partpowLD64;
2107 66721 : move32();
2108 : }
2109 :
2110 5542 : IF( GT_16( shr( delta, s ), delta_cmp ) )
2111 : {
2112 1385 : delta = 0x4000; // 1.Q14
2113 1385 : move16();
2114 1385 : s1 = 1;
2115 1385 : move16();
2116 : }
2117 :
2118 : /* last half partition */
2119 5542 : val = partpow[stopPartM1]; // Qx
2120 5542 : move32();
2121 101031 : FOR( ; i <= part[stopPartM1]; i++ )
2122 : {
2123 95489 : val = L_shl( Mpy_32_16_1( val, delta ), s1 );
2124 95489 : bandpow[i] = val; // Qx
2125 95489 : move32();
2126 : }
2127 : }
2128 5542 : startBand = add( part[stopPartM1], 1 );
2129 5542 : startPart = stopPart;
2130 5542 : move16();
2131 5542 : stopPart = npart;
2132 5542 : move16();
2133 : }
2134 : }
2135 3082 : }
2136 :
2137 98774 : void scalebands_fx(
2138 : const Word32 *partpow, /* i : Power for each partition Qx*/
2139 : Word16 *part, /* i : Partition upper boundaries (band indices starting from 0) */
2140 : const Word16 npart, /* i : Number of partitions */
2141 : Word16 *midband, /* i : Central band of each partition */
2142 : const Word16 nFFTpart, /* i : Number of FFT partitions */
2143 : const Word16 nband, /* i : Number of bands */
2144 : Word32 *bandpow, /* o : Power for each band Qx*/
2145 : const Word16 flag_fft_en )
2146 : {
2147 : Word16 i, j, s, s1, nint, delta, delta_cmp, delta_s;
2148 : Word16 startBand, startPart, stopPart, stopPartM1;
2149 : Word32 tmp, val, partpowLD64, partpowLD64M1;
2150 :
2151 :
2152 98774 : j = 0;
2153 98774 : move16();
2154 98774 : delta = 0;
2155 98774 : move16();
2156 98774 : partpowLD64M1 = 0L; /* to avoid compilation warnings */
2157 98774 : move32();
2158 :
2159 : /* Interpolate the bin/band-wise levels from the partition levels */
2160 98774 : IF( EQ_16( nband, npart ) )
2161 : {
2162 0 : Copy32( partpow, bandpow, npart );
2163 : }
2164 : ELSE
2165 : {
2166 98774 : startBand = 0;
2167 98774 : move16();
2168 98774 : startPart = 0;
2169 98774 : move16();
2170 98774 : stopPart = nFFTpart;
2171 98774 : move16();
2172 :
2173 199565 : WHILE( LT_16( startBand, nband ) )
2174 : {
2175 100791 : stopPartM1 = sub( stopPart, 1 );
2176 100791 : test();
2177 100791 : IF( ( flag_fft_en != 0 ) || ( GE_16( startPart, nFFTpart ) ) )
2178 : {
2179 : /* first half partition */
2180 98796 : j = startPart;
2181 98796 : move16();
2182 :
2183 200838 : FOR( i = startBand; i <= midband[j]; i++ )
2184 : {
2185 102042 : bandpow[i] = partpow[j]; // Qx
2186 102042 : move32();
2187 : }
2188 98796 : j = add( j, 1 );
2189 :
2190 : /* inner partitions */
2191 98796 : IF( j < stopPart )
2192 : {
2193 98571 : partpowLD64M1 = BASOP_Util_Log2( L_add( partpow[j - 1], DELTA_FX ) );
2194 : }
2195 :
2196 : /* Debug values to check this variable is set. */
2197 98796 : delta = 0x4000; // 1.Q14
2198 98796 : move16();
2199 98796 : delta_cmp = 0x4000; // 1.Q14
2200 98796 : move16();
2201 98796 : s1 = 1;
2202 98796 : move16();
2203 98796 : s = 1;
2204 98796 : move16();
2205 :
2206 5951032 : FOR( ; j < stopPart; j++ )
2207 : {
2208 5852236 : nint = sub( midband[j], midband[j - 1] );
2209 :
2210 : /* log-linear interpolation */
2211 5852236 : IF( NE_32( partpow[j], MAX_32 ) )
2212 : {
2213 5852236 : partpowLD64 = BASOP_Util_Log2( L_add( partpow[j], DELTA_FX ) );
2214 : }
2215 : ELSE
2216 : {
2217 0 : partpowLD64 = BASOP_Util_Log2( partpow[j] );
2218 : }
2219 5852236 : tmp = L_sub( partpowLD64, partpowLD64M1 );
2220 5852236 : tmp = Mpy_32_16_1( tmp, getNormReciprocalWord16( nint ) );
2221 :
2222 : /* scale logarithmic value */
2223 5852236 : tmp = L_sub( tmp, DELTA_SHIFT_LD64 );
2224 5852236 : delta_s = DELTA_SHIFT;
2225 5852236 : move16();
2226 :
2227 6013731 : WHILE( tmp > 0 )
2228 : {
2229 161495 : tmp = L_sub( tmp, 33554432l /*0.015625 Q31*/ ); // Q31
2230 161495 : delta_s = add( delta_s, 1 );
2231 : }
2232 5852236 : delta_cmp = shl( 1, s_max( -15, sub( WORD16_BITS - 1, delta_s ) ) );
2233 :
2234 5852236 : tmp = BASOP_Util_InvLog2( tmp );
2235 5852236 : s = norm_l( tmp );
2236 5852236 : s1 = sub( delta_s, s );
2237 :
2238 5852236 : delta = round_fx_sat( L_shl_sat( tmp, s ) ); // Q(14+s)
2239 : /* Choose scale such that the interpolation start and end point both are representable and add 1 additional bit hr. */
2240 5852236 : delta_s = sub( s_min( norm_l( partpow[j - 1] ), norm_l( partpow[j] ) ), 1 );
2241 5852236 : val = L_shl( partpow[j - 1], delta_s );
2242 27176469 : FOR( ; i < midband[j]; i++ )
2243 : {
2244 21324233 : val = L_shl( Mpy_32_16_1( val, delta ), s1 );
2245 21324233 : bandpow[i] = L_shr( val, delta_s );
2246 21324233 : move32();
2247 : }
2248 5852236 : bandpow[i++] = partpow[j]; // Qx
2249 5852236 : move32();
2250 5852236 : partpowLD64M1 = partpowLD64;
2251 5852236 : move32();
2252 : }
2253 :
2254 98796 : IF( GT_16( shr( delta, s ), delta_cmp ) )
2255 : {
2256 11512 : delta = 0x4000; // 1.Q14
2257 11512 : move16();
2258 11512 : s1 = 1;
2259 11512 : move16();
2260 : }
2261 :
2262 : /* last half partition */
2263 98796 : val = partpow[stopPartM1]; // Qx
2264 98796 : move32();
2265 851852 : FOR( ; i <= part[stopPartM1]; i++ )
2266 : {
2267 753056 : test();
2268 753056 : IF( val != 0 && delta != 0 )
2269 : {
2270 528422 : val = L_shl( Mpy_32_16_1( val, delta ), s1 );
2271 528422 : IF( val == 0 )
2272 : {
2273 2941 : val = 1;
2274 2941 : move32();
2275 : }
2276 : }
2277 : ELSE
2278 : {
2279 224634 : val = 0;
2280 224634 : move32();
2281 : }
2282 753056 : bandpow[i] = val; // Qx
2283 753056 : move32();
2284 : }
2285 : }
2286 100791 : startBand = add( part[stopPartM1], 1 );
2287 100791 : startPart = stopPart;
2288 100791 : move16();
2289 100791 : stopPart = npart;
2290 100791 : move16();
2291 : }
2292 : }
2293 98774 : }
2294 :
2295 : /*-------------------------------------------------------------------
2296 : * getmidbands()
2297 : *
2298 : * Get central band for each partition
2299 : *-------------------------------------------------------------------*/
2300 1125334 : static void getmidbands(
2301 : const Word16 *part, /* i : Partition upper boundaries (band indices starting from 0) */
2302 : const Word16 npart, /* i : Number of partitions */
2303 : Word16 *midband, /* o : Central band of each partition */
2304 : Word16 *psize, /* o : Partition sizes */
2305 : Word16 *psize_norm, /* o : Partition sizes, fractional values */
2306 : Word16 *psize_norm_exp, /* o : Exponent for fractional partition sizes */
2307 : Word16 *psize_inv /* o : Inverse of partition sizes */
2308 : )
2309 : {
2310 : Word16 j, max_psize, shift;
2311 : #ifndef ISSUE_1836_replace_overflow_libcom
2312 : #ifdef BASOP_NOGLOB_DECLARE_LOCAL
2313 : Flag Overflow = 0;
2314 : move32();
2315 : #endif
2316 : #endif
2317 :
2318 :
2319 1125334 : max_psize = psize[0];
2320 1125334 : move16();
2321 1125334 : /* first half partition */ move16();
2322 1125334 : midband[0] = part[0];
2323 1125334 : move16();
2324 1125334 : psize[0] = add( part[0], 1 );
2325 1125334 : move16();
2326 1125334 : psize_inv[0] = getNormReciprocalWord16( psize[0] );
2327 1125334 : move16();
2328 : /* inner partitions */
2329 47823420 : FOR( j = 1; j < npart; j++ )
2330 : {
2331 46698086 : midband[j] = shr( add( add( part[j - 1], 1 ), part[j] ), 1 );
2332 46698086 : move16();
2333 46698086 : psize[j] = sub( part[j], part[j - 1] );
2334 46698086 : move16();
2335 46698086 : psize_inv[j] = getNormReciprocalWord16( psize[j] );
2336 46698086 : move16();
2337 46698086 : if ( GT_16( psize[j], max_psize ) )
2338 : {
2339 14062356 : max_psize = psize[j];
2340 14062356 : move16();
2341 : }
2342 : }
2343 :
2344 1125334 : shift = 9;
2345 1125334 : move16();
2346 1125334 : *psize_norm_exp = sub( 15, shift );
2347 1125334 : move16();
2348 48948754 : FOR( j = 0; j < npart; j++ )
2349 : {
2350 : #ifdef ISSUE_1836_replace_overflow_libcom
2351 47823420 : psize_norm[j] = shl_sat( psize[j], shift ); // Q(15 - psize_norm_exp)
2352 : #else
2353 : psize_norm[j] = shl_o( psize[j], shift, &Overflow ); // Q(15 - psize_norm_exp)
2354 : #endif
2355 47823420 : move16();
2356 : }
2357 : /* minimum_statistics needs fixed exponent of 6 */
2358 1125334 : assert( norm_s( -max_psize ) >= 9 );
2359 1125334 : }
2360 :
2361 : /*-------------------------------------------------------------------
2362 : * AnalysisSTFT()
2363 : *
2364 : * STFT analysis filterbank
2365 : *-------------------------------------------------------------------*/
2366 :
2367 15 : void AnalysisSTFT(
2368 : const Word16 *timeDomainInput, /* i : pointer to time signal Q(Q)*/
2369 : Word16 Q,
2370 : Word32 *fftBuffer, /* o : FFT bins */
2371 : Word16 *fftBuffer_exp, /* i : exponent of FFT bins */
2372 : HANDLE_FD_CNG_COM hFdCngCom /* i/o: FD_CNG structure containing all buffers and variables */
2373 : )
2374 : {
2375 : Word16 i, len;
2376 : Word16 len2;
2377 : const PWord16 *olapWin;
2378 : Word16 *olapBuffer;
2379 :
2380 :
2381 15 : assert( ( hFdCngCom->fftlen >> 1 ) == hFdCngCom->frameSize );
2382 :
2383 : /* pointer inititialization */
2384 15 : assert( hFdCngCom->olapBufferAna != NULL );
2385 15 : olapBuffer = hFdCngCom->olapBufferAna;
2386 15 : olapWin = hFdCngCom->olapWinAna;
2387 :
2388 : /* olapWin factor is scaled with one bit */
2389 15 : *fftBuffer_exp = 1;
2390 15 : move16();
2391 15 : len = sub( hFdCngCom->fftlen, hFdCngCom->frameSize );
2392 15 : assert( len <= 320 ); /* see size of olapBuffer */
2393 :
2394 : /* Window the signal */
2395 15 : len2 = shr( len, 1 );
2396 2191 : FOR( i = 0; i < len2; i++ )
2397 : {
2398 2176 : move32();
2399 2176 : move32();
2400 2176 : fftBuffer[i] = L_mult( olapBuffer[i], mult_r( olapWin[i].v.im, 23170 /*1.4142135623730950488016887242097 Q14*/ ) );
2401 2176 : fftBuffer[i + len2] = L_mult( olapBuffer[i + len2], mult_r( olapWin[len2 - 1 - i].v.re, 23170 /*1.4142135623730950488016887242097 Q14*/ ) );
2402 : }
2403 15 : len2 = shr( hFdCngCom->frameSize, 1 );
2404 2191 : FOR( i = 0; i < len2; i++ )
2405 : {
2406 2176 : move32();
2407 2176 : move32();
2408 2176 : fftBuffer[i + len] = L_mult( shr_sat( timeDomainInput[i], Q ), mult_r( olapWin[i].v.re, 23170 /*1.4142135623730950488016887242097 Q14*/ ) );
2409 2176 : fftBuffer[i + len + len2] = L_mult( shr_sat( timeDomainInput[i + len2], Q ), mult_r( olapWin[len2 - 1 - i].v.im, 23170 /*1.4142135623730950488016887242097 Q14*/ ) );
2410 : }
2411 :
2412 : /* Perform FFT */
2413 15 : BASOP_rfft( fftBuffer, hFdCngCom->fftlen, fftBuffer_exp, -1 );
2414 :
2415 4367 : FOR( i = 0; i < len; i++ )
2416 : {
2417 4352 : olapBuffer[i] = shr_sat( timeDomainInput[sub( hFdCngCom->frameSize, len ) + i], Q );
2418 4352 : move16();
2419 : }
2420 15 : }
2421 :
2422 : /*-------------------------------------------------------------------
2423 : * AnalysisSTFT_fx()
2424 : *
2425 : * STFT analysis filterbank
2426 : *-------------------------------------------------------------------*/
2427 90653 : void AnalysisSTFT_fx(
2428 : const Word16 *timeDomainInput, // Q(Q)
2429 : Word16 Q,
2430 : Word32 *fftBuffer, /* o : FFT bins */
2431 : Word16 *fftBuffer_exp, /* i : exponent of FFT bins */
2432 : HANDLE_FD_CNG_COM hFdCngCom /* i/o: FD_CNG structure containing all buffers and variables */
2433 : )
2434 : {
2435 : Word16 i, len;
2436 : const Word32 *olapWin;
2437 : Word16 *olapBuffer;
2438 :
2439 90653 : assert( ( hFdCngCom->fftlen >> 1 ) == hFdCngCom->frameSize );
2440 :
2441 : /* pointer inititialization */
2442 90653 : assert( hFdCngCom->olapBufferAna != NULL );
2443 90653 : olapBuffer = hFdCngCom->olapBufferAna_fx;
2444 90653 : olapWin = hFdCngCom->olapWinAna_fx;
2445 :
2446 : /* olapWin factor is scaled with one bit */
2447 90653 : *fftBuffer_exp = 1;
2448 90653 : move16();
2449 90653 : len = sub( hFdCngCom->fftlen, hFdCngCom->frameSize );
2450 90653 : assert( len <= 320 ); /* see size of olapBuffer */
2451 :
2452 : /* Window the signal */
2453 90653 : Copy( olapBuffer + hFdCngCom->frameSize, olapBuffer, len );
2454 :
2455 26324445 : FOR( i = 0; i < hFdCngCom->frameSize; i++ )
2456 : {
2457 26233792 : olapBuffer[i + len] = shr_sat( timeDomainInput[i], Q ); // Values above MAX_16 seen with 10dB tests in float code
2458 26233792 : move16();
2459 : }
2460 :
2461 : /* Window the signal */
2462 90653 : v_L_mult_3216( olapWin, olapBuffer, fftBuffer, hFdCngCom->fftlen );
2463 :
2464 52558237 : FOR( i = 0; i < hFdCngCom->fftlen; i++ )
2465 : {
2466 52467584 : fftBuffer[i] = L_shr( fftBuffer[i], 11 );
2467 52467584 : move32();
2468 : }
2469 90653 : *fftBuffer_exp = WORD16_BITS + 11;
2470 90653 : move16();
2471 :
2472 : /* Perform FFT */
2473 90653 : RFFTN_fx( fftBuffer, hFdCngCom->fftSineTab_fx, hFdCngCom->fftlen, -1 );
2474 :
2475 :
2476 90653 : return;
2477 : }
2478 :
2479 : /*-------------------------------------------------------------------
2480 : * SynthesisSTFT_enc_ivas_fx()
2481 : *
2482 : * STFT synthesis filterbank
2483 : *-------------------------------------------------------------------*/
2484 :
2485 0 : void SynthesisSTFT_enc_ivas_fx(
2486 : Word32 *fftBuffer, /* i : pointer to FFT bins */
2487 : Word16 fftBufferExp, /* i : exponent of FFT bins */
2488 : Word16 *timeDomainOutput, /* o : pointer to time domain signal */
2489 : Word16 *olapBuffer, /* i/o : pointer to overlap buffer */
2490 : const PWord16 *olapWin, /* i : pointer to overlap window */
2491 : Word16 tcx_transition,
2492 : HANDLE_FD_CNG_COM hFdCngCom, /* i/o : pointer to FD_CNG structure containing all buffers and variables */
2493 : Word16 gen_exc,
2494 : Word16 *Q_new, /* i : Q of generated exc_cng */
2495 : const Word16 element_mode, /* i : element mode */
2496 : const Word16 nchan_out /* i : number of output channels */
2497 : )
2498 : {
2499 : Word16 i, len, scale, tmp;
2500 : Word16 len2, len3, len4;
2501 : Word16 buf[M + 1 + L_FRAME16k];
2502 :
2503 :
2504 : /* Perform IFFT */
2505 0 : scale = 0;
2506 0 : move16();
2507 0 : BASOP_rfft( fftBuffer, hFdCngCom->fftlen, &scale, 1 );
2508 0 : fftBufferExp = add( fftBufferExp, scale );
2509 0 : hFdCngCom->fftBuffer_exp = fftBufferExp;
2510 0 : move16();
2511 :
2512 0 : fftBufferExp = add( fftBufferExp, hFdCngCom->fftlenShift );
2513 :
2514 : /* Perform overlap-add */
2515 : /* Handle overlap in P/S domain for stereo */
2516 0 : test();
2517 0 : test();
2518 0 : IF( ( EQ_16( element_mode, IVAS_CPE_TD ) || EQ_16( element_mode, IVAS_CPE_DFT ) ) && EQ_16( nchan_out, 2 ) )
2519 : {
2520 0 : Copy( olapBuffer + 3 * hFdCngCom->frameSize / 4 - ( M + 1 ), buf, hFdCngCom->frameSize + M + 1 );
2521 0 : set16_fx( olapBuffer, 0, hFdCngCom->fftlen );
2522 : }
2523 : ELSE
2524 : {
2525 0 : Copy( olapBuffer + hFdCngCom->frameSize, olapBuffer, hFdCngCom->frameSize );
2526 0 : set16_fx( olapBuffer + hFdCngCom->frameSize, 0, hFdCngCom->frameSize );
2527 : }
2528 0 : len2 = shr( hFdCngCom->fftlen, 2 );
2529 0 : len4 = shr( hFdCngCom->fftlen, 3 );
2530 0 : len3 = add( len2, len4 );
2531 0 : len = add( hFdCngCom->frameSize, len4 );
2532 0 : IF( tcx_transition )
2533 : {
2534 0 : FOR( i = 0; i < len; i++ )
2535 : {
2536 0 : olapBuffer[i] = round_fx_sat( L_shl_sat( fftBuffer[i], sub( fftBufferExp, 15 ) ) ); // Q(15 - fftBufferExp)
2537 0 : move16();
2538 : }
2539 : }
2540 : ELSE
2541 : {
2542 0 : FOR( i = 0; i < len4; i++ )
2543 : {
2544 0 : olapBuffer[i + 1 * len4] = add_sat( olapBuffer[i + 1 * len4], mult_r( round_fx_sat( L_shl_sat( fftBuffer[i + 1 * len4], sub( fftBufferExp, 15 ) ) ), olapWin[i].v.im ) );
2545 0 : move16();
2546 0 : olapBuffer[i + 2 * len4] = add_sat( olapBuffer[i + 2 * len4], mult_r( round_fx_sat( L_shl_sat( fftBuffer[i + 2 * len4], sub( fftBufferExp, 15 ) ) ), olapWin[len4 - 1 - i].v.re ) );
2547 0 : move16();
2548 : }
2549 0 : FOR( i = len3; i < len; i++ )
2550 : {
2551 0 : olapBuffer[i] = round_fx_sat( L_shl_sat( fftBuffer[i], sub( fftBufferExp, 15 ) ) ); // Q(15 - fftBufferExp)
2552 0 : move16();
2553 : }
2554 : }
2555 :
2556 0 : FOR( i = 0; i < len4; i++ )
2557 : {
2558 0 : olapBuffer[i + 5 * len4] = mult_r( round_fx_sat( L_shl_sat( fftBuffer[i + 5 * len4], sub( fftBufferExp, 15 ) ) ), olapWin[i].v.re );
2559 0 : move16();
2560 0 : olapBuffer[i + 6 * len4] = mult_r( round_fx_sat( L_shl_sat( fftBuffer[i + 6 * len4], sub( fftBufferExp, 15 ) ) ), olapWin[len4 - 1 - i].v.im );
2561 0 : move16();
2562 : }
2563 :
2564 0 : len = add( len, len2 );
2565 0 : FOR( i = len; i < hFdCngCom->fftlen; i++ )
2566 : {
2567 0 : olapBuffer[i] = 0;
2568 0 : move16();
2569 : }
2570 :
2571 : /* Get time-domain signal */
2572 0 : FOR( i = 0; i < hFdCngCom->frameSize; i++ )
2573 : {
2574 0 : timeDomainOutput[i] = mult_r( olapBuffer[i + len4], hFdCngCom->fftlenFac );
2575 0 : move16();
2576 : }
2577 : /* Generate excitation */
2578 0 : test();
2579 0 : test();
2580 0 : IF( ( EQ_16( element_mode, IVAS_CPE_TD ) || EQ_16( element_mode, IVAS_CPE_DFT ) ) && EQ_16( nchan_out, 2 ) )
2581 : {
2582 0 : FOR( i = 0; i < hFdCngCom->frameSize / 2; i++ )
2583 : {
2584 0 : buf[i + ( M + 1 )] = add( buf[i + ( M + 1 )], olapBuffer[i + hFdCngCom->frameSize / 4] );
2585 0 : move16();
2586 : }
2587 :
2588 0 : FOR( i = 0; i < M + 1 + hFdCngCom->frameSize; i++ )
2589 : {
2590 0 : buf[i] = mult_r( buf[i], hFdCngCom->fftlenFac );
2591 0 : move16();
2592 : }
2593 : }
2594 : ELSE
2595 : {
2596 0 : FOR( i = 0; i < M + 1 + hFdCngCom->frameSize; i++ )
2597 : {
2598 0 : buf[i] = mult_r( olapBuffer[i + len4 - M - 1], hFdCngCom->fftlenFac );
2599 0 : move16();
2600 : }
2601 0 : tmp = buf[0];
2602 0 : move16();
2603 : }
2604 0 : IF( EQ_16( gen_exc, 1 ) )
2605 : {
2606 :
2607 0 : E_UTIL_f_preemph2( sub( *Q_new, 1 ), buf + 1, PREEMPH_FAC, M + hFdCngCom->frameSize, &tmp );
2608 0 : Residu3_fx( hFdCngCom->A_cng, buf + 1 + M, hFdCngCom->exc_cng, hFdCngCom->frameSize, 1 );
2609 : }
2610 0 : IF( EQ_16( gen_exc, 2 ) )
2611 : {
2612 0 : *Q_new = E_UTIL_f_preemph3( buf + 1, PREEMPH_FAC, M + hFdCngCom->frameSize, &tmp, 1 );
2613 0 : move16();
2614 0 : Residu3_fx( hFdCngCom->A_cng, buf + 1 + M, hFdCngCom->exc_cng, hFdCngCom->frameSize, 1 );
2615 : }
2616 0 : }
2617 :
2618 : /*-------------------------------------------------------------------
2619 : * SynthesisSTFT()
2620 : *
2621 : * STFT synthesis filterbank
2622 : *-------------------------------------------------------------------*/
2623 :
2624 1021 : void SynthesisSTFT(
2625 : Word32 *fftBuffer, /* i : pointer to FFT bins */
2626 : Word16 fftBufferExp, /* i : exponent of FFT bins */
2627 : Word16 *timeDomainOutput, /* o : pointer to time domain signal Qx*/
2628 : Word16 *olapBuffer, /* i/o : pointer to overlap buffer */
2629 : const PWord16 *olapWin, /* i : pointer to overlap window */
2630 : Word16 tcx_transition,
2631 : HANDLE_FD_CNG_COM hFdCngCom, /* i/o : pointer to FD_CNG structure containing all buffers and variables */
2632 : Word16 gen_exc,
2633 : Word16 *Q_new,
2634 : const Word16 element_mode, /* i : element mode */
2635 : const Word16 nchan_out /* i : number of output channels */
2636 : )
2637 : {
2638 : Word16 i, len, scale, tmp;
2639 : Word16 len2, len3, len4;
2640 : Word16 buf[M + 1 + L_FRAME16k];
2641 :
2642 :
2643 : /* Perform IFFT */
2644 1021 : scale = 0;
2645 1021 : move16();
2646 1021 : BASOP_rfft( fftBuffer, hFdCngCom->fftlen, &scale, 1 );
2647 1021 : fftBufferExp = add( fftBufferExp, scale );
2648 1021 : hFdCngCom->fftBuffer_exp = fftBufferExp;
2649 1021 : move16();
2650 :
2651 1021 : fftBufferExp = add( fftBufferExp, hFdCngCom->fftlenShift );
2652 :
2653 : /* Perform overlap-add */
2654 : /* Handle overlap in P/S domain for stereo */
2655 1021 : IF( ( EQ_16( element_mode, IVAS_CPE_TD ) || EQ_16( element_mode, IVAS_CPE_DFT ) ) && EQ_16( nchan_out, 2 ) )
2656 : {
2657 0 : Copy( olapBuffer + 3 * hFdCngCom->frameSize / 4 - ( M + 1 ), buf, hFdCngCom->frameSize + M + 1 );
2658 0 : set16_fx( olapBuffer, 0, hFdCngCom->fftlen );
2659 : }
2660 : ELSE
2661 : {
2662 1021 : Copy( olapBuffer + hFdCngCom->frameSize, olapBuffer, hFdCngCom->frameSize );
2663 1021 : set16_fx( olapBuffer + hFdCngCom->frameSize, 0, hFdCngCom->frameSize );
2664 : }
2665 1021 : len2 = shr( hFdCngCom->fftlen, 2 );
2666 1021 : len4 = shr( hFdCngCom->fftlen, 3 );
2667 1021 : len3 = add( len2, len4 );
2668 1021 : len = add( hFdCngCom->frameSize, len4 );
2669 1021 : IF( tcx_transition )
2670 : {
2671 0 : FOR( i = 0; i < len; i++ )
2672 : {
2673 0 : olapBuffer[i] = round_fx_sat( L_shl_sat( fftBuffer[i], fftBufferExp - 15 ) ); // Q(15 - fftBufferExp)
2674 0 : move16();
2675 : }
2676 : }
2677 : ELSE
2678 : {
2679 66365 : FOR( i = 0; i < len4; i++ )
2680 : {
2681 65344 : olapBuffer[i + 1 * len4] = add_sat( olapBuffer[i + 1 * len4], mult_r( round_fx_sat( L_shl_sat( fftBuffer[i + 1 * len4], fftBufferExp - 15 ) ), olapWin[i].v.im ) );
2682 65344 : move16();
2683 65344 : olapBuffer[i + 2 * len4] = add_sat( olapBuffer[i + 2 * len4], mult_r( round_fx_sat( L_shl_sat( fftBuffer[i + 2 * len4], fftBufferExp - 15 ) ), olapWin[len4 - 1 - i].v.re ) );
2684 65344 : move16();
2685 : }
2686 131709 : FOR( i = len3; i < len; i++ )
2687 : {
2688 130688 : olapBuffer[i] = round_fx_sat( L_shl_sat( fftBuffer[i], fftBufferExp - 15 ) );
2689 : }
2690 : }
2691 :
2692 66365 : FOR( i = 0; i < len4; i++ )
2693 : {
2694 65344 : olapBuffer[i + 5 * len4] = mult_r( round_fx_sat( L_shl_sat( fftBuffer[i + 5 * len4], fftBufferExp - 15 ) ), olapWin[i].v.re );
2695 65344 : move16();
2696 65344 : olapBuffer[i + 6 * len4] = mult_r( round_fx_sat( L_shl_sat( fftBuffer[i + 6 * len4], fftBufferExp - 15 ) ), olapWin[len4 - 1 - i].v.im );
2697 65344 : move16();
2698 : }
2699 :
2700 1021 : len = add( len, len2 );
2701 66365 : FOR( i = len; i < hFdCngCom->fftlen; i++ )
2702 : {
2703 65344 : olapBuffer[i] = 0;
2704 65344 : move16();
2705 : }
2706 :
2707 : /* Get time-domain signal */
2708 262397 : FOR( i = 0; i < hFdCngCom->frameSize; i++ )
2709 : {
2710 261376 : timeDomainOutput[i] = mult_r( olapBuffer[i + len4], hFdCngCom->fftlenFac );
2711 261376 : move16();
2712 : }
2713 : /* Generate excitation */
2714 : {
2715 279754 : FOR( i = 0; i < M + 1 + hFdCngCom->frameSize; i++ )
2716 : {
2717 278733 : buf[i] = mult_r( olapBuffer[i + len4 - M - 1], hFdCngCom->fftlenFac );
2718 278733 : move16();
2719 : }
2720 1021 : tmp = buf[0];
2721 : }
2722 1021 : IF( EQ_16( gen_exc, 1 ) )
2723 : {
2724 :
2725 0 : E_UTIL_f_preemph2( *Q_new - 1, buf + 1, PREEMPH_FAC, M + hFdCngCom->frameSize, &tmp );
2726 0 : Residu3_fx( hFdCngCom->A_cng, buf + 1 + M, hFdCngCom->exc_cng, hFdCngCom->frameSize, 1 );
2727 : }
2728 1021 : IF( EQ_16( gen_exc, 2 ) )
2729 : {
2730 0 : *Q_new = E_UTIL_f_preemph3( buf + 1, PREEMPH_FAC, M + hFdCngCom->frameSize, &tmp, 1 );
2731 0 : Residu3_fx( hFdCngCom->A_cng, buf + 1 + M, hFdCngCom->exc_cng, hFdCngCom->frameSize, 1 );
2732 : }
2733 1021 : }
2734 :
2735 18381 : void SynthesisSTFT_ivas_fx(
2736 : Word32 *fftBuffer, /* i : pointer to FFT bins */
2737 : Word16 fftBufferExp, /* i : exponent of FFT bins */
2738 : Word16 *timeDomainOutput, /* o : pointer to time domain signal Qx*/
2739 : Word16 *olapBuffer, /* i/o : pointer to overlap buffer */
2740 : const PWord16 *olapWin, /* i : pointer to overlap window */
2741 : Word16 tcx_transition,
2742 : HANDLE_FD_CNG_COM hFdCngCom, /* i/o : pointer to FD_CNG structure containing all buffers and variables */
2743 : Word16 gen_exc,
2744 : Word16 *Q_new,
2745 : const Word16 element_mode, /* i : element mode */
2746 : const Word16 nchan_out /* i : number of output channels */
2747 : )
2748 : {
2749 : Word16 i, len, scale, tmp, sft, exp, tmp1;
2750 : Word16 len2, len3, len4;
2751 : Word16 buf[M + 1 + L_FRAME16k];
2752 : Word32 L_olapBuffer[FFTLEN];
2753 :
2754 :
2755 : /* Perform IFFT */
2756 18381 : scale = 0;
2757 18381 : move16();
2758 18381 : BASOP_rfft( fftBuffer, hFdCngCom->fftlen, &scale, 1 );
2759 18381 : fftBufferExp = add( fftBufferExp, scale );
2760 18381 : hFdCngCom->fftBuffer_exp = fftBufferExp;
2761 18381 : move16();
2762 :
2763 18381 : sft = getScaleFactor32( fftBuffer, hFdCngCom->fftlen );
2764 18381 : scale_sig32( fftBuffer, hFdCngCom->fftlen, sft );
2765 18381 : fftBufferExp = sub( fftBufferExp, sft );
2766 18381 : hFdCngCom->fftBuffer_exp = fftBufferExp;
2767 18381 : move16();
2768 :
2769 18381 : fftBufferExp = add( fftBufferExp, hFdCngCom->fftlenShift );
2770 :
2771 : /* Perform overlap-add */
2772 : /* Handle overlap in P/S domain for stereo */
2773 18381 : test();
2774 18381 : test();
2775 18381 : IF( ( EQ_16( element_mode, IVAS_CPE_TD ) || EQ_16( element_mode, IVAS_CPE_DFT ) ) && EQ_16( nchan_out, 2 ) )
2776 : {
2777 0 : Copy( olapBuffer + 3 * hFdCngCom->frameSize / 4 - ( M + 1 ), buf, hFdCngCom->frameSize + M + 1 );
2778 0 : set16_fx( olapBuffer, 0, hFdCngCom->fftlen );
2779 : }
2780 : ELSE
2781 : {
2782 18381 : Copy_Scale_sig_16_32_no_sat( olapBuffer + hFdCngCom->frameSize, L_olapBuffer, hFdCngCom->frameSize, add( sub( hFdCngCom->olapBufferSynth_exp, sub( 15, hFdCngCom->fftlenShift ) ), 16 ) ); /* exp: hFdCngCom->olapBufferSynth_exp -> (15 - fftlenShift) */
2783 18381 : set16_fx( olapBuffer + hFdCngCom->frameSize, 0, hFdCngCom->frameSize );
2784 : }
2785 18381 : len2 = shr( hFdCngCom->fftlen, 2 );
2786 18381 : len4 = shr( hFdCngCom->fftlen, 3 );
2787 18381 : len3 = add( len2, len4 );
2788 18381 : len = add( hFdCngCom->frameSize, len4 );
2789 18381 : tmp1 = sub( fftBufferExp, 15 );
2790 18381 : IF( tcx_transition )
2791 : {
2792 607052 : FOR( i = 0; i < len; i++ )
2793 : {
2794 605280 : L_olapBuffer[i] = L_shl_sat( fftBuffer[i], tmp1 ); /* Q(31 - (15 - fftlenShift) ) */
2795 605280 : move32();
2796 : }
2797 : }
2798 : ELSE
2799 : {
2800 1256417 : FOR( i = 0; i < len4; i++ )
2801 : {
2802 1239808 : L_olapBuffer[i + 1 * len4] = Madd_32_16( L_olapBuffer[i + 1 * len4], L_shl_sat( fftBuffer[i + 1 * len4], tmp1 ), olapWin[i].v.im ); /* Q(31 - (15 - fftlenShift) ) */
2803 1239808 : move32();
2804 1239808 : L_olapBuffer[i + 2 * len4] = Madd_32_16( L_olapBuffer[i + 2 * len4], L_shl_sat( fftBuffer[i + 2 * len4], tmp1 ), olapWin[len4 - 1 - i].v.re ); /* Q(31 - (15 - fftlenShift) ) */
2805 1239808 : move32();
2806 : }
2807 2496225 : FOR( i = len3; i < len; i++ )
2808 : {
2809 2479616 : L_olapBuffer[i] = L_shl_sat( fftBuffer[i], tmp1 ); /* Q(31 - (15 - fftlenShift) ) */
2810 2479616 : move32();
2811 : }
2812 : }
2813 :
2814 1379245 : FOR( i = 0; i < len4; i++ )
2815 : {
2816 1360864 : L_olapBuffer[i + 5 * len4] = Mpy_32_16_1( L_shl_sat( fftBuffer[i + 5 * len4], tmp1 ), olapWin[i].v.re ); /* Q(31 - (15 - fftlenShift) ) */
2817 1360864 : move32();
2818 1360864 : L_olapBuffer[i + 6 * len4] = Mpy_32_16_1( L_shl_sat( fftBuffer[i + 6 * len4], tmp1 ), olapWin[len4 - 1 - i].v.im ); /* Q(31 - (15 - fftlenShift) ) */
2819 1360864 : move32();
2820 : }
2821 :
2822 18381 : len = add( len, len2 );
2823 1379245 : FOR( i = len; i < hFdCngCom->fftlen; i++ )
2824 : {
2825 1360864 : L_olapBuffer[i] = 0;
2826 1360864 : move32();
2827 : }
2828 :
2829 18381 : sft = L_norm_arr( L_olapBuffer, hFdCngCom->fftlen );
2830 18381 : IF( NE_16( sft, 31 ) )
2831 : {
2832 17096 : Copy_Scale_sig32_16( L_olapBuffer, olapBuffer, hFdCngCom->fftlen, sft );
2833 17096 : hFdCngCom->olapBufferSynth_exp = sub( sub( 15, hFdCngCom->fftlenShift ), sft );
2834 17096 : move16();
2835 : }
2836 :
2837 : /* Get time-domain signal */
2838 5461837 : FOR( i = 0; i < hFdCngCom->frameSize; i++ )
2839 : {
2840 5443456 : timeDomainOutput[i] = mult_r( olapBuffer[i + len4], hFdCngCom->fftlenFac );
2841 5443456 : move16();
2842 : }
2843 : /* Generate excitation */
2844 : {
2845 5774314 : FOR( i = 0; i < M + 1 + hFdCngCom->frameSize; i++ )
2846 : {
2847 5755933 : buf[i] = mult_r( olapBuffer[i + len4 - M - 1], hFdCngCom->fftlenFac );
2848 5755933 : move16();
2849 : }
2850 18381 : exp = add( hFdCngCom->olapBufferSynth_exp, hFdCngCom->fftlenShift );
2851 18381 : tmp = buf[0];
2852 : }
2853 18381 : IF( EQ_16( gen_exc, 1 ) )
2854 : {
2855 18381 : E_UTIL_f_preemph2( 0, buf + 1, PREEMPH_FAC, add( M, hFdCngCom->frameSize ), &tmp );
2856 18381 : Residu3_fx( hFdCngCom->A_cng, buf + 1 + M, hFdCngCom->exc_cng, hFdCngCom->frameSize, 0 );
2857 18381 : *Q_new = sub( 15, exp );
2858 18381 : move16();
2859 : }
2860 18381 : IF( EQ_16( gen_exc, 2 ) )
2861 : {
2862 0 : *Q_new = E_UTIL_f_preemph3( buf + 1, PREEMPH_FAC, add( M, hFdCngCom->frameSize ), &tmp, 0 );
2863 0 : move16();
2864 0 : Residu3_fx( hFdCngCom->A_cng, buf + 1 + M, hFdCngCom->exc_cng, hFdCngCom->frameSize, negate( *Q_new ) );
2865 0 : *Q_new = sub( 15, exp );
2866 0 : move16();
2867 : }
2868 18381 : }
2869 :
2870 : /**************************************************************************************
2871 : * Compute some values used in the bias correction of the minimum statistics algorithm *
2872 : **************************************************************************************/
2873 14270 : void mhvals(
2874 : const Word16 d,
2875 : Word16 *m /*, Word16 * h*/
2876 : )
2877 : {
2878 : Word16 i, j;
2879 14270 : Word16 len = SIZE_SCALE_TABLE_CN;
2880 14270 : move16();
2881 :
2882 :
2883 14270 : assert( d == 72 || d == 12 ); /* function only tested for d==72 and d==12) */
2884 14270 : i = 0;
2885 14270 : move16();
2886 121295 : FOR( i = 0; i < len; i++ )
2887 : {
2888 121295 : IF( LE_16( d, d_array[i] ) )
2889 : {
2890 14270 : BREAK;
2891 : }
2892 : }
2893 14270 : IF( EQ_16( i, len ) )
2894 : {
2895 0 : i = sub( len, 1 );
2896 0 : j = i;
2897 0 : move16();
2898 : }
2899 : ELSE
2900 : {
2901 14270 : j = sub( i, 1 );
2902 : }
2903 14270 : IF( EQ_16( d, d_array[i] ) )
2904 : {
2905 0 : *m = m_array[i];
2906 0 : move16();
2907 : }
2908 : ELSE
2909 : {
2910 : Word32 qi_m, qj_m, q_m, tmp1_m, tmp2_m;
2911 : Word16 qi_e, qj_e, q_e, tmp1_e, tmp2_e, tmp1_w16_m, tmp1_w16_e, shift;
2912 :
2913 :
2914 : /* d_array has exponent 15 */
2915 14270 : qj_e = 15;
2916 14270 : move16();
2917 14270 : qj_m = L_deposit_h( d_array[i - 1] );
2918 :
2919 14270 : qi_e = 15;
2920 14270 : move16();
2921 14270 : qi_m = L_deposit_h( d_array[i] );
2922 :
2923 14270 : q_e = 15;
2924 14270 : move16();
2925 14270 : q_m = L_deposit_h( d );
2926 :
2927 14270 : qj_m = Sqrt32( qj_m, &qj_e );
2928 14270 : qi_m = Sqrt32( qi_m, &qi_e );
2929 14270 : q_m = Sqrt32( q_m, &q_e );
2930 :
2931 14270 : tmp1_m = Mpy_32_32( qi_m, qj_m );
2932 14270 : tmp1_e = add( qi_e, qj_e );
2933 14270 : tmp1_m = L_deposit_h( BASOP_Util_Divide3232_Scale( tmp1_m, q_m, &shift ) );
2934 14270 : tmp1_e = sub( tmp1_e, q_e );
2935 14270 : tmp1_e = add( tmp1_e, shift );
2936 14270 : tmp1_m = BASOP_Util_Add_Mant32Exp( tmp1_m, tmp1_e, L_negate( qj_m ), qj_e, &tmp1_e );
2937 :
2938 14270 : tmp2_m = BASOP_Util_Add_Mant32Exp( qi_m, qi_e, L_negate( qj_m ), qj_e, &tmp2_e );
2939 14270 : tmp1_w16_m = round_fx( tmp2_m );
2940 14270 : tmp1_w16_e = tmp2_e;
2941 14270 : move16();
2942 14270 : BASOP_Util_Divide_MantExp( sub( m_array[j], m_array[i] ), 0, tmp1_w16_m, tmp1_w16_e, &tmp1_w16_m, &tmp1_w16_e );
2943 :
2944 14270 : tmp2_m = Mpy_32_16_1( tmp1_m, tmp1_w16_m );
2945 14270 : tmp2_e = add( tmp1_e, tmp1_w16_e );
2946 :
2947 14270 : tmp2_m = BASOP_Util_Add_Mant32Exp( tmp2_m, tmp2_e, L_deposit_h( m_array[i] ), 0, &tmp2_e );
2948 14270 : assert( tmp2_e == 0 );
2949 14270 : *m = extract_h( tmp2_m );
2950 14270 : move32();
2951 : }
2952 14270 : }
2953 :
2954 : /*-------------------------------------------------------------------
2955 : * rand_gauss()
2956 : *
2957 : * Random generator with Gaussian distribution with mean 0 and std 1
2958 : * Returns:
2959 : * random signal format Q3.29
2960 : *-------------------------------------------------------------------*/
2961 40372388 : Word32 rand_gauss( Word16 *seed )
2962 : {
2963 : Word32 temp;
2964 : Word16 loc_seed;
2965 :
2966 : /* This unrolled version reduces the cycles from 17 to 10 */
2967 40372388 : loc_seed = extract_l( L_mac0( 13849, *seed, 31821 ) );
2968 40372388 : temp = L_deposit_l( loc_seed );
2969 :
2970 40372388 : loc_seed = extract_l( L_mac0( 13849, loc_seed, 31821 ) );
2971 40372388 : temp = L_msu0( temp, loc_seed, -1 );
2972 :
2973 40372388 : loc_seed = extract_l( L_mac0( 13849, loc_seed, 31821 ) );
2974 40372388 : temp = L_msu0( temp, loc_seed, -1 );
2975 :
2976 40372388 : *seed = loc_seed;
2977 40372388 : move16();
2978 40372388 : return L_shl( temp, WORD16_BITS - CNG_RAND_GAUSS_SHIFT );
2979 : }
2980 :
2981 : /*-------------------------------------------------------------------
2982 : * lpc_from_spectrum()
2983 : *
2984 : *
2985 : *-------------------------------------------------------------------*/
2986 :
2987 3802 : void lpc_from_spectrum(
2988 : HANDLE_FD_CNG_COM hFdCngCom,
2989 : const Word16 start, /*i : start band*/
2990 : const Word16 stop, /*i : stop band*/
2991 : const Word16 preemph_fac /*i : preemphase factor format Q1.15*/
2992 : )
2993 : {
2994 : Word16 i, s1, s2, s3, fftlen2, scale, fftlen4, fftlen8, len, step, preemph_fac2;
2995 : Word32 maxVal, r[32], fftBuffer[FFTLEN], *ptr, *pti, nf;
2996 : Word16 tmp, r_h[32], r_l[32];
2997 : const PWord16 *table;
2998 :
2999 3802 : Word32 *powspec = hFdCngCom->cngNoiseLevel; /*i : pointer to noise levels format Q5.27*/
3000 3802 : Word16 powspec_exp = hFdCngCom->cngNoiseLevelExp;
3001 3802 : move16();
3002 3802 : Word16 fftlen = hFdCngCom->fftlen; /*i : size of fft*/
3003 3802 : Word16 *A = hFdCngCom->A_cng; /*o : lpc coefficients format Q3.12*/
3004 3802 : move16();
3005 3802 : Word16 lpcorder = M;
3006 3802 : move16();
3007 :
3008 3802 : scale = 0;
3009 3802 : move16();
3010 3802 : fftlen2 = shr( fftlen, 1 );
3011 3802 : fftlen4 = shr( fftlen, 2 );
3012 3802 : fftlen8 = shr( fftlen, 3 );
3013 :
3014 : /* Power Spectrum */
3015 3802 : maxVal = 0;
3016 3802 : move32();
3017 3802 : len = sub( stop, start );
3018 1121574 : FOR( i = 0; i < len; i++ )
3019 : {
3020 1117772 : maxVal = L_max( maxVal, L_abs( powspec[i] ) );
3021 : }
3022 3802 : s1 = norm_l( maxVal );
3023 3802 : nf = L_shr_r_sat( 1099511680l /*1e-3f Q40*/, add( sub( powspec_exp, s1 ), 9 ) );
3024 3802 : ptr = fftBuffer;
3025 3802 : pti = fftBuffer + 1;
3026 :
3027 11406 : FOR( i = 0; i < start; i++ )
3028 : {
3029 7604 : *ptr = nf;
3030 7604 : move32();
3031 7604 : *pti = L_deposit_l( 0 );
3032 7604 : move32();
3033 7604 : ptr += 2;
3034 7604 : pti += 2;
3035 : }
3036 :
3037 1121574 : FOR( ; i < stop; i++ )
3038 : {
3039 1117772 : *ptr = L_max( nf, L_shl( powspec[i - start], s1 ) );
3040 1117772 : move32();
3041 1117772 : *pti = L_deposit_l( 0 );
3042 1117772 : move32();
3043 1117772 : ptr += 2;
3044 1117772 : pti += 2;
3045 : }
3046 :
3047 3994 : FOR( ; i < fftlen2; i++ )
3048 : {
3049 192 : *ptr = nf;
3050 192 : move32();
3051 192 : *pti = L_deposit_l( 0 );
3052 192 : move32();
3053 192 : ptr += 2;
3054 192 : pti += 2;
3055 : }
3056 :
3057 3802 : fftBuffer[1] = nf;
3058 3802 : move32();
3059 :
3060 : /* Pre-emphasis */
3061 :
3062 3802 : BASOP_getTables( &table, NULL, &step, fftlen4 );
3063 3802 : tmp = round_fx( L_shr( L_add( 0x40000000, L_mult0( preemph_fac, preemph_fac ) ), 1 ) );
3064 3802 : preemph_fac2 = shr( preemph_fac, 1 );
3065 3802 : ptr = fftBuffer;
3066 3802 : *ptr = Mpy_32_16_1( *ptr, sub( tmp, preemph_fac2 ) );
3067 3802 : move32();
3068 3802 : ptr += 2;
3069 281392 : FOR( i = 1; i < fftlen8; i++ )
3070 : {
3071 277590 : move32();
3072 277590 : *ptr = Mpy_32_16_1( *ptr, sub( tmp, mult_r( preemph_fac2, add( shr( table[i - 1].v.re, 1 ), shr( table[i].v.re, 1 ) ) ) ) );
3073 277590 : ptr += 2;
3074 : }
3075 3802 : move32();
3076 3802 : *ptr = Mpy_32_16_1( *ptr, sub( tmp, mult_r( preemph_fac2, add( shr( table[fftlen8 - 1].v.re, 1 ), shr( table[fftlen8 - 1].v.im, 1 ) ) ) ) );
3077 3802 : ptr += 2;
3078 281392 : FOR( i = 1; i < fftlen8; i++ )
3079 : {
3080 277590 : move32();
3081 277590 : *ptr = Mpy_32_16_1( *ptr, sub( tmp, mult_r( preemph_fac2, add( shr( table[fftlen8 - i - 1].v.im, 1 ), shr( table[fftlen8 - i].v.im, 1 ) ) ) ) );
3082 277590 : ptr += 2;
3083 : }
3084 3802 : move32();
3085 3802 : *ptr = Mpy_32_16_1( *ptr, tmp );
3086 3802 : ptr += 2;
3087 281392 : FOR( i = 1; i < fftlen8; i++ )
3088 : {
3089 277590 : move32();
3090 277590 : *ptr = Mpy_32_16_1( *ptr, add( tmp, mult_r( preemph_fac2, add( shr( table[i - 1].v.im, 1 ), shr( table[i].v.im, 1 ) ) ) ) );
3091 277590 : ptr += 2;
3092 : }
3093 3802 : move32();
3094 3802 : *ptr = Mpy_32_16_1( *ptr, add( tmp, mult_r( preemph_fac2, add( shr( table[fftlen8 - 1].v.re, 1 ), shr( table[fftlen8 - 1].v.im, 1 ) ) ) ) );
3095 3802 : ptr += 2;
3096 281392 : FOR( i = 1; i < fftlen8; i++ )
3097 : {
3098 277590 : move32();
3099 277590 : *ptr = Mpy_32_16_1( *ptr, add( tmp, mult_r( preemph_fac2, add( shr( table[fftlen8 - i - 1].v.re, 1 ), shr( table[fftlen8 - i].v.re, 1 ) ) ) ) );
3100 277590 : ptr += 2;
3101 : }
3102 3802 : move32();
3103 3802 : fftBuffer[1] = Mpy_32_16_1( fftBuffer[1], add( tmp, preemph_fac2 ) );
3104 3802 : maxVal = 0;
3105 3802 : move32();
3106 2254938 : FOR( i = 0; i < fftlen; i++ )
3107 : {
3108 2251136 : maxVal = L_max( maxVal, L_abs( fftBuffer[i] ) );
3109 : }
3110 3802 : s2 = norm_l( maxVal );
3111 2254938 : FOR( i = 0; i < fftlen; i++ )
3112 : {
3113 2251136 : fftBuffer[i] = L_shl( fftBuffer[i], s2 );
3114 2251136 : move32();
3115 : }
3116 :
3117 : /* Autocorrelation */
3118 :
3119 3802 : BASOP_rfft( fftBuffer, fftlen, &scale, 1 );
3120 :
3121 3802 : s3 = getScaleFactor32( fftBuffer, add( lpcorder, 1 ) );
3122 :
3123 68436 : FOR( i = 0; i <= lpcorder; i++ )
3124 : {
3125 64634 : r[i] = L_shl( fftBuffer[i], s3 );
3126 64634 : move32();
3127 : }
3128 :
3129 3802 : r[0] = Mpy_32_32( r[0], 1074278656l /*1.0005f Q30*/ );
3130 3802 : move32();
3131 64634 : FOR( i = 1; i <= lpcorder; i++ )
3132 : {
3133 60832 : r[i] = Mpy_32_32( r[i], 1073741824l /*1.f Q30*/ );
3134 60832 : move32();
3135 : }
3136 3802 : s3 = getScaleFactor32( r, add( lpcorder, 1 ) );
3137 :
3138 68436 : FOR( i = 0; i <= lpcorder; i++ )
3139 : {
3140 64634 : r[i] = L_shl( r[i], s3 );
3141 64634 : move32();
3142 : }
3143 :
3144 68436 : FOR( i = 0; i <= lpcorder; i++ )
3145 : {
3146 64634 : L_Extract( r[i], &r_h[i], &r_l[i] );
3147 : }
3148 :
3149 : /* LPC */
3150 :
3151 3802 : E_LPC_lev_dur( r_h, r_l, A, NULL, lpcorder, NULL );
3152 3802 : }
3153 :
3154 : /*
3155 : msvq_decoder
3156 :
3157 : Parameters:
3158 :
3159 : cb i : Codebook (indexed cb[stages][levels][p]) format Q9.7
3160 : stages i : Number of stages
3161 : N i : Vector dimension
3162 : maxN i : Codebook vector dimension
3163 : Idx o : Indices
3164 : uq[] i : Quantized vector format Q9.7
3165 :
3166 :
3167 : Function:
3168 : multi stage vector dequantisation
3169 :
3170 : Returns:
3171 : void
3172 : */
3173 0 : void msvq_decoder(
3174 : const Word16 *const cb[], /* i : Codebook (indexed cb[*stages][levels][p]) */
3175 : const Word16 stages, /* i : Number of stages */
3176 : const Word16 N, /* i : Vector dimension */
3177 : const Word16 maxN, /* i : Codebook vector dimension */
3178 : const Word16 Idx[], /* i : Indices */
3179 : Word16 *uq /* o : quantized vector */
3180 : )
3181 : {
3182 : Word16 s, i, offset;
3183 :
3184 : // PMT("msvq_decoder Not verified")
3185 :
3186 0 : offset = i_mult( Idx[0], maxN );
3187 0 : FOR( i = 0; i < N; i++ )
3188 : {
3189 0 : uq[i] = cb[0][offset + i];
3190 0 : move16();
3191 : }
3192 :
3193 0 : FOR( s = 1; s < stages; s++ )
3194 : {
3195 0 : offset = i_mult( Idx[s], maxN );
3196 :
3197 0 : FOR( i = 0; i < N; i++ )
3198 : {
3199 0 : uq[i] = add( uq[i], cb[s][offset + i] );
3200 0 : move16();
3201 : }
3202 : }
3203 0 : }
3204 : /*-------------------------------------------------------------------
3205 : * FdCng_exc()
3206 : *
3207 : * Generate FD-CNG as LP excitation
3208 : *-------------------------------------------------------------------*/
3209 :
3210 18381 : void FdCng_exc(
3211 : HANDLE_FD_CNG_COM hFdCngCom,
3212 : Word16 *CNG_mode,
3213 : Word16 L_frame,
3214 : Word16 *lsp_old,
3215 : Word16 first_CNG,
3216 : Word16 *lspCNG,
3217 : Word16 *Aq, /* o: LPC coeffs Q12*/
3218 : Word16 *lsp_new, /* o: lsp Q15 */
3219 : Word16 *lsf_new, /* o: lsf Qlog2(2.56) */
3220 : Word16 *exc, /* o: LP excitation Q12 */
3221 : Word16 *exc2, /* o: LP excitation Q12 */
3222 : Word16 *bwe_exc /* o: LP excitation for BWE Q12*/
3223 : )
3224 : {
3225 : Word16 i;
3226 :
3227 18381 : *CNG_mode = -1;
3228 18381 : move16();
3229 :
3230 103435 : FOR( i = 0; i < L_frame / L_SUBFR; i++ )
3231 : {
3232 85054 : Copy( hFdCngCom->A_cng, Aq + i * ( M + 1 ), M + 1 );
3233 : }
3234 :
3235 18381 : E_LPC_a_lsp_conversion( Aq, lsp_new, lsp_old, M );
3236 :
3237 18381 : IF( first_CNG == 0 )
3238 : {
3239 7447 : Copy( lsp_old, lspCNG, M );
3240 : }
3241 312477 : FOR( i = 0; i < M; i++ )
3242 : {
3243 : /* AR low-pass filter */
3244 294096 : lspCNG[i] = mac_r( L_mult( CNG_ISF_FACT_FX, lspCNG[i] ), 32768 - CNG_ISF_FACT_FX, lsp_new[i] );
3245 294096 : move16(); /* Q15 (15+15+1-16) */
3246 : }
3247 :
3248 18381 : IF( EQ_16( L_frame, L_FRAME16k ) )
3249 : {
3250 11530 : lsp2lsf_fx( lsp_new, lsf_new, M, INT_FS_16k_FX );
3251 : }
3252 : ELSE
3253 : {
3254 6851 : E_LPC_lsp_lsf_conversion( lsp_new, lsf_new, M );
3255 : }
3256 18381 : Copy( hFdCngCom->exc_cng, exc, L_frame );
3257 18381 : Copy( hFdCngCom->exc_cng, exc2, L_frame );
3258 18381 : IF( bwe_exc != NULL )
3259 : {
3260 11158 : IF( EQ_16( L_frame, L_FRAME ) )
3261 : {
3262 6851 : interp_code_5over2_fx( exc2, bwe_exc, L_frame );
3263 : }
3264 : ELSE
3265 : {
3266 4307 : interp_code_4over2_fx( exc2, bwe_exc, L_frame );
3267 : }
3268 : }
3269 18381 : }
3270 :
3271 : /*-------------------------------------------------------------------
3272 : * SynthesisSTFT_fx()
3273 : *
3274 : * STFT synthesis filterbank
3275 : *-------------------------------------------------------------------*/
3276 :
3277 59265 : void SynthesisSTFT_fx(
3278 : Word32 *fftBuffer,
3279 : /* i : FFT bins */ // Q15
3280 : Word32 *timeDomainOutput,
3281 : Word32 *olapBuffer, // Qin
3282 : const Word16 *olapWin,
3283 : const Word16 tcx_transition,
3284 : HANDLE_FD_CNG_COM hFdCngCom, /* i/o: FD_CNG structure containing all buffers and variables */
3285 : const Word16 element_mode, /* i : element mode */
3286 : const Word16 nchan_out /* i : number of output channels */
3287 : )
3288 : {
3289 : Word16 i;
3290 : Word32 buf_fx[M + 1 + 320], tmp_fx;
3291 :
3292 : /* Perform IFFT */
3293 59265 : RFFTN_fx( fftBuffer, hFdCngCom->fftSineTab_fx, hFdCngCom->fftlen, 1 );
3294 :
3295 : /* Handle overlap in P/S domain for stereo */
3296 59265 : test();
3297 59265 : test();
3298 59265 : IF( ( EQ_16( element_mode, IVAS_CPE_TD ) || EQ_16( element_mode, IVAS_CPE_DFT ) ) && EQ_16( nchan_out, 2 ) )
3299 : {
3300 890 : Copy32( olapBuffer + sub( i_mult( 3, shr( hFdCngCom->frameSize, 2 ) ), ( M + 1 ) ), buf_fx, add( hFdCngCom->frameSize, M + 1 ) ); // Qin
3301 890 : set32_fx( olapBuffer, 0, hFdCngCom->fftlen );
3302 : }
3303 : ELSE
3304 : {
3305 58375 : Copy32( olapBuffer + hFdCngCom->frameSize, olapBuffer, hFdCngCom->frameSize );
3306 58375 : set32_fx( olapBuffer + hFdCngCom->frameSize, 0, hFdCngCom->frameSize ); /*olapBuffer, fftBuffer, olapWin*/
3307 : }
3308 :
3309 59265 : IF( tcx_transition )
3310 : {
3311 0 : FOR( i = 0; i < ( 5 * ( hFdCngCom->frameSize >> 2 ) ); i++ )
3312 : {
3313 0 : olapBuffer[i] = fftBuffer[i]; // Q15
3314 0 : move32();
3315 : }
3316 : }
3317 : ELSE
3318 : {
3319 8522689 : FOR( i = hFdCngCom->frameSize / 4; i < ( 3 * ( hFdCngCom->frameSize >> 2 ) ); i++ )
3320 : {
3321 8463424 : olapBuffer[i] = L_add( olapBuffer[i], Mpy_32_16_1( fftBuffer[i], olapWin[( i - ( hFdCngCom->frameSize >> 2 ) )] ) );
3322 8463424 : move32();
3323 : }
3324 8522689 : FOR( ; i < ( 5 * ( hFdCngCom->frameSize >> 2 ) ); i++ )
3325 : {
3326 8463424 : olapBuffer[i] = fftBuffer[i];
3327 8463424 : move32();
3328 : }
3329 : }
3330 8522689 : FOR( ; i < ( 7 * ( hFdCngCom->frameSize >> 2 ) ); i++ )
3331 : {
3332 8463424 : olapBuffer[i] = Mpy_32_16_1( fftBuffer[i], olapWin[( i - ( 3 * ( hFdCngCom->frameSize >> 2 ) ) )] );
3333 8463424 : move32();
3334 : }
3335 :
3336 4290977 : FOR( ; i < hFdCngCom->fftlen; i++ )
3337 : {
3338 4231712 : olapBuffer[i] = 0;
3339 4231712 : move32();
3340 : }
3341 :
3342 59265 : Word32 fftScale = 0;
3343 59265 : SWITCH( hFdCngCom->fftlen )
3344 : {
3345 27422 : case 640:
3346 27422 : fftScale = FFT_SCALING_640;
3347 27422 : move32();
3348 27422 : BREAK;
3349 31843 : case 512:
3350 31843 : fftScale = FFT_SCALING_512;
3351 31843 : move32();
3352 31843 : BREAK;
3353 0 : default:
3354 0 : assert( !"Not supported FFT length!" );
3355 : }
3356 : /* Get time-domain signal */
3357 : // v_multc(olapBuffer + hFdCngCom->frameSize / 4, (float)(hFdCngCom->fftlen / 2), timeDomainOutput, hFdCngCom->frameSize);
3358 59265 : v_multc_fixed( olapBuffer + hFdCngCom->frameSize / 4, fftScale, timeDomainOutput, hFdCngCom->frameSize ); // Q_in - 9
3359 : /* Get excitation */
3360 59265 : test();
3361 59265 : test();
3362 59265 : IF( ( EQ_16( element_mode, IVAS_CPE_TD ) || EQ_16( element_mode, IVAS_CPE_DFT ) ) && EQ_16( nchan_out, 2 ) )
3363 : {
3364 114810 : FOR( i = 0; i < hFdCngCom->frameSize / 2; i++ )
3365 : {
3366 113920 : buf_fx[i + ( M + 1 )] = L_add( buf_fx[i + ( M + 1 )], olapBuffer[( i + ( hFdCngCom->frameSize >> 2 ) )] );
3367 113920 : move32();
3368 : }
3369 : // v_multc(buf, (float)(hFdCngCom->fftlen / 2), buf, M + 1 + hFdCngCom->frameSize);
3370 890 : v_multc_fixed( buf_fx, fftScale, buf_fx, add( M + 1, hFdCngCom->frameSize ) );
3371 : }
3372 : ELSE
3373 : {
3374 : // v_multc(olapBuffer + hFdCngCom->frameSize / 4 - (M + 1), (float)(hFdCngCom->fftlen / 2), buf, M + 1 + hFdCngCom->frameSize);
3375 58375 : v_multc_fixed( olapBuffer + sub( shr( hFdCngCom->frameSize, 2 ), ( M + 1 ) ), fftScale, buf_fx, add( M + 1, hFdCngCom->frameSize ) );
3376 : }
3377 :
3378 59265 : tmp_fx = buf_fx[0];
3379 59265 : move32();
3380 : // preemph(buf + 1, PREEMPH_FAC_FLT, M + hFdCngCom->frameSize, &tmp);
3381 59265 : preemph_ivas_fx( buf_fx + 1, PREEMPH_FAC, add( M, hFdCngCom->frameSize ), &tmp_fx );
3382 : // residu(hFdCngCom->A_cng_flt, M, buf + 1 + M, hFdCngCom->exc_cng_flt, hFdCngCom->frameSize);
3383 :
3384 : // residu_ivas_fx( hFdCngCom->A_cng, Q13, M, buf_fx + 1 + M, hFdCngCom->exc_cng_32fx, hFdCngCom->frameSize );
3385 59265 : residu_ivas_fx( hFdCngCom->A_cng, ( 15 - norm_s( hFdCngCom->A_cng[0] - 1 ) ), M, buf_fx + 1 + M, hFdCngCom->exc_cng_32fx, hFdCngCom->frameSize );
3386 :
3387 59265 : return;
3388 : }
3389 :
3390 4395 : void SynthesisSTFT_dirac_fx(
3391 : Word32 *fftBuffer,
3392 : /* i : FFT bins */ // hFdCngCom->fftBuffer_exp
3393 : Word32 *timeDomainOutput,
3394 : Word32 *olapBuffer, // Q_in
3395 : const Word16 *olapWin,
3396 : const Word16 samples_out,
3397 : HANDLE_FD_CNG_COM hFdCngCom /* i/o: FD_CNG structure containing all buffers and variables */
3398 : )
3399 : {
3400 : Word16 i;
3401 : Word32 buf[M + 1 + 320], tmp;
3402 :
3403 : /* Perform IFFT */
3404 4395 : RFFTN_fx( fftBuffer, hFdCngCom->fftSineTab_fx, hFdCngCom->fftlen, 1 );
3405 :
3406 : /* Handle overlap in P/S domain for stereo */
3407 4395 : Copy32( olapBuffer + hFdCngCom->frameSize, olapBuffer, hFdCngCom->frameSize );
3408 4395 : set32_fx( olapBuffer + hFdCngCom->frameSize, 0, hFdCngCom->frameSize ); /*olapBuffer, fftBuffer, olapWin*/
3409 :
3410 650411 : FOR( i = ( hFdCngCom->frameSize >> 2 ); i < ( 3 * ( hFdCngCom->frameSize >> 2 ) ); i++ )
3411 : {
3412 646016 : olapBuffer[i] = L_add( olapBuffer[i], Mpy_32_16_1( fftBuffer[i], olapWin[( i - ( hFdCngCom->frameSize >> 2 ) )] ) );
3413 646016 : move32();
3414 : }
3415 650411 : FOR( ; i < ( 5 * ( hFdCngCom->frameSize >> 2 ) ); i++ )
3416 : {
3417 646016 : olapBuffer[i] = fftBuffer[i];
3418 646016 : move32();
3419 : }
3420 :
3421 650411 : FOR( ; i < ( 7 * ( hFdCngCom->frameSize >> 2 ) ); i++ )
3422 : {
3423 646016 : olapBuffer[i] = fftBuffer[i];
3424 646016 : move32();
3425 : }
3426 :
3427 327403 : FOR( ; i < hFdCngCom->fftlen; i++ )
3428 : {
3429 323008 : olapBuffer[i] = 0;
3430 323008 : move32();
3431 : }
3432 :
3433 4395 : Word32 fftScale = 0;
3434 4395 : move32();
3435 4395 : SWITCH( hFdCngCom->fftlen )
3436 : {
3437 2608 : case 640:
3438 2608 : fftScale = FFT_SCALING_640;
3439 2608 : move32();
3440 2608 : break;
3441 1787 : case 512:
3442 1787 : fftScale = FFT_SCALING_512;
3443 1787 : move32();
3444 1787 : break;
3445 0 : default:
3446 0 : assert( !"Not supported FFT length!" );
3447 : }
3448 :
3449 : /* Get time-domain signal */
3450 4395 : v_multc_fixed( olapBuffer + shr( hFdCngCom->frameSize, 2 ), fftScale, timeDomainOutput, samples_out ); // Q_in - 9
3451 :
3452 : /* Get excitation */
3453 4395 : v_multc_fixed( olapBuffer + sub( shr( hFdCngCom->frameSize, 2 ), ( M + 1 ) ), fftScale, buf, add( M + 1, hFdCngCom->frameSize ) );
3454 4395 : tmp = buf[0];
3455 4395 : move32();
3456 4395 : preemph_ivas_fx( buf + 1, PREEMPH_FAC, M + hFdCngCom->frameSize, &tmp );
3457 : // residu_ivas_fx( hFdCngCom->A_cng, Q13, M, buf + 1 + M, hFdCngCom->exc_cng_32fx, hFdCngCom->frameSize );
3458 4395 : residu_ivas_fx( hFdCngCom->A_cng, sub( 15, norm_s( hFdCngCom->A_cng[0] - 1 ) ), M, buf + 1 + M, hFdCngCom->exc_cng_32fx, hFdCngCom->frameSize );
3459 :
3460 : /* update and window olapBuf if we have a output frame that is shorter than the default frame size...*/
3461 4395 : IF( LT_16( samples_out, hFdCngCom->frameSize ) )
3462 : {
3463 0 : Copy32( olapBuffer + samples_out, olapBuffer + hFdCngCom->frameSize, i_mult( 3, shr( hFdCngCom->frameSize, 2 ) ) );
3464 : }
3465 650411 : FOR( i = ( 5 * ( hFdCngCom->frameSize >> 2 ) ); i < ( 7 * ( hFdCngCom->frameSize >> 2 ) ); i++ )
3466 : {
3467 646016 : olapBuffer[i] = Mpy_32_16_1( olapBuffer[i], olapWin[( i - ( 3 * ( hFdCngCom->frameSize >> 2 ) ) )] );
3468 646016 : move32();
3469 : }
3470 :
3471 4395 : return;
3472 : }
3473 :
3474 46665862 : Word32 rand_gauss_fx(
3475 : Word32 *x,
3476 : Word16 *seed,
3477 : Word16 q )
3478 : {
3479 : Word32 temp;
3480 :
3481 46665862 : temp = own_random( seed );
3482 46665862 : temp = L_add( temp, own_random( seed ) );
3483 46665862 : temp = L_add( temp, own_random( seed ) );
3484 46665862 : temp = L_shr( temp, sub( 15, q ) );
3485 :
3486 46665862 : *x = temp;
3487 46665862 : move32();
3488 :
3489 46665862 : return temp;
3490 : }
3491 :
3492 : /*-------------------------------------------------------------------
3493 : * rand_gauss_fix()
3494 : *
3495 : * Random generator with Gaussian distribution with mean 0 and std 1
3496 : *-------------------------------------------------------------------*/
3497 :
3498 50857104 : Word16 rand_gauss_fix(
3499 : Word16 *x,
3500 : Word16 *seed )
3501 : {
3502 : Word32 temp;
3503 :
3504 50857104 : temp = shr( own_random( seed ), Q2 );
3505 50857104 : temp = L_add( temp, shr( own_random( seed ), Q2 ) );
3506 50857104 : temp = L_add( temp, shr( own_random( seed ), Q2 ) );
3507 :
3508 50857104 : *x = (Word16) temp;
3509 50857104 : move32();
3510 :
3511 50857104 : return (Word16) temp;
3512 : }
|
