Line data Source code
1 : /******************************************************************************************************
2 :
3 : (C) 2022-2025 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB,
4 : Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,
5 : Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,
6 : Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other
7 : contributors to this repository. All Rights Reserved.
8 :
9 : This software is protected by copyright law and by international treaties.
10 : The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB,
11 : Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,
12 : Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,
13 : Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other
14 : contributors to this repository retain full ownership rights in their respective contributions in
15 : the software. This notice grants no license of any kind, including but not limited to patent
16 : license, nor is any license granted by implication, estoppel or otherwise.
17 :
18 : Contributors are required to enter into the IVAS codec Public Collaboration agreement before making
19 : contributions.
20 :
21 : This software is provided "AS IS", without any express or implied warranties. The software is in the
22 : development stage. It is intended exclusively for experts who have experience with such software and
23 : solely for the purpose of inspection. All implied warranties of non-infringement, merchantability
24 : and fitness for a particular purpose are hereby disclaimed and excluded.
25 :
26 : Any dispute, controversy or claim arising under or in relation to providing this software shall be
27 : submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in
28 : accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and
29 : the United Nations Convention on Contracts on the International Sales of Goods.
30 :
31 : *******************************************************************************************************/
32 :
33 : /*====================================================================================
34 : EVS Codec 3GPP TS26.443 Nov 04, 2021. Version 12.14.0 / 13.10.0 / 14.6.0 / 15.4.0 / 16.3.0
35 : ====================================================================================*/
36 :
37 : #include <stdint.h>
38 : #include "options.h"
39 : #ifdef DEBUGGING
40 : #include "debug.h"
41 : #endif
42 : #include <math.h>
43 : #include "rom_dec.h"
44 : #include "rom_com.h"
45 : #include "cnst.h"
46 : #include "prot.h"
47 : #include "wmc_auto.h"
48 :
49 : /*---------------------------------------------------------------------*
50 : * Local constants
51 : *---------------------------------------------------------------------*/
52 :
53 : #define FEC_MAX 512
54 : #define FEC_NB_PULSE_MAX 20
55 : #define FEC_FFT_MAX_SIZE 512
56 : #define FEC_DCIM_FILT_SIZE_MAX 60
57 :
58 : #define PHASE_DITH ( PI2 )
59 :
60 : #define DELTA_CORR 6 /* Range for phase correction around peak */
61 : #define THRESH_TR_dB 10.0f
62 : #define THRESH_TR_LIN (float) pow( 10.0f, THRESH_TR_dB / 10.0f )
63 : #define THRESH_TR_LIN_INV (float) pow( 10.0f, -THRESH_TR_dB / 10.0f )
64 : #define MAX_INCREASE_GRPOW 0.0f /* maximum amplification in case of transients */
65 : #define MAX_INCREASE_GRPOW_LIN (float) pow( 10.0f, MAX_INCREASE_GRPOW / 10.0f )
66 :
67 : #define PHASE_DITH_SCALE (float) pow( 2.0, -16.0 ) /* for scaling random short values to +/- pi */
68 :
69 : #define BURST_PHDITH_THRESH ( 4 - 1 ) /* speech start phase dither with <burst_phdith_thresh> losses in a row */
70 : #define BURST_PHDITH_RAMPUP_LEN 2 /* speech ramp up degree of phase dither over a length of <burst_phdith_rampup_len> frames */
71 : #define BURST_ATT_THRESH ( 3 - 1 ) /* speech start attenuate with <burst_att_thresh> losses in a row */
72 : #define ATT_PER_FRAME 4 /* speech attenuation in dB */
73 : #define BETA_MUTE_THR 10 /* time threshold to start beta-noise attenuation */
74 : #define BETA_MUTE_FAC 0.5f /* attenuation factor per additional bad frame */
75 :
76 : #define LGW32k 7
77 : #define LGW16k 6
78 : #define LGW48k LGW32k + 1 /* Use the same frequency groups as for SWB + 1 */
79 :
80 : #define L_TRANA_LOG32k 8
81 : #define L_TRANA_LOG16k 7
82 :
83 : #define DELTA_CORR_F0_INT 2 /* Constant controls the bin range where Jacobsen is used */
84 : #define ST_PFIND_SENS 0.93f /* peakfinder sensitivity */
85 : #define L_PROT_NS 32000000L /* Prototype frame length in nanoseconds (32 ms) */
86 : #define PH_ECU_CORR_LIMIT 0.85f /* Correlation limit for IVAS Phase ECU activation */
87 : #define PH_ECU_N_LIMIT 56 /* fec_alg analysis frame limit for IVAS Phase ECU activation */
88 : #define PFIND_SENS 0.97f /* peakfinder sensitivity */
89 :
90 : /*---------------------------------------------------------------------*
91 : * Local functions
92 : *---------------------------------------------------------------------*/
93 :
94 : static int16_t rand_phase( const int16_t seed, float *sin_F, float *cos_F );
95 : static float imax2_jacobsen_mag( const float *y_re, const float *y_im );
96 :
97 : /*-------------------------------------------------------------------*
98 : * mult_rev2()
99 : *
100 : * Multiplication of two vectors second vector is multiplied in reverse order
101 : *-------------------------------------------------------------------*/
102 :
103 376 : static void mult_rev2(
104 : const float x1[], /* i : Input vector 1 */
105 : const float x2[], /* i : Input vector 2 */
106 : float y[], /* o : Output vector that contains vector 1 .* vector 2 */
107 : const int16_t N /* i : Vector length */
108 : )
109 : {
110 : int16_t i, j;
111 :
112 65912 : for ( i = 0, j = N - 1; i < N; i++, j-- )
113 : {
114 65536 : y[i] = x1[i] * x2[j];
115 : }
116 :
117 376 : return;
118 : }
119 :
120 :
121 : /*-------------------------------------------------------------------*
122 : * fft_spec2()
123 : *
124 : * Square magnitude of fft spectrum
125 : *-------------------------------------------------------------------*/
126 :
127 564 : static void fft_spec2(
128 : float x[], /* i/o: Input vector: complex spectrum -> square magnitude spectrum */
129 : const int16_t N /* i : Vector length */
130 : )
131 : {
132 : int16_t i, j;
133 :
134 196608 : for ( i = 1, j = N - 1; i < N / 2; i++, j-- )
135 : {
136 196044 : x[i] = x[i] * x[i] + x[j] * x[j];
137 : }
138 :
139 564 : x[0] *= x[0];
140 564 : x[N / 2] *= x[N / 2];
141 :
142 564 : return;
143 : }
144 :
145 : /*------------------------------------------------------------------*
146 : * rand_phase()
147 : *
148 : * randomized phase in form of sin and cos components
149 : *------------------------------------------------------------------*/
150 :
151 : /*! r: Updated seed from RNG */
152 57050 : static int16_t rand_phase(
153 : const int16_t seed, /* i : RNG seed */
154 : float *sin_F, /* o : random phase sin value */
155 : float *cos_F /* o : random phase cos value */
156 : )
157 : {
158 57050 : const float *sincos = sincos_t_ext + 128;
159 57050 : int16_t seed2 = seed;
160 57050 : own_random( &seed2 );
161 :
162 57050 : if ( seed2 & 0x40 )
163 : {
164 28524 : *sin_F = sincos[seed2 >> 8];
165 : }
166 : else
167 : {
168 28526 : *sin_F = -sincos[seed2 >> 8];
169 : }
170 :
171 57050 : if ( seed2 & 0x80 )
172 : {
173 28510 : *cos_F = sincos[-( seed2 >> 8 )];
174 : }
175 : else
176 : {
177 28540 : *cos_F = -sincos[-( seed2 >> 8 )];
178 : }
179 :
180 57050 : return seed2;
181 : }
182 :
183 : /*-----------------------------------------------------------------------------
184 : * imax2_jacobsen_mag()
185 : *
186 : * refine peak interpolation using jacobsen and periodic speca ana windows
187 : *----------------------------------------------------------------------------*/
188 :
189 : /*! r: The location, relative to the middle of the 3 given data point, of the maximum. (Q15)*/
190 2700 : float imax2_jacobsen_mag(
191 : const float *y_re, /* i : The 3 given data points. real part order -1 0 1 */
192 : const float *y_im /* i : The 3 given data points. imag part order 1 0 -1 (from FFT) */
193 : )
194 : {
195 : float posi;
196 : const float *pY;
197 : float y_m1_re, y_0_re, y_p1_re;
198 : float y_m1_im, y_0_im, y_p1_im;
199 : float N_re, N_im;
200 : float D_re, D_im;
201 : float numer, denom;
202 :
203 : /* Jacobsen estimates peak offset relative y_0 using
204 : * X_m1 - X_p1
205 : * d = REAL ( ------------------- ) * c_jacob
206 : * 2*X_0 - X_m1 -Xp1
207 : *
208 : * Where c_jacob is a window dependent constant
209 : */
210 : #define C_JACOB 1.1453f /* % assume 0.1875 hammrect window 'symmetric' */
211 :
212 : /* Get the bin parameters into variables */
213 2700 : pY = y_re;
214 2700 : y_m1_re = *pY++;
215 2700 : y_0_re = *pY++;
216 2700 : y_p1_re = *pY++;
217 :
218 : /* Same for imaginary parts - note reverse order from FFT */
219 2700 : pY = y_im;
220 2700 : y_p1_im = *pY++;
221 2700 : y_0_im = *pY++;
222 2700 : y_m1_im = *pY++;
223 :
224 : /* prepare numerator real and imaginary parts*/
225 2700 : N_re = y_m1_re - y_p1_re;
226 2700 : N_im = y_m1_im - y_p1_im;
227 :
228 : /* prepare denominator real and imaginary parts */
229 :
230 2700 : D_re = 2 * y_0_re - y_m1_re - y_p1_re;
231 2700 : D_im = 2 * y_0_im - y_m1_im - y_p1_im;
232 :
233 : /* REAL part of complex division */
234 2700 : numer = N_re * D_re + N_im * D_im;
235 2700 : denom = D_re * D_re + D_im * D_im;
236 :
237 2700 : test();
238 2700 : if ( numer != 0 && denom != 0 )
239 : {
240 2700 : posi = numer / denom * C_JACOB;
241 : }
242 : else
243 : {
244 0 : posi = 0; /* flat top, division is not possible choose center freq */
245 : }
246 :
247 2700 : return posi;
248 : }
249 :
250 :
251 : /*------------------------------------------------------------------*
252 : * trans_ana()
253 : *
254 : * Transient analysis
255 : *------------------------------------------------------------------*/
256 :
257 276 : static void trans_ana(
258 : const float *xfp, /* i : Input signal */
259 : float *mag_chg, /* i/o: Magnitude modification */
260 : float *ph_dith, /* i/o: Phase dither */
261 : float *mag_chg_1st, /* i/o: per band magnitude modifier for transients */
262 : const int16_t output_frame, /* i : Frame length */
263 : const int16_t time_offs, /* i : Time offset */
264 : const float est_mus_content, /* i : 0.0=speech_like ... 1.0=Music (==st->env_stab ) */
265 : const int16_t last_fec, /* i : signal that previous frame was concealed with fec_alg*/
266 : const int16_t element_mode, /* i : element_mode req to handle EVS_MONO specific BE path*/
267 : float *alpha, /* o : Magnitude modification factors for fade to average */
268 : float *beta, /* o : Magnitude modification factors for fade to average */
269 : float *beta_mute, /* o : Factor for long-term mute */
270 : float Xavg[LGW_MAX] /* o : Frequency group average gain to fade to */
271 : )
272 : {
273 : const float *w_hamm;
274 : float grp_pow_chg, att_val, att_degree;
275 : float xfp_left[L_TRANA48k], xfp_right[L_TRANA48k];
276 : float gr_pow_left[LGW_MAX], gr_pow_right[LGW_MAX];
277 : const float *xfp_;
278 276 : int16_t Ltrana, Ltrana_2, Lprot, LtranaLog = 0, Lgw, k, burst_len;
279 : int16_t att_always[LGW_MAX]; /* fixed attenuation per frequency group if set to 1*/
280 : int16_t burst_phdith_thresh;
281 : int16_t burst_att_thresh;
282 : float att_per_frame;
283 : int16_t tr_dec[LGW_MAX];
284 :
285 : /* check burst error */
286 276 : burst_len = time_offs / output_frame + 1;
287 :
288 276 : set_s( att_always, 0, LGW_MAX );
289 276 : *ph_dith = 0.0f;
290 :
291 : /* softly shift attenuation just a bit later for estimated "stable" music_content */
292 276 : burst_phdith_thresh = BURST_PHDITH_THRESH + (int16_t) ( est_mus_content * 1.0f + 0.5f );
293 276 : burst_att_thresh = BURST_ATT_THRESH + (int16_t) ( est_mus_content * 1.0f + 0.5f );
294 276 : att_per_frame = (float) ( ATT_PER_FRAME - (int16_t) ( est_mus_content * 1.0f + 0.5f ) ); /* only slighty less att for music */
295 276 : att_per_frame *= 0.1f;
296 :
297 276 : if ( burst_len > burst_phdith_thresh )
298 : {
299 : /* increase degree of dither */
300 72 : *ph_dith = PHASE_DITH * min( 1.0f, ( (float) burst_len - (float) burst_phdith_thresh ) / (float) BURST_PHDITH_RAMPUP_LEN );
301 : }
302 :
303 276 : att_degree = 0;
304 276 : if ( burst_len > burst_att_thresh )
305 : {
306 74 : set_s( att_always, 1, LGW_MAX );
307 :
308 : /* increase degree of attenuation */
309 74 : if ( burst_len - burst_att_thresh <= PH_ECU_MUTE_START )
310 : {
311 30 : att_degree = (float) ( burst_len - burst_att_thresh ) * att_per_frame;
312 : }
313 : else
314 : {
315 44 : att_degree = (float) PH_ECU_MUTE_START * att_per_frame + ( burst_len - burst_att_thresh - PH_ECU_MUTE_START ) * 6.0206f;
316 : }
317 : }
318 :
319 276 : Lprot = ( 2 * output_frame * 4 ) / 5; /* 4/5==1024/1280, keep mult within short */
320 276 : Ltrana = Lprot / QUOT_LPR_LTR;
321 276 : Ltrana_2 = Ltrana / 2;
322 :
323 276 : if ( output_frame == L_FRAME48k )
324 : {
325 221 : w_hamm = w_hamm48k_2;
326 221 : Lgw = LGW48k;
327 : }
328 55 : else if ( output_frame == L_FRAME32k )
329 : {
330 55 : w_hamm = w_hamm32k_2;
331 55 : LtranaLog = L_TRANA_LOG32k;
332 55 : Lgw = LGW32k;
333 : }
334 : else
335 : {
336 0 : w_hamm = w_hamm16k_2;
337 0 : LtranaLog = L_TRANA_LOG16k;
338 0 : Lgw = LGW16k;
339 : }
340 :
341 276 : if ( burst_len <= 1 || ( burst_len == 2 && last_fec ) )
342 : {
343 188 : set_f( alpha, 1.0f, LGW_MAX );
344 188 : set_f( beta, 0.0f, LGW_MAX );
345 188 : *beta_mute = BETA_MUTE_FAC_INI;
346 :
347 : /* apply hamming window */
348 188 : v_mult( xfp, w_hamm, xfp_left, Ltrana_2 );
349 188 : mult_rev2( xfp + Ltrana_2, w_hamm, xfp_left + Ltrana_2, Ltrana_2 );
350 :
351 188 : xfp_ = xfp + Lprot - Ltrana;
352 188 : v_mult( xfp_, w_hamm, xfp_right, Ltrana_2 );
353 188 : mult_rev2( xfp_ + Ltrana_2, w_hamm, xfp_right + Ltrana_2, Ltrana_2 );
354 :
355 : /* spectrum */
356 188 : if ( output_frame == L_FRAME48k )
357 : {
358 136 : fft3( xfp_left, xfp_left, Ltrana );
359 136 : fft3( xfp_right, xfp_right, Ltrana );
360 : }
361 : else
362 : {
363 52 : fft_rel( xfp_left, Ltrana, LtranaLog );
364 52 : fft_rel( xfp_right, Ltrana, LtranaLog );
365 : }
366 :
367 : /* square representation */
368 188 : fft_spec2( xfp_left, Ltrana );
369 188 : fft_spec2( xfp_right, Ltrana );
370 :
371 : /* band powers in frequency groups
372 : exclude bin at 0 and at EVS_PI from calculation */
373 188 : xfp_left[Ltrana_2] = 0.0f;
374 188 : xfp_right[Ltrana_2] = 0.0f;
375 : }
376 :
377 2429 : for ( k = 0; k < Lgw; k++ )
378 : {
379 2153 : if ( burst_len <= 1 || ( burst_len == 2 && last_fec ) )
380 : {
381 1452 : gr_pow_left[k] = sum_f( xfp_left + gw[k], gw[k + 1] - gw[k] );
382 1452 : gr_pow_right[k] = sum_f( xfp_right + gw[k], gw[k + 1] - gw[k] );
383 :
384 : /* check if transient in any of the bands */
385 1452 : gr_pow_left[k] += FLT_MIN; /* otherwise div by zero may occur */
386 1452 : gr_pow_right[k] += FLT_MIN;
387 :
388 1452 : Xavg[k] = (float) ( sqrt( 0.5f * ( gr_pow_left[k] + gr_pow_right[k] ) / (float) ( gw[k + 1] - gw[k] ) ) );
389 :
390 1452 : grp_pow_chg = gr_pow_right[k] / gr_pow_left[k];
391 :
392 : /* dither phase in case of transient */
393 : /* separate transition detection and application of forced burst dithering */
394 1452 : tr_dec[k] = ( grp_pow_chg > THRESH_TR_LIN ) || ( grp_pow_chg < THRESH_TR_LIN_INV );
395 :
396 : /* magnitude modification */
397 1452 : if ( tr_dec[k] || att_always[k] )
398 : {
399 68 : att_val = min( MAX_INCREASE_GRPOW_LIN, grp_pow_chg );
400 68 : att_val = (float) sqrt( att_val );
401 68 : mag_chg_1st[k] = att_val;
402 68 : mag_chg[k] = att_val;
403 : }
404 : else
405 : {
406 1384 : mag_chg_1st[k] = 1.0f;
407 1384 : mag_chg[k] = 1.0f;
408 : }
409 : }
410 : else
411 : {
412 701 : if ( burst_len < OFF_FRAMES_LIMIT )
413 : {
414 525 : mag_chg[k] = mag_chg_1st[k] * (float) pow( 10.0, -att_degree / 20.0 );
415 : }
416 : else
417 : {
418 176 : mag_chg[k] = 0;
419 : }
420 :
421 701 : if ( element_mode != EVS_MONO )
422 : {
423 701 : if ( k == 0 && burst_len > BETA_MUTE_THR ) /* beta_mute final long term attenuation adjusted only once per frame in the first sub-band, Ref Eq(184) in 26.447 */
424 : {
425 60 : *beta_mute *= BETA_MUTE_FAC;
426 : }
427 : }
428 : else
429 : {
430 0 : if ( burst_len > BETA_MUTE_THR ) /* legacy incorrect(too fast) EVS_MONO operation, still kept for BE. To be updated after EVS CR, ref Eq (184) in 26.447 */
431 : {
432 0 : *beta_mute *= BETA_MUTE_FAC;
433 : }
434 : }
435 :
436 701 : alpha[k] = mag_chg[k];
437 701 : beta[k] = (float) ( sqrt( 1.0f - SQR( alpha[k] ) ) * *beta_mute );
438 701 : if ( k >= LGW32k - 1 )
439 : {
440 173 : beta[k] *= 0.1f;
441 : }
442 528 : else if ( k >= LGW16k - 1 )
443 : {
444 88 : beta[k] *= 0.5f;
445 : }
446 : }
447 : }
448 :
449 276 : return;
450 : }
451 :
452 :
453 : /*------------------------------------------------------------------*
454 : * peakfinder()
455 : *
456 : * Peak-picking algorithm
457 : *------------------------------------------------------------------*/
458 :
459 874 : void peakfinder(
460 : const float *x0, /* i : vector from which the maxima will be found */
461 : const int16_t len0, /* i : length of input vector */
462 : int16_t *plocs, /* o : the indicies of the identified peaks in x0 */
463 : int16_t *cInd, /* o : number of identified peaks */
464 : const float sel, /* i : The amount above surrounding data for a peak to be identified */
465 : const int16_t endpoints /* i : Flag to include endpoints in peak search */
466 : )
467 : {
468 : float minMag, tempMag, leftMin;
469 : float dx0[L_PROT48k_2], x[L_PROT48k_2 + 1], peakMag[MAX_PLOCS];
470 874 : int16_t k, i, len, tempLoc = 0, foundPeak, ii, xInd;
471 : int16_t *ind, indarr[L_PROT48k_2 + 1], peakLoc[MAX_PLOCS];
472 :
473 874 : ind = indarr;
474 :
475 : /* Find derivative */
476 874 : v_sub( x0 + 1, x0, dx0, len0 - 1 );
477 :
478 : /* This is so we find the first of repeated values */
479 145666 : for ( i = 0; i < len0 - 1; i++ )
480 : {
481 144792 : if ( dx0[i] == 0.0f )
482 : {
483 6205 : dx0[i] = -1.0e-12f;
484 : }
485 : }
486 :
487 : /* Find where the derivative changes sign
488 : Include endpoints in potential peaks and valleys */
489 874 : k = 0;
490 :
491 874 : if ( endpoints )
492 : {
493 188 : x[k] = x0[0];
494 188 : ind[k++] = 0;
495 : }
496 :
497 144792 : for ( i = 1; i < len0 - 1; i++ )
498 : {
499 143918 : if ( dx0[i - 1] * dx0[i] < 0 )
500 : {
501 77297 : ind[k] = i;
502 77297 : x[k++] = x0[i];
503 : }
504 : }
505 :
506 874 : if ( endpoints )
507 : {
508 188 : ind[k] = len0 - 1;
509 188 : x[k++] = x0[len0 - 1];
510 : }
511 : /* x only has the peaks, valleys, and endpoints */
512 874 : len = k;
513 874 : minimum( x, len, &minMag );
514 :
515 874 : if ( ( len > 2 ) || ( !endpoints && ( len > 0 ) ) )
516 : {
517 : /* Set initial parameters for loop */
518 870 : tempMag = minMag;
519 870 : foundPeak = 0;
520 870 : leftMin = minMag;
521 :
522 870 : if ( endpoints )
523 : {
524 : /* Deal with first point a little differently since tacked it on
525 : Calculate the sign of the derivative since we taked the first point
526 : on it does not necessarily alternate like the rest. */
527 :
528 : /* The first point is larger or equal to the second */
529 188 : if ( x[0] >= x[1] )
530 : {
531 69 : ii = -1;
532 69 : if ( x[1] >= x[2] ) /* x[1] is not extremum -> overwrite with x[0] */
533 : {
534 0 : x[1] = x[0];
535 0 : ind[1] = ind[0];
536 0 : ind++;
537 0 : len--;
538 : }
539 : }
540 : else /* First point is smaller than the second */
541 : {
542 119 : ii = 0;
543 119 : if ( x[1] < x[2] ) /* x[1] is not extremum -> overwrite with x[0] */
544 : {
545 0 : x[1] = x[0];
546 0 : ind[1] = ind[0];
547 0 : ind++;
548 0 : len--;
549 : }
550 : }
551 : }
552 : else
553 : {
554 682 : ii = -1; /* First point is a peak */
555 682 : if ( len >= 2 )
556 : {
557 682 : if ( x[1] >= x[0] )
558 : {
559 119 : ii = 0; /* First point is a valley, skip it */
560 : }
561 : }
562 : }
563 870 : *cInd = 0;
564 :
565 : /* Loop through extrema which should be peaks and then valleys */
566 39360 : while ( ii < len - 1 )
567 : {
568 38945 : ii++; /* This is a peak */
569 :
570 : /*Reset peak finding if we had a peak and the next peak is bigger
571 : than the last or the left min was small enough to reset.*/
572 38945 : if ( foundPeak )
573 : {
574 4535 : tempMag = minMag;
575 4535 : foundPeak = 0;
576 : }
577 :
578 : /* Make sure we don't iterate past the length of our vector */
579 38945 : if ( ii == len - 1 )
580 : {
581 455 : break; /* We assign the last point differently out of the loop */
582 : }
583 :
584 : /* Found new peak that was larger than temp mag and selectivity larger
585 : than the minimum to its left. */
586 38490 : if ( ( x[ii] > tempMag ) && ( x[ii] > leftMin + sel ) )
587 : {
588 5071 : tempLoc = ii;
589 5071 : tempMag = x[ii];
590 : }
591 :
592 38490 : ii++; /* Move onto the valley */
593 :
594 : /* Come down at least sel from peak */
595 38490 : if ( !foundPeak && ( tempMag > sel + x[ii] ) )
596 : {
597 4652 : foundPeak = 1; /* We have found a peak */
598 4652 : leftMin = x[ii];
599 4652 : peakLoc[*cInd] = tempLoc; /* Add peak to index */
600 4652 : peakMag[*cInd] = tempMag;
601 4652 : ( *cInd )++;
602 : }
603 33838 : else if ( x[ii] < leftMin ) /* New left minimum */
604 : {
605 4648 : leftMin = x[ii];
606 : }
607 : }
608 :
609 : /* Check end point */
610 870 : if ( x[len - 1] > tempMag && x[len - 1] > leftMin + sel )
611 : {
612 138 : peakLoc[*cInd] = len - 1;
613 138 : peakMag[*cInd] = x[len - 1];
614 138 : ( *cInd )++;
615 : }
616 732 : else if ( !foundPeak && tempMag > minMag ) /* Check if we still need to add the last point */
617 : {
618 4 : peakLoc[*cInd] = tempLoc;
619 4 : peakMag[*cInd] = tempMag;
620 4 : ( *cInd )++;
621 : }
622 :
623 : /* Create output */
624 5664 : for ( i = 0; i < *cInd; i++ )
625 : {
626 4794 : plocs[i] = ind[peakLoc[i]];
627 : }
628 : }
629 : else
630 : {
631 4 : if ( endpoints )
632 : {
633 : /* This is a monotone function where an endpoint is the only peak */
634 0 : xInd = ( x[0] > x[1] ) ? 0 : 1;
635 0 : peakMag[0] = x[xInd];
636 0 : if ( peakMag[0] > minMag + sel )
637 : {
638 0 : plocs[0] = ind[xInd];
639 0 : *cInd = 1;
640 : }
641 : else
642 : {
643 0 : *cInd = 0;
644 : }
645 : }
646 : else
647 : {
648 : /* Input constant or all zeros -- no peaks found */
649 4 : *cInd = 0;
650 : }
651 : }
652 :
653 874 : return;
654 : }
655 :
656 :
657 : /*-------------------------------------------------------------------*
658 : * imax_pos()
659 : *
660 : * Get interpolated maximum position
661 : *-------------------------------------------------------------------*/
662 :
663 : /*! r: interpolated maximum position */
664 1721 : float imax_pos(
665 : const float *y /* i : Input vector for peak interpolation */
666 : )
667 : {
668 : float posi, y1, y2, y3, y3_y1, y2i;
669 : float ftmp_den1, ftmp_den2;
670 : /* Seek the extrema of the parabola P(x) defined by 3 consecutive points so that P([-1 0 1]) = [y1 y2 y3] */
671 1721 : y1 = y[0];
672 1721 : y2 = y[1];
673 1721 : y3 = y[2];
674 1721 : y3_y1 = y3 - y1;
675 1721 : ftmp_den1 = ( y1 + y3 - 2 * y2 );
676 1721 : ftmp_den2 = ( 4 * y2 - 2 * y1 - 2 * y3 );
677 :
678 1721 : if ( ftmp_den2 == 0.0f || ftmp_den1 == 0.0f )
679 : {
680 0 : return ( 0.0f ); /* early exit with left-most value */
681 : }
682 :
683 1721 : y2i = -0.125f * SQR( y3_y1 ) / ( ftmp_den1 ) + y2;
684 : /* their corresponding normalized locations */
685 1721 : posi = y3_y1 / ( ftmp_den2 );
686 : /* Interpolated maxima if locations are not within [-1,1], calculated extrema are ignored */
687 1721 : if ( posi >= 1.0f || posi <= -1.0f )
688 : {
689 1 : posi = y3 > y1 ? 1.0f : -1.0f;
690 : }
691 : else
692 : {
693 1720 : if ( y1 >= y2i )
694 : {
695 1 : posi = ( y1 > y3 ) ? -1.0f : 1.0f;
696 : }
697 1719 : else if ( y3 >= y2i )
698 : {
699 0 : posi = 1.0f;
700 : }
701 : }
702 :
703 1721 : return posi + 1.0f;
704 : }
705 :
706 :
707 : /*-------------------------------------------------------------------*
708 : * spec_ana()
709 : *
710 : * Spectral analysis
711 : *-------------------------------------------------------------------*/
712 :
713 188 : static void spec_ana(
714 : const float *prevsynth, /* i : Input signal */
715 : int16_t *plocs, /* o : The indicies of the identified peaks */
716 : float *plocsi, /* o : Interpolated positions of the identified peaks */
717 : int16_t *num_plocs, /* o : Number of identified peaks */
718 : float *X_sav, /* o : Stored fft spectrum */
719 : const int16_t output_frame, /* i : Frame length */
720 : const int16_t bwidth, /* i : Encoded bandwidth */
721 : const int16_t element_mode, /* i : IVAS element mode */
722 : float *noise_fac, /* o : for few peaks zeroing valleys decision making */
723 : const float pcorr )
724 : {
725 188 : int16_t i, Lprot, LprotLog2 = 0, hamm_len2 = 0, Lprot2_1, m;
726 : float *pPlocsi;
727 : int16_t *pPlocs;
728 : int16_t currPlocs, endPlocs, Lprot2p1, nJacob;
729 : int16_t n, k;
730 : int16_t st_point;
731 : int16_t end_point;
732 :
733 : float sig, noise, nsr;
734 : float window_corr_step, window_corr;
735 188 : const float *w_hamm = NULL;
736 : float xfp[L_PROT48k];
737 : float Xmax, Xmin, sel;
738 : int16_t stop_band_start;
739 : int16_t stop_band_length;
740 :
741 188 : Lprot = 2 * output_frame * L_PROT32k / 1280;
742 188 : Lprot2_1 = Lprot / 2 + 1;
743 :
744 188 : if ( output_frame == L_FRAME48k )
745 : {
746 136 : w_hamm = w_hamm_sana48k_2;
747 136 : hamm_len2 = L_PROT_HAMM_LEN2_48k;
748 : }
749 52 : else if ( output_frame == L_FRAME32k )
750 : {
751 52 : w_hamm = w_hamm_sana32k_2;
752 52 : hamm_len2 = L_PROT_HAMM_LEN2_32k;
753 52 : LprotLog2 = 10;
754 : }
755 : else
756 : {
757 0 : w_hamm = w_hamm_sana16k_2;
758 0 : hamm_len2 = L_PROT_HAMM_LEN2_16k;
759 0 : LprotLog2 = 9;
760 : }
761 :
762 : /* Apply hamming-rect window */
763 188 : mvr2r( prevsynth + hamm_len2, xfp + hamm_len2, Lprot - 2 * hamm_len2 );
764 188 : if ( element_mode == EVS_MONO )
765 : {
766 0 : v_mult( prevsynth, w_hamm, xfp, hamm_len2 );
767 0 : mult_rev2( prevsynth + Lprot - hamm_len2, w_hamm, xfp + Lprot - hamm_len2, hamm_len2 );
768 : }
769 : else
770 : {
771 188 : window_corr = w_hamm[0];
772 188 : window_corr_step = w_hamm[0] / hamm_len2;
773 49340 : for ( i = 0; i < hamm_len2; i++ )
774 : {
775 49152 : xfp[i] = prevsynth[i] * ( w_hamm[i] - window_corr );
776 49152 : xfp[Lprot - i - 1] = prevsynth[Lprot - i - 1] * ( w_hamm[i] - window_corr );
777 49152 : window_corr -= window_corr_step;
778 : }
779 : }
780 :
781 : /* Spectrum */
782 188 : if ( output_frame == L_FRAME48k )
783 : {
784 136 : fft3( xfp, xfp, Lprot );
785 : }
786 : else
787 : {
788 52 : fft_rel( xfp, Lprot, LprotLog2 );
789 : }
790 :
791 : /* Apply zeroing of non-coded FFT spectrum */
792 188 : if ( output_frame > inner_frame_tbl[bwidth] )
793 : {
794 19 : stop_band_start = 128 << bwidth;
795 19 : stop_band_length = Lprot - ( stop_band_start << 1 );
796 19 : stop_band_start = stop_band_start + 1;
797 19 : set_f( xfp + stop_band_start, 0, stop_band_length );
798 : }
799 :
800 188 : mvr2r( xfp, X_sav, Lprot );
801 :
802 : /* Magnitude representation */
803 188 : fft_spec2( xfp, Lprot );
804 :
805 131448 : for ( i = 0; i < Lprot2_1; i++ )
806 : {
807 131260 : xfp[i] = (float) sqrt( (double) xfp[i] );
808 : }
809 :
810 : /* Find maxima */
811 188 : maximum( xfp, Lprot2_1, &Xmax );
812 188 : minimum( xfp, Lprot2_1, &Xmin );
813 188 : if ( element_mode == EVS_MONO )
814 : {
815 0 : sel = ( Xmax - Xmin ) * ( 1.0f - PFIND_SENS );
816 : }
817 : else
818 : {
819 188 : sel = ( Xmax - Xmin ) * ( 1.0f - ST_PFIND_SENS );
820 : }
821 :
822 188 : peakfinder( xfp, Lprot2_1, plocs, num_plocs, sel, TRUE ); /* NB peak at xfp[0] and xfp Lprot2_1-1 may occur */
823 :
824 : /* Currently not the pitch correlation but some LF correlation */
825 188 : if ( element_mode != EVS_MONO && *num_plocs > 50 && pcorr < 0.6f )
826 : {
827 7 : *num_plocs = 0;
828 : }
829 :
830 188 : if ( element_mode == EVS_MONO )
831 : {
832 : /* Refine peaks */
833 0 : for ( m = 0; m < *num_plocs; m++ )
834 : {
835 0 : if ( plocs[m] == 0 )
836 : {
837 0 : plocsi[m] = plocs[m] + imax_pos( &xfp[plocs[m]] );
838 : }
839 0 : else if ( plocs[m] == Lprot / 2 )
840 : {
841 0 : plocsi[m] = plocs[m] - 2 + imax_pos( &xfp[plocs[m] - 2] );
842 : }
843 : else
844 : {
845 0 : plocsi[m] = plocs[m] - 1 + imax_pos( &xfp[plocs[m] - 1] );
846 : }
847 : }
848 : }
849 : else
850 : {
851 :
852 188 : Lprot2p1 = Lprot / 2 + 1;
853 :
854 : /* Refine peaks */
855 188 : pPlocsi = plocsi;
856 188 : pPlocs = plocs;
857 188 : n = *num_plocs; /* number of peaks to process */
858 :
859 : /* Special case-- The very 1st peak if it is at 0 index position (DC) */
860 : /* With DELTA_CORR_F0_INT == 2 one needs to handle both *pPlocs==0 and *pPlocs==1 */
861 188 : if ( n > 0 && *pPlocs == 0 ) /* Very 1st peak position possible to have a peak at 0/DC index position. */
862 : {
863 25 : *pPlocsi++ = *pPlocs + imax_pos( &xfp[*pPlocs] );
864 25 : pPlocs++;
865 25 : n = n - 1;
866 : }
867 :
868 188 : if ( n > 0 && *pPlocs == 1 ) /* Also 2nd peak position uses DC which makes jacobsen unsuitable. */
869 : {
870 24 : *pPlocsi++ = *pPlocs - 1 + imax_pos( &xfp[*pPlocs - 1] );
871 24 : currPlocs = *pPlocs++;
872 24 : n = n - 1;
873 : }
874 :
875 : /* All remaining peaks except the very last two possible integer positions */
876 188 : currPlocs = *pPlocs++;
877 188 : endPlocs = Lprot2p1 - DELTA_CORR_F0_INT; /* last *pPlocs position for Jacobsen */
878 :
879 : /* precompute number of turns based on endpoint integer location and make into a proper for loop */
880 188 : if ( n > 0 )
881 : {
882 180 : nJacob = n;
883 180 : if ( sub( endPlocs, plocs[sub( *num_plocs, 1 )] ) <= 0 )
884 : {
885 0 : nJacob = sub( nJacob, 1 );
886 : }
887 :
888 2880 : for ( k = 0; k < nJacob; k++ )
889 : {
890 2700 : *pPlocsi++ = currPlocs + imax2_jacobsen_mag( &( X_sav[currPlocs - 1] ), &( X_sav[Lprot - 1 - currPlocs] ) );
891 2700 : currPlocs = *pPlocs++;
892 : }
893 180 : n = n - nJacob;
894 : }
895 :
896 : /* At this point there should at most two plocs left to process */
897 : /* the position before fs/2 and fs/2 both use the same magnitude points */
898 188 : if ( n > 0 )
899 : {
900 : /* [ . . . . . . . ] Lprot/2+1 positions */
901 : /* | | | */
902 : /* 0 (Lprot/2-2) (Lprot/2) */
903 :
904 0 : if ( currPlocs == ( Lprot2p1 - DELTA_CORR_F0_INT ) ) /* Also 2nd last peak position uses fs/2 which makes jacobsen less suitable. */
905 : {
906 0 : *pPlocsi++ = currPlocs - 1 + imax_pos( &xfp[currPlocs - 1] );
907 0 : currPlocs = *pPlocs++;
908 0 : n = n - 1;
909 : }
910 :
911 : /* Here the only remaining point would be a fs/2 plocs */
912 : /* pXfp = xfp + sub(Lprot2,1); already set just a reminder where it
913 : * whould point */
914 0 : if ( n > 0 ) /* fs/2 which makes special case . */
915 : {
916 0 : *pPlocsi++ = currPlocs - 2 + imax_pos( &xfp[currPlocs - 2] );
917 0 : currPlocs = *pPlocs++;
918 0 : n = n - 1;
919 : }
920 : }
921 :
922 : /* For few peaks decide noise floor attenuation */
923 188 : if ( *num_plocs < 3 && *num_plocs > 0 )
924 : {
925 5 : sig = sum_f( xfp, Lprot2_1 ) + EPSILON;
926 :
927 : /*excluding peaks and neighboring bins*/
928 12 : for ( i = 0; i < *num_plocs; i++ )
929 : {
930 7 : st_point = max( 0, plocs[i] - DELTA_CORR );
931 7 : end_point = min( Lprot2_1 - 1, plocs[i] + DELTA_CORR );
932 7 : set_f( &xfp[st_point], 0.0f, end_point - st_point + 1 );
933 : }
934 5 : noise = sum_f( xfp, Lprot2_1 ) + EPSILON;
935 5 : nsr = noise / sig;
936 :
937 5 : if ( nsr < 0.03f )
938 : {
939 0 : *noise_fac = 0.5f;
940 : }
941 : else
942 : {
943 5 : *noise_fac = 1.0f;
944 : }
945 : }
946 : }
947 :
948 188 : return;
949 : }
950 :
951 : /*-------------------------------------------------------------------*
952 : * subst_spec()
953 : *
954 : * Substitution spectrum calculation
955 : *-------------------------------------------------------------------*/
956 :
957 276 : static void subst_spec(
958 : const int16_t *plocs, /* i : The indicies of the identified peaks */
959 : const float *plocsi, /* i : Interpolated positions of the identified peaks */
960 : int16_t *num_plocs, /* i/o: Number of identified peaks */
961 : const int16_t time_offs, /* i : Time offset */
962 : float *X, /* i/o: FFT spectrum */
963 : const float *mag_chg, /* i : Magnitude modification */
964 : const float ph_dith, /* i : Phase dither */
965 : const int16_t *is_trans, /* i : Transient flags */
966 : const int16_t output_frame, /* i : Frame length */
967 : int16_t *seed, /* i/o: Random seed */
968 : const float *alpha, /* i : Magnitude modification factors for fade to average */
969 : const float *beta, /* i : Magnitude modification factors for fade to average */
970 : float beta_mute, /* i : Factor for long-term mute */
971 : const float Xavg[LGW_MAX], /* i : Frequency group averages to fade to */
972 : const int16_t element_mode, /* i : IVAS element mode */
973 : const int16_t ph_ecu_lookahead, /* i : Phase ECU lookahead */
974 : const float noise_fac /* i : noise factor */
975 :
976 : )
977 : {
978 : const float *sincos;
979 : int16_t Xph_short;
980 : float corr_phase[MAX_PLOCS], Xph;
981 : float Lprot_1, cos_F, sin_F, tmp;
982 : int16_t Lprot, Lecu, m, i, e, im_ind, delta_corr_up, delta_corr_dn, delta_tmp;
983 : float mag_chg_local; /* for peak attenuation in burst */
984 : int16_t k;
985 : float one_peak_flag_mask;
986 : float alpha_local;
987 : float beta_local;
988 :
989 276 : sincos = sincos_t_ext + 128;
990 276 : Lprot = (int16_t) ( L_PROT32k * output_frame / 640 );
991 276 : Lprot_1 = 1.0f / Lprot;
992 276 : Lecu = output_frame * 2;
993 :
994 : /* Correction phase of the identified peaks */
995 276 : if ( is_trans[0] || is_trans[1] )
996 : {
997 13 : *num_plocs = 0;
998 : }
999 : else
1000 : {
1001 263 : tmp = PI2 * ( Lecu - ( Lecu - Lprot ) / 2 + NS2SA( output_frame * FRAMES_PER_SEC, PH_ECU_ALDO_OLP2_NS ) - ph_ecu_lookahead - output_frame / 2 + time_offs ) * Lprot_1;
1002 7887 : for ( m = 0; m < *num_plocs; m++ )
1003 : {
1004 7624 : corr_phase[m] = plocsi[m] * tmp;
1005 : }
1006 : }
1007 276 : one_peak_flag_mask = 1; /* all ones mask -> keep */
1008 276 : if ( element_mode != EVS_MONO )
1009 : {
1010 276 : if ( ( *num_plocs > 0 ) && sub( *num_plocs, 3 ) < 0 )
1011 : {
1012 5 : one_peak_flag_mask = noise_fac; /* all zeroes mask -> zero */
1013 : }
1014 276 : if ( *num_plocs == 0 )
1015 : {
1016 20 : X[0] = 0; /* reset DC if there are no peaks */
1017 20 : X[shr( Lprot, 1 )] = 0; /* also reset fs/2 if there are no peaks */
1018 : }
1019 : }
1020 :
1021 276 : i = 1;
1022 276 : k = 0;
1023 276 : im_ind = Lprot - 1;
1024 7900 : for ( m = 0; m < *num_plocs; m++ )
1025 : {
1026 7624 : delta_corr_dn = DELTA_CORR;
1027 7624 : delta_corr_up = DELTA_CORR;
1028 :
1029 7624 : if ( m > 0 )
1030 : {
1031 7368 : delta_tmp = ( plocs[m] - plocs[m - 1] - 1 ) / 2;
1032 7368 : if ( delta_tmp < DELTA_CORR )
1033 : {
1034 6221 : delta_corr_dn = delta_tmp;
1035 : }
1036 : }
1037 :
1038 7624 : if ( m < *num_plocs - 1 )
1039 : {
1040 7368 : delta_tmp = ( plocs[m + 1] - plocs[m] - 1 ) / 2;
1041 7368 : if ( delta_tmp < DELTA_CORR )
1042 : {
1043 6221 : delta_corr_up = delta_tmp;
1044 : }
1045 : }
1046 :
1047 : /* Input Xph */
1048 26474 : while ( i < plocs[m] - delta_corr_dn )
1049 : {
1050 18850 : *seed = own_random( seed );
1051 :
1052 18850 : if ( *seed & 0x40 )
1053 : {
1054 9393 : sin_F = sincos[*seed >> 8];
1055 : }
1056 : else
1057 : {
1058 9457 : sin_F = -sincos[*seed >> 8];
1059 : }
1060 :
1061 18850 : if ( *seed & 0x80 )
1062 : {
1063 9551 : cos_F = sincos[-( *seed >> 8 )];
1064 : }
1065 : else
1066 : {
1067 9299 : cos_F = -sincos[-( *seed >> 8 )];
1068 : }
1069 :
1070 18850 : if ( element_mode == EVS_MONO )
1071 : {
1072 0 : tmp = ( X[i] * cos_F - X[im_ind] * sin_F );
1073 0 : X[im_ind] = ( X[i] * sin_F + X[im_ind] * cos_F );
1074 : }
1075 : else
1076 : {
1077 18850 : tmp = one_peak_flag_mask * ( X[i] * cos_F - X[im_ind] * sin_F );
1078 18850 : X[im_ind] = one_peak_flag_mask * ( X[i] * sin_F + X[im_ind] * cos_F );
1079 : }
1080 :
1081 18850 : if ( alpha[k] < 1.0f )
1082 : {
1083 6410 : *seed = rand_phase( *seed, &sin_F, &cos_F );
1084 6410 : X[i] = alpha[k] * tmp + beta[k] * Xavg[k] * cos_F;
1085 6410 : X[im_ind] = alpha[k] * X[im_ind] + beta[k] * Xavg[k] * sin_F;
1086 : }
1087 : else
1088 : {
1089 12440 : X[i] = mag_chg[k] * tmp;
1090 12440 : X[im_ind] *= mag_chg[k];
1091 : }
1092 18850 : i++;
1093 18850 : im_ind--;
1094 18850 : if ( i >= ivas_gwlpr[k + 1] )
1095 : {
1096 185 : k++;
1097 : }
1098 : }
1099 :
1100 7624 : e = plocs[m] + delta_corr_up;
1101 7624 : if ( e > Lprot / 2 - 1 )
1102 : {
1103 0 : e = Lprot / 2 - 1;
1104 : }
1105 :
1106 7624 : Xph = corr_phase[m];
1107 : /* extract fractional phase integer index in the range [0...1023] */
1108 7624 : Xph_short = (int16_t) ( 0x000003ff & (int32_t) ( ( Xph * 512 ) / EVS_PI ) );
1109 :
1110 :
1111 7624 : if ( Xph_short >= 512 )
1112 : {
1113 3834 : sin_F = -sincos_t_ext[Xph_short - 512];
1114 3834 : if ( Xph_short < 768 )
1115 : {
1116 1865 : cos_F = -sincos_t_ext[Xph_short - 512 + 256];
1117 : }
1118 : else
1119 : {
1120 1969 : cos_F = sincos_t_ext[-Xph_short + 1024 + 256];
1121 : }
1122 : }
1123 : else
1124 : {
1125 3790 : sin_F = sincos_t_ext[Xph_short];
1126 3790 : if ( Xph_short < 256 )
1127 : {
1128 2022 : cos_F = sincos_t_ext[Xph_short + 256];
1129 : }
1130 : else
1131 : {
1132 1768 : cos_F = -sincos_t_ext[-Xph_short + 256 + 512];
1133 : }
1134 : }
1135 :
1136 54107 : while ( i <= e )
1137 : {
1138 46483 : mag_chg_local = mag_chg[k];
1139 :
1140 46483 : if ( ph_dith != 0.0f )
1141 : {
1142 : /* Call phase randomization only when needed */
1143 24732 : Xph = corr_phase[m];
1144 24732 : *seed = own_random( seed );
1145 24732 : Xph += *seed * ph_dith * PHASE_DITH_SCALE; /* where ph_dith is 0..2PI, or -2PI (in transient), bin phase scaling factor from trans_ana */
1146 :
1147 24732 : if ( ph_dith > 0.0f )
1148 : {
1149 : /* up to 6 dB additional att of peaks in non_transient longer bursts, (when peak phase is randomized) */
1150 : /* 0.5~= sqrt((float)pow(10.0,-6/10.0)); ph_dith= 0..2pi,--> scale=1.0 ...5 */
1151 24732 : mag_chg_local *= 0.5f + ( 1.0f - ( 1.0f / PHASE_DITH ) * ph_dith ) * 0.5f;
1152 : }
1153 :
1154 24732 : Xph_short = (int16_t) ( ( (int32_t) ( ( Xph * 512 ) / EVS_PI ) ) & 0x000003ff );
1155 :
1156 24732 : if ( Xph_short >= 512 )
1157 : {
1158 12354 : sin_F = -sincos_t_ext[Xph_short - 512];
1159 12354 : if ( Xph_short < 768 )
1160 : {
1161 6074 : cos_F = -sincos_t_ext[Xph_short - 512 + 256];
1162 : }
1163 : else
1164 : {
1165 6280 : cos_F = sincos_t_ext[-Xph_short + 1024 + 256];
1166 : }
1167 : }
1168 : else
1169 : {
1170 12378 : sin_F = sincos_t_ext[Xph_short];
1171 12378 : if ( Xph_short < 256 )
1172 : {
1173 6167 : cos_F = sincos_t_ext[Xph_short + 256];
1174 : }
1175 : else
1176 : {
1177 6211 : cos_F = -sincos_t_ext[-Xph_short + 256 + 512];
1178 : }
1179 : }
1180 : }
1181 :
1182 46483 : tmp = ( X[i] * cos_F - X[im_ind] * sin_F );
1183 46483 : X[im_ind] = ( X[i] * sin_F + X[im_ind] * cos_F );
1184 46483 : if ( alpha[k] < 1.0f )
1185 : {
1186 25446 : alpha_local = mag_chg_local;
1187 25446 : beta_local = (float) ( beta_mute * sqrt( 1.0f - SQR( alpha_local ) ) );
1188 25446 : if ( k >= LGW32k - 1 )
1189 : {
1190 9879 : beta_local *= 0.1f;
1191 : }
1192 15567 : else if ( k >= LGW16k - 1 )
1193 : {
1194 7733 : beta_local *= 0.5f;
1195 : }
1196 :
1197 25446 : *seed = rand_phase( *seed, &sin_F, &cos_F );
1198 25446 : X[i] = alpha_local * tmp + beta_local * Xavg[k] * cos_F;
1199 25446 : X[im_ind] = alpha_local * X[im_ind] + beta_local * Xavg[k] * sin_F;
1200 : }
1201 : else
1202 : {
1203 21037 : X[i] = mag_chg_local * tmp;
1204 21037 : X[im_ind] *= mag_chg_local;
1205 : }
1206 :
1207 46483 : i++;
1208 46483 : im_ind--;
1209 46483 : if ( i >= ivas_gwlpr[k + 1] )
1210 : {
1211 1100 : k++;
1212 : }
1213 : }
1214 : }
1215 :
1216 132555 : while ( i < Lprot / 2 )
1217 : {
1218 132279 : *seed = own_random( seed );
1219 :
1220 132279 : if ( *seed & 0x40 )
1221 : {
1222 65829 : sin_F = sincos[*seed >> 8];
1223 : }
1224 : else
1225 : {
1226 66450 : sin_F = -sincos[*seed >> 8];
1227 : }
1228 :
1229 132279 : if ( *seed & 0x80 )
1230 : {
1231 66103 : cos_F = sincos[-( *seed >> 8 )];
1232 : }
1233 : else
1234 : {
1235 66176 : cos_F = -sincos[-( *seed >> 8 )];
1236 : }
1237 :
1238 :
1239 132279 : if ( element_mode == EVS_MONO )
1240 : {
1241 0 : tmp = ( X[i] * cos_F - X[im_ind] * sin_F );
1242 0 : X[im_ind] = ( X[i] * sin_F + X[im_ind] * cos_F );
1243 : }
1244 : else
1245 : {
1246 132279 : tmp = one_peak_flag_mask * ( X[i] * cos_F - X[im_ind] * sin_F );
1247 132279 : X[im_ind] = one_peak_flag_mask * ( X[i] * sin_F + X[im_ind] * cos_F );
1248 : }
1249 :
1250 132279 : if ( alpha[k] < 1.0f )
1251 : {
1252 25194 : *seed = rand_phase( *seed, &sin_F, &cos_F );
1253 25194 : X[i] = alpha[k] * tmp + beta[k] * Xavg[k] * cos_F;
1254 25194 : X[im_ind] = alpha[k] * X[im_ind] + beta[k] * Xavg[k] * sin_F;
1255 25194 : im_ind--;
1256 : }
1257 : else
1258 : {
1259 107085 : X[i] = mag_chg[k] * tmp;
1260 107085 : X[im_ind--] *= mag_chg[k];
1261 : }
1262 132279 : i++;
1263 :
1264 132279 : if ( i >= ivas_gwlpr[k + 1] )
1265 : {
1266 868 : k++;
1267 : }
1268 : }
1269 :
1270 276 : return;
1271 : }
1272 :
1273 : /*--------------------------------------------------------------------------
1274 : * rec_wtda()
1275 : *
1276 : * Windowing and TDA of reconstructed frame
1277 : *--------------------------------------------------------------------------*/
1278 :
1279 276 : static void rec_wtda(
1280 : float *X, /* i/o: ECU frame / unwindowed ECU frame */
1281 : float *ecu_rec, /* o : Reconstructed frame in tda domain */
1282 : const int16_t output_frame, /* i : Frame length */
1283 : const int16_t Lprot, /* i : Prototype frame length */
1284 : const float old_dec[270], /* i : end of last decoded for OLA before tda and itda */
1285 : const int16_t element_mode, /* i : IVAS element mode */
1286 : const int16_t *num_p, /* i : Number of peaks */
1287 : const int16_t *plocs /* i : Peak locations */
1288 : )
1289 : {
1290 : int16_t timesh;
1291 : int16_t xf_len;
1292 : int16_t i;
1293 : float *p_ecu;
1294 : float g;
1295 : float tbl_delta;
1296 : float xsubst_[2 * L_FRAME48k];
1297 : const float *w_hamm;
1298 : float *pX_start, *pX_end;
1299 : float tmp;
1300 : int16_t hamm_len2;
1301 : float *pNew;
1302 : const float *pOldW, *pNewW;
1303 : float xfwin[NS2SA( L_FRAME48k * FRAMES_PER_SEC, N_ZERO_MDCT_NS - ( 2 * FRAME_SIZE_NS - L_PROT_NS ) / 2 )];
1304 : const float *pOld;
1305 : int16_t copy_len;
1306 : int16_t ola_len;
1307 :
1308 276 : copy_len = NS2SA( output_frame * FRAMES_PER_SEC, ( 2 * FRAME_SIZE_NS - L_PROT_NS ) / 2 ); /* prototype fill on each side of xsubst to fill MDCT Frame */
1309 276 : ola_len = NS2SA( output_frame * FRAMES_PER_SEC, N_ZERO_MDCT_NS - ( 2 * FRAME_SIZE_NS - L_PROT_NS ) / 2 ); /* remaining lengt of LA_ZEROS to overlap add decoded with xsubst */
1310 :
1311 276 : if ( output_frame == L_FRAME48k )
1312 : {
1313 221 : w_hamm = w_hamm_sana48k_2;
1314 221 : hamm_len2 = L_PROT_HAMM_LEN2_48k;
1315 : }
1316 55 : else if ( output_frame == L_FRAME32k )
1317 : {
1318 55 : w_hamm = w_hamm_sana32k_2;
1319 55 : hamm_len2 = L_PROT_HAMM_LEN2_32k;
1320 : }
1321 : else
1322 : {
1323 0 : w_hamm = w_hamm_sana16k_2;
1324 0 : hamm_len2 = L_PROT_HAMM_LEN2_16k;
1325 : }
1326 :
1327 276 : if ( element_mode != EVS_MONO && *num_p > 0 && plocs[0] > 3 )
1328 : {
1329 : /* Perform inverse windowing of hammrect */
1330 105 : pX_start = X;
1331 105 : pX_end = X + Lprot - 1;
1332 27753 : for ( i = 0; i < hamm_len2; i++ )
1333 : {
1334 27648 : tmp = 1.0f / *w_hamm;
1335 27648 : *pX_start *= tmp;
1336 27648 : *pX_end *= tmp;
1337 27648 : pX_start++;
1338 27648 : pX_end--;
1339 27648 : w_hamm++;
1340 : }
1341 : }
1342 :
1343 : /* extract reconstructed frame with aldo window */
1344 276 : timesh = NS2SA( output_frame * FRAMES_PER_SEC, N_ZERO_MDCT_NS ) - ( 2 * output_frame - Lprot ) / 2;
1345 :
1346 276 : set_f( xsubst_, 0.0f, 2 * output_frame - Lprot + timesh );
1347 276 : mvr2r( X, xsubst_ + 2 * output_frame - Lprot + timesh, Lprot - timesh );
1348 :
1349 : /* Copy and OLA look ahead zero part of MDCT window from decoded signal */
1350 276 : if ( element_mode != EVS_MONO )
1351 : {
1352 276 : mvr2r( old_dec, xsubst_ + NS2SA( output_frame * FRAMES_PER_SEC, N_ZERO_MDCT_NS ), copy_len ); /* also need to scale to Q0 ?? */
1353 276 : pOld = old_dec + copy_len;
1354 276 : pNew = xsubst_ + copy_len + NS2SA( output_frame * FRAMES_PER_SEC, N_ZERO_MDCT_NS );
1355 276 : sinq( EVS_PI / ( ola_len * 2 ), 0.0f, ola_len, xfwin );
1356 276 : v_mult( xfwin, xfwin, xfwin, ola_len ); /* xfwin = sin^2 of 0..pi/4 */
1357 276 : pOldW = xfwin + ola_len - 1;
1358 276 : pNewW = xfwin;
1359 20374 : for ( i = 0; i < ola_len; i++ )
1360 : {
1361 20098 : *pNew = *pOld * *pOldW + *pNew * *pNewW;
1362 20098 : pOld += 1;
1363 20098 : pNew += 1;
1364 20098 : pOldW -= 1;
1365 20098 : pNewW += 1;
1366 : }
1367 : }
1368 : else
1369 : {
1370 : /* Smoothen onset of ECU frame */
1371 0 : xf_len = (int16_t) ( (float) output_frame * N_ZERO_MDCT_NS / FRAME_SIZE_NS ) - ( output_frame - Lprot / 2 );
1372 0 : p_ecu = xsubst_ + 2 * output_frame - Lprot + timesh;
1373 0 : tbl_delta = 64.f / xf_len; /* 64 samples = 1/4 cycle in sincos_t */
1374 0 : for ( i = 0; i < xf_len; i++, p_ecu++ )
1375 : {
1376 0 : g = sincos_t[( (int16_t) ( i * tbl_delta ) )];
1377 0 : g *= g;
1378 0 : *p_ecu = g * ( *p_ecu );
1379 : }
1380 : }
1381 :
1382 : /* Apply TDA and windowing to ECU frame */
1383 276 : wtda( xsubst_ + output_frame, ecu_rec, NULL, ALDO_WINDOW, ALDO_WINDOW, output_frame );
1384 :
1385 276 : return;
1386 : }
1387 :
1388 :
1389 : /*--------------------------------------------------------------------------
1390 : * rec_frame()
1391 : *
1392 : * Frame reconstruction
1393 : *--------------------------------------------------------------------------*/
1394 :
1395 276 : static void rec_frame(
1396 : float *X, /* i/o: FFT spectrum / IFFT of spectrum */
1397 : float *ecu_rec, /* o : Reconstructed frame in tda domain */
1398 : const int16_t output_frame, /* i : Frame length */
1399 : const float *old_dec, /* i : end of last decoded for OLA before tda and itda */
1400 : const int16_t element_mode, /* i : IVAS element mode */
1401 : const int16_t *num_p, /* i : Number of peaks */
1402 : const int16_t *plocs /* i : Peak locations */
1403 : )
1404 : {
1405 276 : int16_t Lprot, LprotLog2 = 0;
1406 :
1407 276 : Lprot = 2 * output_frame * L_PROT32k / 1280;
1408 :
1409 276 : if ( output_frame == L_FRAME48k )
1410 : {
1411 221 : LprotLog2 = 9;
1412 : }
1413 55 : else if ( output_frame == L_FRAME32k )
1414 : {
1415 55 : LprotLog2 = 10;
1416 : }
1417 : else
1418 : {
1419 0 : LprotLog2 = 9;
1420 : }
1421 :
1422 : /* extend spectrum and IDFT */
1423 276 : if ( output_frame == L_FRAME48k )
1424 : {
1425 221 : ifft3( X, X, Lprot );
1426 : }
1427 : else
1428 : {
1429 55 : ifft_rel( X, Lprot, LprotLog2 );
1430 : }
1431 :
1432 276 : rec_wtda( X, ecu_rec, output_frame, Lprot, old_dec, element_mode, num_p, plocs );
1433 :
1434 276 : return;
1435 : }
1436 :
1437 :
1438 : /*--------------------------------------------------------------------------
1439 : * fir_dwn()
1440 : *
1441 : * FIR downsampling filter
1442 : *--------------------------------------------------------------------------*/
1443 :
1444 277 : static void fir_dwn(
1445 : const float x[], /* i : input vector */
1446 : const float h[], /* i : impulse response of the FIR filter */
1447 : float y[], /* o : output vector (result of filtering) */
1448 : const int16_t L, /* i : input vector size */
1449 : const int16_t K, /* i : order of the FIR filter (K+1 coefs.) */
1450 : const int16_t decimation /* i : decimation */
1451 : )
1452 : {
1453 : float s;
1454 : int16_t i, j, k, mmax;
1455 :
1456 277 : k = 0;
1457 :
1458 : /* do the filtering */
1459 87532 : for ( i = K / 2; i < L; i += decimation )
1460 : {
1461 87255 : s = x[i] * h[0];
1462 :
1463 87255 : if ( i < K )
1464 : {
1465 1385 : mmax = i;
1466 : }
1467 : else
1468 : {
1469 85870 : mmax = K;
1470 : }
1471 :
1472 4946505 : for ( j = 1; j <= mmax; j++ )
1473 : {
1474 4859250 : s += h[j] * x[i - j];
1475 : }
1476 :
1477 87255 : y[k] = s;
1478 87255 : k++;
1479 : }
1480 :
1481 1662 : for ( ; i < L + K / 2; i += decimation )
1482 : {
1483 1385 : s = 0;
1484 :
1485 63385 : for ( j = i - L + 1; j <= K; j++ )
1486 : {
1487 62000 : s += h[j] * x[i - j];
1488 : }
1489 :
1490 1385 : y[k] = s;
1491 1385 : k++;
1492 : }
1493 :
1494 277 : return;
1495 : }
1496 :
1497 :
1498 : /*--------------------------------------------------------------------------
1499 : * fec_ecu_pitch()
1500 : *
1501 : * Pitch/correlation analysis and adaptive analysis frame length calculation
1502 : *--------------------------------------------------------------------------*/
1503 :
1504 277 : static void fec_ecu_pitch(
1505 : const float *prevsynth, /* i : previous synthesis */
1506 : float *prevsynth_LP, /* o : down-sampled synthesis */
1507 : const int16_t L, /* i : Output frame length */
1508 : int16_t *N, /* o : Analysis frame length in 8kHz */
1509 : float *min_corr, /* o : correlation */
1510 : int16_t *decimatefator, /* o : Decimation factor */
1511 : const int16_t HqVoicing /* i : Hq voicing flag */
1512 : )
1513 : {
1514 : int16_t i, filt_size;
1515 : float accA, accB, accC, Ryy;
1516 : int16_t delay_ind, k, cb_start, cb_end, tmp_short, Lon20;
1517 : float Asr_LP[FEC_DCIM_FILT_SIZE_MAX + 1];
1518 :
1519 277 : switch ( L )
1520 : {
1521 221 : case L_FRAME48k:
1522 221 : *decimatefator = 6;
1523 221 : filt_size = 60;
1524 221 : mvr2r( Asr_LP48, Asr_LP, filt_size + 1 );
1525 221 : break;
1526 56 : case L_FRAME32k:
1527 56 : *decimatefator = 4;
1528 56 : filt_size = 40;
1529 56 : mvr2r( Asr_LP32, Asr_LP, filt_size + 1 );
1530 56 : break;
1531 0 : case L_FRAME16k:
1532 0 : *decimatefator = 2;
1533 0 : filt_size = 20;
1534 0 : mvr2r( Asr_LP16, Asr_LP, filt_size + 1 );
1535 0 : break;
1536 0 : default:
1537 0 : *decimatefator = 2;
1538 0 : filt_size = 40;
1539 0 : mvr2r( Asr_LP16, Asr_LP, filt_size + 1 );
1540 0 : break;
1541 : }
1542 :
1543 : /* We need to inverse the ALDO window */
1544 :
1545 : /* Resampling to work at 8Khz */
1546 277 : fir_dwn( prevsynth, Asr_LP, prevsynth_LP, 2 * L, filt_size, *decimatefator ); /* resampling without delay */
1547 :
1548 277 : Lon20 = (int16_t) ( ( L / 20 ) / *decimatefator );
1549 :
1550 : /* Correlation analysis */
1551 277 : *min_corr = 0;
1552 277 : accC = 0;
1553 13573 : for ( k = 0; k < 6 * Lon20; k++ )
1554 : {
1555 13296 : accC += prevsynth_LP[34 * Lon20 + k] * prevsynth_LP[34 * Lon20 + k];
1556 : }
1557 :
1558 277 : if ( HqVoicing == 1 )
1559 : {
1560 0 : cb_start = 0;
1561 0 : cb_end = 33 * Lon20;
1562 : }
1563 : else
1564 : {
1565 277 : cb_start = 0;
1566 277 : cb_end = 28 * Lon20;
1567 : }
1568 :
1569 277 : tmp_short = 34 * Lon20;
1570 277 : accB = 0;
1571 277 : delay_ind = cb_start;
1572 47993 : for ( i = cb_start; i < cb_end; i++ ) /* cb_end = 35 let 6 ms min of loop size */
1573 : {
1574 47802 : accA = 0;
1575 47802 : if ( i == cb_start )
1576 : {
1577 277 : accB = 0;
1578 13573 : for ( k = 0; k < 6 * Lon20; k++ )
1579 : {
1580 13296 : accA += prevsynth_LP[i + k] * prevsynth_LP[tmp_short + k];
1581 13296 : accB += prevsynth_LP[i + k] * prevsynth_LP[i + k];
1582 : }
1583 : }
1584 : else
1585 : {
1586 47525 : accB = accB - prevsynth_LP[i - 1] * prevsynth_LP[i - 1] + prevsynth_LP[i + 6 * Lon20 - 1] * prevsynth_LP[i + 6 * Lon20 - 1];
1587 2328725 : for ( k = 0; k < 6 * Lon20; k++ )
1588 : {
1589 2281200 : accA += prevsynth_LP[i + k] * prevsynth_LP[tmp_short + k];
1590 : }
1591 : }
1592 :
1593 : /* tests to avoid division by zero */
1594 47802 : if ( accB <= 0 )
1595 : {
1596 0 : accB = 1.0f;
1597 : }
1598 :
1599 47802 : if ( accC <= 0 )
1600 : {
1601 0 : accC = 1.0f;
1602 : }
1603 :
1604 47802 : if ( accB * accC <= 0 )
1605 : {
1606 0 : Ryy = accA;
1607 : }
1608 : else
1609 : {
1610 47802 : Ryy = accA / (float) sqrt( ( accB * accC ) );
1611 : }
1612 :
1613 47802 : if ( Ryy > *min_corr )
1614 : {
1615 2316 : *min_corr = Ryy;
1616 2316 : delay_ind = i;
1617 : }
1618 :
1619 47802 : if ( HqVoicing == 0 && *min_corr > 0.95f )
1620 : {
1621 86 : break;
1622 : }
1623 : }
1624 :
1625 277 : *N = 40 * Lon20 - delay_ind - 6 * Lon20;
1626 :
1627 277 : return;
1628 : }
1629 :
1630 :
1631 : /*--------------------------------------------------------------------------
1632 : * fec_ecu_dft()
1633 : *
1634 : * DFT analysis on adaptive frame length. Analysis frame stretched to
1635 : * next power of 2 using linear interpolation.
1636 : *--------------------------------------------------------------------------*/
1637 :
1638 1 : static void fec_ecu_dft(
1639 : const float *prevsynth_LP, /* i : Downsampled past synthesis (2*160 samples) */
1640 : const int16_t N, /* i : Analysis frame length in 8 kHz (corr. max) */
1641 : float *Tfr, /* o : DFT coefficients, real part */
1642 : float *Tfi, /* o : DFT coefficients, imag part */
1643 : float *sum_Tf_abs, /* o : Sum of magnitude spectrum */
1644 : float *Tf_abs, /* o : Magnitude spectrum */
1645 : int16_t *Nfft, /* o : DFT analysis length 2^(nextpow2(N)) */
1646 : const int16_t element_mode /* i : IVAS element mode */
1647 : )
1648 : {
1649 : float target[2 * L_FRAME8k], tmp;
1650 : int16_t i, Lon20, tmp_short, N_LP, k;
1651 : int16_t alignment_point;
1652 :
1653 1 : Lon20 = (int16_t) 160 / 20;
1654 1 : if ( element_mode == EVS_MONO )
1655 : {
1656 0 : alignment_point = 2 * 160 - 3 * Lon20;
1657 : }
1658 : else
1659 : {
1660 1 : alignment_point = 2 * 160;
1661 : }
1662 :
1663 :
1664 55 : for ( i = 0; i < N; i++ )
1665 : {
1666 54 : target[i] = prevsynth_LP[alignment_point - N + i];
1667 : }
1668 :
1669 : /* DFT */
1670 1 : *sum_Tf_abs = 0;
1671 1 : *Nfft = (int16_t) pow( 2, (int16_t) ceil( log( N ) / log( 2 ) ) );
1672 1 : tmp = ( (float) N - 1.0f ) / ( (float) *Nfft - 1.0f );
1673 :
1674 1 : set_f( Tfr, 0.0f, *Nfft );
1675 1 : set_f( Tfi, 0.0f, *Nfft );
1676 1 : Tfr[0] = target[0];
1677 1 : Tfr[*Nfft - 1] = target[N - 1];
1678 63 : for ( i = 1; i < *Nfft - 1; i++ ) /* interpolation for FFT */
1679 : {
1680 62 : tmp_short = (int16_t) floor( i * tmp );
1681 62 : Tfr[i] = target[tmp_short] + ( (float) i * tmp - ( (float) tmp_short ) ) * ( target[tmp_short + 1] - target[tmp_short] );
1682 : }
1683 :
1684 1 : DoRTFTn( Tfr, Tfi, *Nfft );
1685 1 : N_LP = (int16_t) ceil( *Nfft / 2 );
1686 :
1687 33 : for ( k = 0; k < N_LP; k++ )
1688 : {
1689 32 : Tf_abs[k] = (float) sqrt( Tfr[k] * Tfr[k] + Tfi[k] * Tfi[k] );
1690 32 : *sum_Tf_abs += Tf_abs[k];
1691 : }
1692 :
1693 1 : return;
1694 : }
1695 :
1696 : /*--------------------------------------------------------------------------*
1697 : * singenerator()
1698 : *
1699 : * fast cosinus generator Amp*cos(2*pi*freq+phi)
1700 : *--------------------------------------------------------------------------*/
1701 :
1702 :
1703 8 : static void singenerator(
1704 : const int16_t L, /* i : size of output */
1705 : const float cosfreq, /* i : cosine of 1-sample dephasing at the given frequency */
1706 : const float sinfreq, /* i : sine of 1-sample dephasing at the given frequency */
1707 : const float a_re, /* i : real part of complex spectral coefficient at the given frequency */
1708 : const float a_im, /* i : imag part of complex spectral coefficient at the given frequency */
1709 : float xx[] /* o : output vector */
1710 : )
1711 : {
1712 :
1713 : float *ptr, L_C0, L_S0, L_C1, L_S1;
1714 : float C0, S0, C1, S1;
1715 : int16_t i;
1716 :
1717 8 : L_C0 = a_re;
1718 8 : S0 = a_im;
1719 :
1720 8 : ptr = xx;
1721 :
1722 8 : *ptr = *ptr + L_C0;
1723 8 : ptr++;
1724 :
1725 5120 : for ( i = 0; i < L / 2 - 1; i++ )
1726 : {
1727 5112 : C0 = L_C0;
1728 5112 : L_C1 = C0 * cosfreq;
1729 5112 : L_C1 = L_C1 - S0 * sinfreq;
1730 5112 : L_S1 = C0 * sinfreq;
1731 5112 : S1 = L_S1 + S0 * cosfreq;
1732 5112 : *ptr = *ptr + L_C1;
1733 5112 : ptr++;
1734 :
1735 5112 : C1 = L_C1;
1736 5112 : L_C0 = C1 * cosfreq;
1737 5112 : L_C0 = L_C0 - S1 * sinfreq;
1738 5112 : L_S0 = C1 * sinfreq;
1739 5112 : S0 = L_S0 + S1 * cosfreq;
1740 5112 : *ptr = *ptr + L_C0;
1741 5112 : ptr++;
1742 : }
1743 :
1744 8 : C0 = L_C0;
1745 8 : L_C1 = C0 * cosfreq;
1746 8 : L_C1 = L_C1 - S0 * sinfreq;
1747 8 : *ptr = *ptr + L_C1;
1748 8 : ptr++;
1749 :
1750 8 : return;
1751 : }
1752 :
1753 :
1754 : /*--------------------------------------------------------------------------
1755 : * sinusoidal_synthesis()
1756 : *
1757 : * ECU frame sinusoid generation
1758 : *--------------------------------------------------------------------------*/
1759 :
1760 1 : static void sinusoidal_synthesis(
1761 : const float *Tfr, /* i : DFT coefficients, real part */
1762 : const float *Tfi, /* i : DFT coefficients, imag part */
1763 : float *Tf_abs, /* i : Magnitude spectrum */
1764 : const int16_t N, /* i : Analysis frame length in 8 kHz (corr. max) */
1765 : const int16_t L, /* i : Output frame length */
1766 : const int16_t decimate_factor, /* i : Downsampling factor */
1767 : const int16_t Nfft, /* i : DFT analysis length 2^(nextpow2(N)) */
1768 : const float sum_Tf_abs, /* i : Sum of magnitude spectrum */
1769 : float *synthesis, /* o : ECU sinusoidal synthesis */
1770 : const int16_t HqVoicing /* i : HQ voicing flag */
1771 : )
1772 : {
1773 1 : int16_t i, k, nb_pulses, indmax = 0, nb_pulses_final;
1774 : int16_t pulses[FEC_MAX / 2];
1775 : float a_re[FEC_NB_PULSE_MAX], a_im[FEC_NB_PULSE_MAX], cosfreq, sinfreq;
1776 : float freq[FEC_NB_PULSE_MAX], tmp;
1777 : float mmax, cumsum;
1778 1 : int16_t Lon20 = 8;
1779 :
1780 : /* peak selection */
1781 : int16_t PL, cpt;
1782 : float old, new_s;
1783 : int16_t *p_pulses;
1784 : int16_t glued;
1785 :
1786 1 : p_pulses = pulses;
1787 1 : nb_pulses = 0;
1788 1 : new_s = Tf_abs[1];
1789 1 : glued = 1;
1790 1 : cpt = 0;
1791 1 : old = 0;
1792 :
1793 1 : PL = 0;
1794 1 : if ( N > Lon20 * 10 || HqVoicing )
1795 : {
1796 0 : PL = 1;
1797 : }
1798 26 : while ( cpt <= N / 2 - 1 - 2 )
1799 : {
1800 25 : if ( Tf_abs[cpt] > old && Tf_abs[cpt] > new_s )
1801 : {
1802 11 : glued = cpt;
1803 :
1804 22 : for ( i = glued; i < cpt + PL + 1; i++ )
1805 : {
1806 11 : *p_pulses++ = i;
1807 11 : nb_pulses++;
1808 : }
1809 11 : old = Tf_abs[cpt + PL];
1810 11 : new_s = Tf_abs[cpt + 2 + PL];
1811 11 : cpt = cpt + PL + 1;
1812 11 : glued = 1;
1813 : }
1814 : else
1815 : {
1816 14 : old = Tf_abs[cpt];
1817 14 : new_s = Tf_abs[cpt + 2];
1818 14 : cpt++;
1819 14 : glued = 0;
1820 : }
1821 : }
1822 :
1823 :
1824 1 : nb_pulses_final = 0;
1825 :
1826 : /* peak selection : keep the more energetics (max 20) */
1827 1 : tmp = 1.0f / (float) ( Nfft / 2 );
1828 1 : cumsum = 0;
1829 8 : for ( i = 0; i < min( FEC_NB_PULSE_MAX, nb_pulses ); i++ )
1830 : {
1831 8 : mmax = 0;
1832 96 : for ( k = 0; k < nb_pulses; k++ )
1833 : {
1834 88 : if ( Tf_abs[pulses[k]] > mmax )
1835 : {
1836 24 : mmax = Tf_abs[pulses[k]];
1837 24 : indmax = pulses[k];
1838 : }
1839 : }
1840 8 : cumsum += Tf_abs[indmax];
1841 :
1842 8 : if ( HqVoicing || cumsum < sum_Tf_abs * 0.7f )
1843 : {
1844 7 : a_re[i] = Tfr[indmax] * tmp;
1845 7 : a_im[i] = Tfi[indmax] * tmp;
1846 7 : freq[i] = (float) indmax * 2 / ( (float) N );
1847 7 : nb_pulses_final++;
1848 7 : Tf_abs[indmax] = -1;
1849 : }
1850 : else
1851 : {
1852 1 : a_re[i] = Tfr[indmax] * tmp;
1853 1 : a_im[i] = Tfi[indmax] * tmp;
1854 1 : freq[i] = (float) indmax * 2 / ( (float) N );
1855 1 : nb_pulses_final++;
1856 1 : break;
1857 : }
1858 : }
1859 :
1860 1 : nb_pulses = nb_pulses_final;
1861 :
1862 : /* sinusoidal synthesis */
1863 1 : set_f( synthesis, 0.0f, 40 * L / 20 );
1864 :
1865 1 : if ( HqVoicing )
1866 : {
1867 0 : k = 40 * L / 20;
1868 : }
1869 : else
1870 : {
1871 1 : k = 40 * L / 20;
1872 : }
1873 :
1874 9 : for ( i = 0; i < nb_pulses; i++ )
1875 : {
1876 8 : cosfreq = (float) cos( EVS_PI * freq[i] / (float) decimate_factor );
1877 8 : sinfreq = (float) sin( EVS_PI * freq[i] / (float) decimate_factor );
1878 8 : singenerator( k, cosfreq, sinfreq, a_re[i], a_im[i], synthesis );
1879 : }
1880 :
1881 1 : return;
1882 : }
1883 :
1884 : /*--------------------------------------------------------------------------
1885 : * fec_noise_filling()
1886 : *
1887 : * Adds noise component to ECU frame by
1888 : * - subtracting the sinusoidal synthesis from the analysis frame
1889 : * - copying noise segments at random indices and adding them together
1890 : * with an overlap-add operation
1891 : * Copying the beginning of the frame from the past synthesis and aligning
1892 : * it to be inserted into wtda
1893 : *--------------------------------------------------------------------------*/
1894 :
1895 1 : static void fec_noise_filling(
1896 : const float *prevsynth, /* i : Past synthesis buffer (length 2*L) */
1897 : float *synthesis, /* i/o: Sinusoidal ECU / Sinusoidal ECU + noise */
1898 : const int16_t L, /* i : Output frame length */
1899 : const int16_t N, /* i : Analysis frame length in 8 kHz (corr. max) */
1900 : const int16_t HqVoicing, /* i : HQ voicing flag */
1901 : float *gapsynth, /* o : Buffer for crossfade to good frame */
1902 : int16_t *ni_seed_forfec, /* i/o: Random seed for picking noise segments */
1903 : const int16_t element_mode, /* i : IVAS element mode */
1904 : const float *old_out )
1905 : {
1906 : float SS[L_FRAME48k / 2], tmp;
1907 : int16_t Rnd_N_noise;
1908 : int16_t k, kk, i;
1909 : int16_t N_noise;
1910 : float noisevect[34 * L_FRAME48k / 20];
1911 : const float *p_mdct_ola;
1912 : int16_t alignment_point;
1913 :
1914 1 : if ( element_mode == EVS_MONO )
1915 : {
1916 0 : alignment_point = 2 * L - 3 * L / 20;
1917 : }
1918 : else
1919 : {
1920 1 : alignment_point = 2 * L;
1921 : }
1922 1 : mvr2r( prevsynth + alignment_point - N, noisevect, N );
1923 :
1924 :
1925 : /* Noise addition on full band */
1926 : /* residual */
1927 1 : if ( N < L )
1928 : {
1929 1 : N_noise = ( N / 2 );
1930 217 : for ( k = 0; k < N; k++ )
1931 : {
1932 216 : noisevect[k] = ( noisevect[k] - synthesis[k] );
1933 : }
1934 : }
1935 : else
1936 : {
1937 0 : N_noise = L / 2;
1938 0 : for ( k = 0; k < L; k++ )
1939 : {
1940 0 : noisevect[k] = ( noisevect[N - L + k] - synthesis[N - L + k] );
1941 : }
1942 : }
1943 :
1944 1 : if ( HqVoicing )
1945 : {
1946 0 : for ( i = 0; i < N; i++ )
1947 : {
1948 0 : noisevect[i] *= 0.25f;
1949 : }
1950 : }
1951 :
1952 1 : kk = 0;
1953 1 : k = 0;
1954 1 : Rnd_N_noise = N_noise;
1955 :
1956 21 : while ( k < 2 * L )
1957 : {
1958 20 : if ( kk == 0 )
1959 : {
1960 10 : tmp = ( ( own_random( ni_seed_forfec ) / (PCM16_TO_FLT_FAC) *0.2f ) + 0.5f );
1961 10 : kk = 1;
1962 : }
1963 : else
1964 : {
1965 10 : tmp = ( ( own_random( ni_seed_forfec ) / (PCM16_TO_FLT_FAC) *0.3f ) + 0.7f );
1966 10 : kk = 0;
1967 : }
1968 :
1969 20 : Rnd_N_noise = (int16_t) ( (float) N_noise * tmp );
1970 :
1971 20 : sinq( (const float) EVS_PI / ( 2.0f * (float) Rnd_N_noise ), (const float) EVS_PI / ( 4.0f * (float) Rnd_N_noise ), (const int16_t) Rnd_N_noise, SS );
1972 :
1973 1305 : for ( i = 0; i < Rnd_N_noise; i++ )
1974 : {
1975 1285 : if ( k < 2 * L )
1976 : {
1977 1280 : synthesis[k] += ( noisevect[N_noise - Rnd_N_noise + i] * SS[i] + noisevect[N_noise + i] * SS[Rnd_N_noise - i - 1] ); /* *noisefact; */
1978 : }
1979 1285 : k++;
1980 : }
1981 : }
1982 :
1983 1 : if ( element_mode == EVS_MONO )
1984 : {
1985 0 : kk = 7 * L / 20;
1986 0 : p_mdct_ola = prevsynth + 37 * L / 20;
1987 : }
1988 : else
1989 : {
1990 1 : kk = NS2SA( L * FRAMES_PER_SEC, N_ZERO_MDCT_NS );
1991 1 : p_mdct_ola = old_out + kk;
1992 : }
1993 :
1994 : /* overlappadd with the ms of valid mdct of the last frame */
1995 1 : sinq( EVS_PI / ( 6.0f * L / 20.0f ), EVS_PI / ( 12.0f * L / 20.0f ), 3 * L / 20, SS );
1996 :
1997 97 : for ( k = 0; k < 3 * L / 20; k++ )
1998 : {
1999 96 : tmp = SS[k] * SS[k];
2000 96 : synthesis[k] = p_mdct_ola[k] * ( 1 - tmp ) + synthesis[k] * tmp;
2001 : }
2002 :
2003 1 : mvr2r( synthesis, synthesis + kk, 2 * L - kk );
2004 1 : mvr2r( synthesis + L, gapsynth, L );
2005 :
2006 1 : mvr2r( prevsynth + alignment_point - kk, synthesis, kk );
2007 :
2008 1 : return;
2009 : }
2010 :
2011 :
2012 : /*--------------------------------------------------------------------------
2013 : * fec_alg()
2014 : *
2015 : * Pitch based error-concealment algorithm with adaptive analysis frame
2016 : * length
2017 : *--------------------------------------------------------------------------*/
2018 :
2019 1 : static void fec_alg(
2020 : const float *prevsynth, /* i : previous synthesis */
2021 : const float *prevsynth_LP, /* i : down-sampled synthesis */
2022 : float *ecu_rec, /* o : ECU frame with Windowing/TDA */
2023 : const int16_t output_frame, /* i : Output frame length */
2024 : const int16_t N, /* i : Analysis frame length in 8 kHz (corr. max) */
2025 : const int16_t decimatefactor, /* i : Downsampling factor */
2026 : const int16_t HqVoicing, /* i : HQ voicing flag */
2027 : float *gapsynth, /* o : Buffer for crossfade to good frame */
2028 : int16_t *ni_seed_forfec, /* i/o: Random seed for picking noise segments */
2029 : const int16_t element_mode, /* i : IVAS element mode */
2030 : const float *old_out
2031 :
2032 : )
2033 : {
2034 : int16_t n, Nfft;
2035 : float sum_Tf_abs;
2036 : float Tfr[FEC_FFT_MAX_SIZE];
2037 : float Tfi[FEC_FFT_MAX_SIZE];
2038 : float Tf_abs[FEC_FFT_MAX_SIZE / 2];
2039 : float synthesis[2 * L_FRAME48k];
2040 :
2041 1 : fec_ecu_dft( prevsynth_LP, N, Tfr, Tfi, &sum_Tf_abs, Tf_abs, &Nfft, element_mode );
2042 :
2043 1 : sinusoidal_synthesis( Tfr, Tfi, Tf_abs, N, output_frame, decimatefactor, Nfft, sum_Tf_abs, synthesis, HqVoicing );
2044 :
2045 1 : fec_noise_filling( prevsynth, synthesis, output_frame, N * decimatefactor, HqVoicing, gapsynth, ni_seed_forfec, element_mode, old_out );
2046 :
2047 1 : n = (int16_t) ( (float) output_frame * N_ZERO_MDCT_NS / FRAME_SIZE_NS );
2048 1 : wtda( synthesis + ( output_frame - n ), ecu_rec, NULL, ALDO_WINDOW, ALDO_WINDOW, output_frame );
2049 :
2050 1 : return;
2051 : }
2052 :
2053 :
2054 : /*--------------------------------------------------------------------------
2055 : * hq_phase_ecu()
2056 : *
2057 : * Main routine for HQ phase ECU
2058 : *--------------------------------------------------------------------------*/
2059 :
2060 276 : static void hq_phase_ecu(
2061 : const float *prevsynth, /* i : buffer of previously synthesized signal */
2062 : float *ecu_rec, /* o : reconstructed frame in tda domain */
2063 : int16_t *time_offs, /* i/o: Sample offset for consecutive frame losses*/
2064 : float *X_sav, /* i/o: Stored spectrum of prototype frame */
2065 : int16_t *num_p, /* i/o: Number of identified peaks */
2066 : int16_t *plocs, /* i/o: Peak locations */
2067 : float *plocsi, /* i/o: Interpolated peak locations */
2068 : const float env_stab, /* i : Envelope stability parameter */
2069 : int16_t *last_fec, /* i/o: Flag for usage of pitch dependent ECU */
2070 : const int16_t prev_bfi, /* i : indicating burst frame error */
2071 : const int16_t old_is_transient[2], /* i : flags indicating previous transient frames*/
2072 : float *mag_chg_1st, /* i/o: per band magnitude modifier for transients*/
2073 : float Xavg[LGW_MAX], /* i/o: Frequency group average gain to fade to */
2074 : float *beta_mute, /* o : Factor for long-term mute */
2075 : const int16_t bwidth, /* i : Encoded bandwidth */
2076 : const int16_t output_frame, /* i : frame length */
2077 : const float pcorr, /* i : pitch correlation */
2078 : const int16_t element_mode /* i : IVAS element mode */
2079 : )
2080 : {
2081 : int16_t Lprot;
2082 : float mag_chg[LGW_MAX], ph_dith, X[L_PROT48k];
2083 : int16_t seed;
2084 : float alpha[LGW_MAX], beta[LGW_MAX];
2085 : const float *old_dec;
2086 : float noise_fac;
2087 : int16_t ph_ecu_lookahead;
2088 :
2089 276 : noise_fac = 1.0f;
2090 :
2091 276 : if ( element_mode == EVS_MONO )
2092 : {
2093 0 : ph_ecu_lookahead = NS2SA( output_frame * FRAMES_PER_SEC, PH_ECU_LOOKAHEAD_NS );
2094 : }
2095 : else
2096 : {
2097 276 : ph_ecu_lookahead = 0;
2098 : }
2099 :
2100 :
2101 276 : Lprot = ( 2 * output_frame * 4 ) / 5;
2102 :
2103 276 : if ( !prev_bfi || ( prev_bfi && *last_fec && ( *time_offs == output_frame ) ) )
2104 : {
2105 188 : if ( !( prev_bfi && *last_fec && ( element_mode == EVS_MONO ) ) )
2106 : {
2107 188 : *time_offs = 0; /* IVAS reset of offset time counter, timeoffset variable later also used to calculate burst length */
2108 : }
2109 :
2110 188 : trans_ana( prevsynth + 2 * output_frame - Lprot - *time_offs + ph_ecu_lookahead, mag_chg, &ph_dith, mag_chg_1st, output_frame, *time_offs, env_stab, *last_fec, element_mode, alpha, beta, beta_mute, Xavg );
2111 188 : spec_ana( prevsynth + 2 * output_frame - Lprot - *time_offs + ph_ecu_lookahead, plocs, plocsi, num_p, X_sav, output_frame, bwidth, element_mode, &noise_fac, pcorr );
2112 :
2113 188 : if ( prev_bfi && *last_fec )
2114 : {
2115 0 : if ( element_mode != EVS_MONO )
2116 : {
2117 0 : *time_offs = (int16_t) ( *time_offs + output_frame ); /* USAN avoid risk of internal int32_t in "+=" */
2118 0 : if ( *time_offs <= 0 )
2119 : { /* detected wrap around of st->time_offs */
2120 0 : *time_offs = (int16_t) INT16_MAX; /* keep a very high value so that the long term muting stays on */
2121 : }
2122 : }
2123 : else
2124 : {
2125 0 : *time_offs = (int16_t) ( *time_offs + output_frame ); /* EVS_MONO BE compatible, but EVS CR needed as wrap will cause burst length muting envelope instability issues */
2126 : }
2127 : }
2128 : }
2129 : else
2130 : {
2131 88 : *time_offs = (int16_t) ( *time_offs + output_frame ); /* cast added for USAN, "+=" avoided as it may creat a truncation from int to int16_t */
2132 88 : if ( *time_offs <= 0 )
2133 : {
2134 : /* detect wrap around of st->time_offs */
2135 10 : *time_offs = (int16_t) INT16_MAX; /* high value --> continued muting will ensure that the now saturated seed is not creating tones */
2136 : }
2137 :
2138 88 : trans_ana( prevsynth + 2 * output_frame - Lprot, mag_chg, &ph_dith, mag_chg_1st, output_frame, *time_offs, env_stab, 0, element_mode, alpha, beta, beta_mute, Xavg ); /* 1.0 stable-music, 0.0 speech-like */
2139 : }
2140 :
2141 276 : mvr2r( X_sav, X, Lprot );
2142 :
2143 : /* seed for own_rand2 */
2144 276 : seed = *time_offs;
2145 276 : if ( *num_p > 0 )
2146 : {
2147 268 : seed = (int16_t) ( seed + plocs[*num_p - 1] ); /* explicit cast and "+=" not used, as it "+=" may create a cast from 32 bit to 16 bit, triggering USAN */
2148 : }
2149 :
2150 276 : subst_spec( plocs, plocsi, num_p, *time_offs, X, mag_chg, ph_dith, old_is_transient, output_frame, &seed, alpha, beta, *beta_mute, Xavg, element_mode, ph_ecu_lookahead, noise_fac );
2151 :
2152 : /* reconstructed frame in tda domain */
2153 276 : old_dec = prevsynth + 2 * output_frame - NS2SA( output_frame * FRAMES_PER_SEC, N_ZERO_MDCT_NS );
2154 276 : rec_frame( X, ecu_rec, output_frame, old_dec, element_mode, num_p, plocs );
2155 :
2156 276 : return;
2157 : }
2158 :
2159 :
2160 : /*--------------------------------------------------------------------------
2161 : * hq_ecu()
2162 : *
2163 : * Main routine for HQ ECU
2164 : *--------------------------------------------------------------------------*/
2165 :
2166 277 : void hq_ecu(
2167 : const float *prevsynth, /* i : buffer of previously synthesized signal */
2168 : float *ecu_rec, /* o : reconstructed frame in tda domain */
2169 : int16_t *time_offs, /* i/o: Sample offset for consecutive frame losses*/
2170 : float *X_sav, /* i/o: Stored spectrum of prototype frame */
2171 : int16_t *num_p, /* i/o: Number of identified peaks */
2172 : int16_t *plocs, /* i/o: Peak locations */
2173 : float *plocsi, /* i/o: Interpolated peak locations */
2174 : const float env_stab, /* i : Envelope stability parameter */
2175 : int16_t *last_fec, /* i/o: Flag for usage of pitch dependent ECU */
2176 : const int16_t ph_ecu_HqVoicing, /* i : HQ Voicing flag */
2177 : int16_t *ph_ecu_active, /* i : Phase ECU active flag */
2178 : float *gapsynth, /* o : Gap synthesis for crossfade to good frame */
2179 : const int16_t prev_bfi, /* i : indicating burst frame error */
2180 : const int16_t old_is_transient[2], /* i : flags indicating previous transient frames*/
2181 : float *mag_chg_1st, /* i/o: per band magnitude modifier for transients*/
2182 : float Xavg[LGW_MAX], /* i/o: Frequency group average gain to fade to */
2183 : float *beta_mute, /* o : Factor for long-term mute */
2184 : const int16_t output_frame, /* i : frame length */
2185 : Decoder_State *st /* i/o: decoder state structure */
2186 : )
2187 : {
2188 : int16_t N;
2189 277 : float corr = 0.0f;
2190 : int16_t decimatefactor;
2191 : float prevsynth_LP[2 * L_FRAME8k];
2192 : HQ_DEC_HANDLE hHQ_core;
2193 : const float *fec_alg_input;
2194 : int16_t evs_mode_selection;
2195 : int16_t ivas_mode_selection;
2196 :
2197 277 : hHQ_core = st->hHQ_core;
2198 277 : if ( st->element_mode == EVS_MONO )
2199 : {
2200 0 : fec_alg_input = prevsynth + NS2SA( output_frame * FRAMES_PER_SEC, ACELP_LOOK_NS / 2 - PH_ECU_LOOKAHEAD_NS );
2201 : }
2202 : else
2203 : {
2204 277 : fec_alg_input = prevsynth - NS2SA( output_frame * FRAMES_PER_SEC, PH_ECU_LOOKAHEAD_NS );
2205 : }
2206 :
2207 : /* find pitch and R value */
2208 277 : if ( !( output_frame < L_FRAME16k ) )
2209 : {
2210 277 : fec_ecu_pitch( fec_alg_input, prevsynth_LP, output_frame, &N, &corr, &decimatefactor, ph_ecu_HqVoicing );
2211 : }
2212 : else
2213 : {
2214 0 : corr = 0.0f;
2215 0 : decimatefactor = 4;
2216 0 : N = output_frame / 4;
2217 : }
2218 :
2219 1 : evs_mode_selection = ( st->total_brate >= 48000 && ( output_frame >= L_FRAME16k && !prev_bfi && ( !old_is_transient[0] || old_is_transient[1] ) &&
2220 555 : ( ph_ecu_HqVoicing || ( ( ( hHQ_core->env_stab_plc > 0.5 ) && ( corr < 0.6 ) ) || ( hHQ_core->env_stab_plc < 0.5 && ( corr > 0.85 ) ) ) ) ) ) ||
2221 277 : ( st->total_brate < 48000 && ( ( ph_ecu_HqVoicing || corr > 0.85 ) && !prev_bfi && ( !old_is_transient[0] || old_is_transient[1] ) ) );
2222 :
2223 277 : ivas_mode_selection = ( N < PH_ECU_N_LIMIT ) || ( corr < PH_ECU_CORR_LIMIT );
2224 :
2225 277 : if ( ( ( st->element_mode == EVS_MONO ) && evs_mode_selection ) ||
2226 277 : ( ( st->element_mode != EVS_MONO ) && evs_mode_selection && ivas_mode_selection ) )
2227 :
2228 : {
2229 1 : fec_alg( fec_alg_input, prevsynth_LP, ecu_rec, output_frame, N, decimatefactor, ph_ecu_HqVoicing, gapsynth, &hHQ_core->ni_seed_forfec, st->element_mode, st->hHQ_core->old_out );
2230 1 : *last_fec = 1;
2231 1 : *ph_ecu_active = 0;
2232 1 : *time_offs = output_frame;
2233 : }
2234 : else
2235 : {
2236 276 : hq_phase_ecu( prevsynth - NS2SA( output_frame * FRAMES_PER_SEC, PH_ECU_LOOKAHEAD_NS ), ecu_rec, time_offs, X_sav, num_p, plocs, plocsi, env_stab, last_fec, prev_bfi, old_is_transient, mag_chg_1st, Xavg, beta_mute, st->bwidth, output_frame, corr, st->element_mode );
2237 :
2238 276 : *last_fec = 0;
2239 276 : *ph_ecu_active = 1;
2240 : }
2241 :
2242 277 : return;
2243 : }
|
