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 720 : 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 118992 : for ( i = 0, j = N - 1; i < N; i++, j-- )
113 : {
114 118272 : y[i] = x1[i] * x2[j];
115 : }
116 :
117 720 : return;
118 : }
119 :
120 :
121 : /*-------------------------------------------------------------------*
122 : * fft_spec2()
123 : *
124 : * Square magnitude of fft spectrum
125 : *-------------------------------------------------------------------*/
126 :
127 1080 : 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 354816 : for ( i = 1, j = N - 1; i < N / 2; i++, j-- )
135 : {
136 353736 : x[i] = x[i] * x[i] + x[j] * x[j];
137 : }
138 :
139 1080 : x[0] *= x[0];
140 1080 : x[N / 2] *= x[N / 2];
141 :
142 1080 : 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 171150 : 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 171150 : const float *sincos = sincos_t_ext + 128;
159 171150 : int16_t seed2 = seed;
160 171150 : own_random( &seed2 );
161 :
162 171150 : if ( seed2 & 0x40 )
163 : {
164 85572 : *sin_F = sincos[seed2 >> 8];
165 : }
166 : else
167 : {
168 85578 : *sin_F = -sincos[seed2 >> 8];
169 : }
170 :
171 171150 : if ( seed2 & 0x80 )
172 : {
173 85530 : *cos_F = sincos[-( seed2 >> 8 )];
174 : }
175 : else
176 : {
177 85620 : *cos_F = -sincos[-( seed2 >> 8 )];
178 : }
179 :
180 171150 : 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 6315 : 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 6315 : pY = y_re;
214 6315 : y_m1_re = *pY++;
215 6315 : y_0_re = *pY++;
216 6315 : y_p1_re = *pY++;
217 :
218 : /* Same for imaginary parts - note reverse order from FFT */
219 6315 : pY = y_im;
220 6315 : y_p1_im = *pY++;
221 6315 : y_0_im = *pY++;
222 6315 : y_m1_im = *pY++;
223 :
224 : /* prepare numerator real and imaginary parts*/
225 6315 : N_re = y_m1_re - y_p1_re;
226 6315 : N_im = y_m1_im - y_p1_im;
227 :
228 : /* prepare denominator real and imaginary parts */
229 :
230 6315 : D_re = 2 * y_0_re - y_m1_re - y_p1_re;
231 6315 : D_im = 2 * y_0_im - y_m1_im - y_p1_im;
232 :
233 : /* REAL part of complex division */
234 6315 : numer = N_re * D_re + N_im * D_im;
235 6315 : denom = D_re * D_re + D_im * D_im;
236 :
237 6315 : test();
238 6315 : if ( numer != 0 && denom != 0 )
239 : {
240 6315 : 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 6315 : return posi;
248 : }
249 :
250 :
251 : /*------------------------------------------------------------------*
252 : * trans_ana()
253 : *
254 : * Transient analysis
255 : *------------------------------------------------------------------*/
256 :
257 624 : 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 624 : 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 624 : burst_len = time_offs / output_frame + 1;
287 :
288 624 : set_s( att_always, 0, LGW_MAX );
289 624 : *ph_dith = 0.0f;
290 :
291 : /* softly shift attenuation just a bit later for estimated "stable" music_content */
292 624 : burst_phdith_thresh = BURST_PHDITH_THRESH + (int16_t) ( est_mus_content * 1.0f + 0.5f );
293 624 : burst_att_thresh = BURST_ATT_THRESH + (int16_t) ( est_mus_content * 1.0f + 0.5f );
294 624 : att_per_frame = (float) ( ATT_PER_FRAME - (int16_t) ( est_mus_content * 1.0f + 0.5f ) ); /* only slighty less att for music */
295 624 : att_per_frame *= 0.1f;
296 :
297 624 : if ( burst_len > burst_phdith_thresh )
298 : {
299 : /* increase degree of dither */
300 216 : *ph_dith = PHASE_DITH * min( 1.0f, ( (float) burst_len - (float) burst_phdith_thresh ) / (float) BURST_PHDITH_RAMPUP_LEN );
301 : }
302 :
303 624 : att_degree = 0;
304 624 : if ( burst_len > burst_att_thresh )
305 : {
306 222 : set_s( att_always, 1, LGW_MAX );
307 :
308 : /* increase degree of attenuation */
309 222 : if ( burst_len - burst_att_thresh <= PH_ECU_MUTE_START )
310 : {
311 90 : att_degree = (float) ( burst_len - burst_att_thresh ) * att_per_frame;
312 : }
313 : else
314 : {
315 132 : 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 624 : Lprot = ( 2 * output_frame * 4 ) / 5; /* 4/5==1024/1280, keep mult within short */
320 624 : Ltrana = Lprot / QUOT_LPR_LTR;
321 624 : Ltrana_2 = Ltrana / 2;
322 :
323 624 : if ( output_frame == L_FRAME48k )
324 : {
325 459 : w_hamm = w_hamm48k_2;
326 459 : Lgw = LGW48k;
327 : }
328 165 : else if ( output_frame == L_FRAME32k )
329 : {
330 165 : w_hamm = w_hamm32k_2;
331 165 : LtranaLog = L_TRANA_LOG32k;
332 165 : 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 624 : if ( burst_len <= 1 || ( burst_len == 2 && last_fec ) )
342 : {
343 360 : set_f( alpha, 1.0f, LGW_MAX );
344 360 : set_f( beta, 0.0f, LGW_MAX );
345 360 : *beta_mute = BETA_MUTE_FAC_INI;
346 :
347 : /* apply hamming window */
348 360 : v_mult( xfp, w_hamm, xfp_left, Ltrana_2 );
349 360 : mult_rev2( xfp + Ltrana_2, w_hamm, xfp_left + Ltrana_2, Ltrana_2 );
350 :
351 360 : xfp_ = xfp + Lprot - Ltrana;
352 360 : v_mult( xfp_, w_hamm, xfp_right, Ltrana_2 );
353 360 : mult_rev2( xfp_ + Ltrana_2, w_hamm, xfp_right + Ltrana_2, Ltrana_2 );
354 :
355 : /* spectrum */
356 360 : if ( output_frame == L_FRAME48k )
357 : {
358 204 : fft3( xfp_left, xfp_left, Ltrana );
359 204 : fft3( xfp_right, xfp_right, Ltrana );
360 : }
361 : else
362 : {
363 156 : fft_rel( xfp_left, Ltrana, LtranaLog );
364 156 : fft_rel( xfp_right, Ltrana, LtranaLog );
365 : }
366 :
367 : /* square representation */
368 360 : fft_spec2( xfp_left, Ltrana );
369 360 : fft_spec2( xfp_right, Ltrana );
370 :
371 : /* band powers in frequency groups
372 : exclude bin at 0 and at EVS_PI from calculation */
373 360 : xfp_left[Ltrana_2] = 0.0f;
374 360 : xfp_right[Ltrana_2] = 0.0f;
375 : }
376 :
377 5451 : for ( k = 0; k < Lgw; k++ )
378 : {
379 4827 : if ( burst_len <= 1 || ( burst_len == 2 && last_fec ) )
380 : {
381 2724 : gr_pow_left[k] = sum_f( xfp_left + gw[k], gw[k + 1] - gw[k] );
382 2724 : 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 2724 : gr_pow_left[k] += FLT_MIN; /* otherwise div by zero may occur */
386 2724 : gr_pow_right[k] += FLT_MIN;
387 :
388 2724 : Xavg[k] = (float) ( sqrt( 0.5f * ( gr_pow_left[k] + gr_pow_right[k] ) / (float) ( gw[k + 1] - gw[k] ) ) );
389 :
390 2724 : 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 2724 : tr_dec[k] = ( grp_pow_chg > THRESH_TR_LIN ) || ( grp_pow_chg < THRESH_TR_LIN_INV );
395 :
396 : /* magnitude modification */
397 2724 : if ( tr_dec[k] || att_always[k] )
398 : {
399 204 : att_val = min( MAX_INCREASE_GRPOW_LIN, grp_pow_chg );
400 204 : att_val = (float) sqrt( att_val );
401 204 : mag_chg_1st[k] = att_val;
402 204 : mag_chg[k] = att_val;
403 : }
404 : else
405 : {
406 2520 : mag_chg_1st[k] = 1.0f;
407 2520 : mag_chg[k] = 1.0f;
408 : }
409 : }
410 : else
411 : {
412 2103 : if ( burst_len < OFF_FRAMES_LIMIT )
413 : {
414 1575 : mag_chg[k] = mag_chg_1st[k] * (float) pow( 10.0, -att_degree / 20.0 );
415 : }
416 : else
417 : {
418 528 : mag_chg[k] = 0;
419 : }
420 :
421 2103 : if ( element_mode != EVS_MONO )
422 : {
423 2103 : 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 180 : *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 2103 : alpha[k] = mag_chg[k];
437 2103 : beta[k] = (float) ( sqrt( 1.0f - SQR( alpha[k] ) ) * *beta_mute );
438 2103 : if ( k >= LGW32k - 1 )
439 : {
440 519 : beta[k] *= 0.1f;
441 : }
442 1584 : else if ( k >= LGW16k - 1 )
443 : {
444 264 : beta[k] *= 0.5f;
445 : }
446 : }
447 : }
448 :
449 624 : return;
450 : }
451 :
452 :
453 : /*------------------------------------------------------------------*
454 : * peakfinder()
455 : *
456 : * Peak-picking algorithm
457 : *------------------------------------------------------------------*/
458 :
459 2418 : 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 2418 : int16_t k, i, len, tempLoc = 0, foundPeak, ii, xInd;
471 : int16_t *ind, indarr[L_PROT48k_2 + 1], peakLoc[MAX_PLOCS];
472 :
473 2418 : ind = indarr;
474 :
475 : /* Find derivative */
476 2418 : v_sub( x0 + 1, x0, dx0, len0 - 1 );
477 :
478 : /* This is so we find the first of repeated values */
479 280122 : for ( i = 0; i < len0 - 1; i++ )
480 : {
481 277704 : if ( dx0[i] == 0.0f )
482 : {
483 18615 : dx0[i] = -1.0e-12f;
484 : }
485 : }
486 :
487 : /* Find where the derivative changes sign
488 : Include endpoints in potential peaks and valleys */
489 2418 : k = 0;
490 :
491 2418 : if ( endpoints )
492 : {
493 360 : x[k] = x0[0];
494 360 : ind[k++] = 0;
495 : }
496 :
497 277704 : for ( i = 1; i < len0 - 1; i++ )
498 : {
499 275286 : if ( dx0[i - 1] * dx0[i] < 0 )
500 : {
501 139581 : ind[k] = i;
502 139581 : x[k++] = x0[i];
503 : }
504 : }
505 :
506 2418 : if ( endpoints )
507 : {
508 360 : ind[k] = len0 - 1;
509 360 : x[k++] = x0[len0 - 1];
510 : }
511 : /* x only has the peaks, valleys, and endpoints */
512 2418 : len = k;
513 2418 : minimum( x, len, &minMag );
514 :
515 2418 : if ( ( len > 2 ) || ( !endpoints && ( len > 0 ) ) )
516 : {
517 : /* Set initial parameters for loop */
518 2406 : tempMag = minMag;
519 2406 : foundPeak = 0;
520 2406 : leftMin = minMag;
521 :
522 2406 : 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 360 : if ( x[0] >= x[1] )
530 : {
531 105 : ii = -1;
532 105 : 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 255 : ii = 0;
543 255 : 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 2046 : ii = -1; /* First point is a peak */
555 2046 : if ( len >= 2 )
556 : {
557 2046 : if ( x[1] >= x[0] )
558 : {
559 357 : ii = 0; /* First point is a valley, skip it */
560 : }
561 : }
562 : }
563 2406 : *cInd = 0;
564 :
565 : /* Loop through extrema which should be peaks and then valleys */
566 71568 : while ( ii < len - 1 )
567 : {
568 70527 : 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 70527 : if ( foundPeak )
573 : {
574 11718 : tempMag = minMag;
575 11718 : foundPeak = 0;
576 : }
577 :
578 : /* Make sure we don't iterate past the length of our vector */
579 70527 : if ( ii == len - 1 )
580 : {
581 1365 : 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 69162 : if ( ( x[ii] > tempMag ) && ( x[ii] > leftMin + sel ) )
587 : {
588 13326 : tempLoc = ii;
589 13326 : tempMag = x[ii];
590 : }
591 :
592 69162 : ii++; /* Move onto the valley */
593 :
594 : /* Come down at least sel from peak */
595 69162 : if ( !foundPeak && ( tempMag > sel + x[ii] ) )
596 : {
597 12069 : foundPeak = 1; /* We have found a peak */
598 12069 : leftMin = x[ii];
599 12069 : peakLoc[*cInd] = tempLoc; /* Add peak to index */
600 12069 : peakMag[*cInd] = tempMag;
601 12069 : ( *cInd )++;
602 : }
603 57093 : else if ( x[ii] < leftMin ) /* New left minimum */
604 : {
605 10119 : leftMin = x[ii];
606 : }
607 : }
608 :
609 : /* Check end point */
610 2406 : if ( x[len - 1] > tempMag && x[len - 1] > leftMin + sel )
611 : {
612 414 : peakLoc[*cInd] = len - 1;
613 414 : peakMag[*cInd] = x[len - 1];
614 414 : ( *cInd )++;
615 : }
616 1992 : else if ( !foundPeak && tempMag > minMag ) /* Check if we still need to add the last point */
617 : {
618 12 : peakLoc[*cInd] = tempLoc;
619 12 : peakMag[*cInd] = tempMag;
620 12 : ( *cInd )++;
621 : }
622 :
623 : /* Create output */
624 14901 : for ( i = 0; i < *cInd; i++ )
625 : {
626 12495 : plocs[i] = ind[peakLoc[i]];
627 : }
628 : }
629 : else
630 : {
631 12 : 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 12 : *cInd = 0;
650 : }
651 : }
652 :
653 2418 : return;
654 : }
655 :
656 :
657 : /*-------------------------------------------------------------------*
658 : * imax_pos()
659 : *
660 : * Get interpolated maximum position
661 : *-------------------------------------------------------------------*/
662 :
663 : /*! r: interpolated maximum position */
664 5061 : 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 5061 : y1 = y[0];
672 5061 : y2 = y[1];
673 5061 : y3 = y[2];
674 5061 : y3_y1 = y3 - y1;
675 5061 : ftmp_den1 = ( y1 + y3 - 2 * y2 );
676 5061 : ftmp_den2 = ( 4 * y2 - 2 * y1 - 2 * y3 );
677 :
678 5061 : if ( ftmp_den2 == 0.0f || ftmp_den1 == 0.0f )
679 : {
680 0 : return ( 0.0f ); /* early exit with left-most value */
681 : }
682 :
683 5061 : y2i = -0.125f * SQR( y3_y1 ) / ( ftmp_den1 ) + y2;
684 : /* their corresponding normalized locations */
685 5061 : posi = y3_y1 / ( ftmp_den2 );
686 : /* Interpolated maxima if locations are not within [-1,1], calculated extrema are ignored */
687 5061 : if ( posi >= 1.0f || posi <= -1.0f )
688 : {
689 3 : posi = y3 > y1 ? 1.0f : -1.0f;
690 : }
691 : else
692 : {
693 5058 : if ( y1 >= y2i )
694 : {
695 3 : posi = ( y1 > y3 ) ? -1.0f : 1.0f;
696 : }
697 5055 : else if ( y3 >= y2i )
698 : {
699 0 : posi = 1.0f;
700 : }
701 : }
702 :
703 5061 : return posi + 1.0f;
704 : }
705 :
706 :
707 : /*-------------------------------------------------------------------*
708 : * spec_ana()
709 : *
710 : * Spectral analysis
711 : *-------------------------------------------------------------------*/
712 :
713 360 : 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 360 : 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 360 : 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 360 : Lprot = 2 * output_frame * L_PROT32k / 1280;
742 360 : Lprot2_1 = Lprot / 2 + 1;
743 :
744 360 : if ( output_frame == L_FRAME48k )
745 : {
746 204 : w_hamm = w_hamm_sana48k_2;
747 204 : hamm_len2 = L_PROT_HAMM_LEN2_48k;
748 : }
749 156 : else if ( output_frame == L_FRAME32k )
750 : {
751 156 : w_hamm = w_hamm_sana32k_2;
752 156 : hamm_len2 = L_PROT_HAMM_LEN2_32k;
753 156 : 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 360 : mvr2r( prevsynth + hamm_len2, xfp + hamm_len2, Lprot - 2 * hamm_len2 );
764 360 : 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 360 : window_corr = w_hamm[0];
772 360 : window_corr_step = w_hamm[0] / hamm_len2;
773 89064 : for ( i = 0; i < hamm_len2; i++ )
774 : {
775 88704 : xfp[i] = prevsynth[i] * ( w_hamm[i] - window_corr );
776 88704 : xfp[Lprot - i - 1] = prevsynth[Lprot - i - 1] * ( w_hamm[i] - window_corr );
777 88704 : window_corr -= window_corr_step;
778 : }
779 : }
780 :
781 : /* Spectrum */
782 360 : if ( output_frame == L_FRAME48k )
783 : {
784 204 : fft3( xfp, xfp, Lprot );
785 : }
786 : else
787 : {
788 156 : fft_rel( xfp, Lprot, LprotLog2 );
789 : }
790 :
791 : /* Apply zeroing of non-coded FFT spectrum */
792 360 : if ( output_frame > inner_frame_tbl[bwidth] )
793 : {
794 57 : stop_band_start = 128 << bwidth;
795 57 : stop_band_length = Lprot - ( stop_band_start << 1 );
796 57 : stop_band_start = stop_band_start + 1;
797 57 : set_f( xfp + stop_band_start, 0, stop_band_length );
798 : }
799 :
800 360 : mvr2r( xfp, X_sav, Lprot );
801 :
802 : /* Magnitude representation */
803 360 : fft_spec2( xfp, Lprot );
804 :
805 237264 : for ( i = 0; i < Lprot2_1; i++ )
806 : {
807 236904 : xfp[i] = (float) sqrt( (double) xfp[i] );
808 : }
809 :
810 : /* Find maxima */
811 360 : maximum( xfp, Lprot2_1, &Xmax );
812 360 : minimum( xfp, Lprot2_1, &Xmin );
813 360 : if ( element_mode == EVS_MONO )
814 : {
815 0 : sel = ( Xmax - Xmin ) * ( 1.0f - PFIND_SENS );
816 : }
817 : else
818 : {
819 360 : sel = ( Xmax - Xmin ) * ( 1.0f - ST_PFIND_SENS );
820 : }
821 :
822 360 : 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 360 : if ( element_mode != EVS_MONO && *num_plocs > 50 && pcorr < 0.6f )
826 : {
827 21 : *num_plocs = 0;
828 : }
829 :
830 360 : 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 360 : Lprot2p1 = Lprot / 2 + 1;
853 :
854 : /* Refine peaks */
855 360 : pPlocsi = plocsi;
856 360 : pPlocs = plocs;
857 360 : 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 360 : if ( n > 0 && *pPlocs == 0 ) /* Very 1st peak position possible to have a peak at 0/DC index position. */
862 : {
863 24 : *pPlocsi++ = *pPlocs + imax_pos( &xfp[*pPlocs] );
864 24 : pPlocs++;
865 24 : n = n - 1;
866 : }
867 :
868 360 : if ( n > 0 && *pPlocs == 1 ) /* Also 2nd peak position uses DC which makes jacobsen unsuitable. */
869 : {
870 21 : *pPlocsi++ = *pPlocs - 1 + imax_pos( &xfp[*pPlocs - 1] );
871 21 : currPlocs = *pPlocs++;
872 21 : n = n - 1;
873 : }
874 :
875 : /* All remaining peaks except the very last two possible integer positions */
876 360 : currPlocs = *pPlocs++;
877 360 : 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 360 : if ( n > 0 )
881 : {
882 336 : nJacob = n;
883 336 : if ( sub( endPlocs, plocs[sub( *num_plocs, 1 )] ) <= 0 )
884 : {
885 0 : nJacob = sub( nJacob, 1 );
886 : }
887 :
888 6651 : for ( k = 0; k < nJacob; k++ )
889 : {
890 6315 : *pPlocsi++ = currPlocs + imax2_jacobsen_mag( &( X_sav[currPlocs - 1] ), &( X_sav[Lprot - 1 - currPlocs] ) );
891 6315 : currPlocs = *pPlocs++;
892 : }
893 336 : 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 360 : 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 360 : if ( *num_plocs < 3 && *num_plocs > 0 )
924 : {
925 15 : sig = sum_f( xfp, Lprot2_1 ) + EPSILON;
926 :
927 : /*excluding peaks and neighboring bins*/
928 36 : for ( i = 0; i < *num_plocs; i++ )
929 : {
930 21 : st_point = max( 0, plocs[i] - DELTA_CORR );
931 21 : end_point = min( Lprot2_1 - 1, plocs[i] + DELTA_CORR );
932 21 : set_f( &xfp[st_point], 0.0f, end_point - st_point + 1 );
933 : }
934 15 : noise = sum_f( xfp, Lprot2_1 ) + EPSILON;
935 15 : nsr = noise / sig;
936 :
937 15 : if ( nsr < 0.03f )
938 : {
939 0 : *noise_fac = 0.5f;
940 : }
941 : else
942 : {
943 15 : *noise_fac = 1.0f;
944 : }
945 : }
946 : }
947 :
948 360 : return;
949 : }
950 :
951 : /*-------------------------------------------------------------------*
952 : * subst_spec()
953 : *
954 : * Substitution spectrum calculation
955 : *-------------------------------------------------------------------*/
956 :
957 624 : 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 624 : sincos = sincos_t_ext + 128;
990 624 : Lprot = (int16_t) ( L_PROT32k * output_frame / 640 );
991 624 : Lprot_1 = 1.0f / Lprot;
992 624 : Lecu = output_frame * 2;
993 :
994 : /* Correction phase of the identified peaks */
995 624 : if ( is_trans[0] || is_trans[1] )
996 : {
997 39 : *num_plocs = 0;
998 : }
999 : else
1000 : {
1001 585 : 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 21570 : for ( m = 0; m < *num_plocs; m++ )
1003 : {
1004 20985 : corr_phase[m] = plocsi[m] * tmp;
1005 : }
1006 : }
1007 624 : one_peak_flag_mask = 1; /* all ones mask -> keep */
1008 624 : if ( element_mode != EVS_MONO )
1009 : {
1010 624 : if ( ( *num_plocs > 0 ) && sub( *num_plocs, 3 ) < 0 )
1011 : {
1012 15 : one_peak_flag_mask = noise_fac; /* all zeroes mask -> zero */
1013 : }
1014 624 : if ( *num_plocs == 0 )
1015 : {
1016 60 : X[0] = 0; /* reset DC if there are no peaks */
1017 60 : X[shr( Lprot, 1 )] = 0; /* also reset fs/2 if there are no peaks */
1018 : }
1019 : }
1020 :
1021 624 : i = 1;
1022 624 : k = 0;
1023 624 : im_ind = Lprot - 1;
1024 21609 : for ( m = 0; m < *num_plocs; m++ )
1025 : {
1026 20985 : delta_corr_dn = DELTA_CORR;
1027 20985 : delta_corr_up = DELTA_CORR;
1028 :
1029 20985 : if ( m > 0 )
1030 : {
1031 20421 : delta_tmp = ( plocs[m] - plocs[m - 1] - 1 ) / 2;
1032 20421 : if ( delta_tmp < DELTA_CORR )
1033 : {
1034 17235 : delta_corr_dn = delta_tmp;
1035 : }
1036 : }
1037 :
1038 20985 : if ( m < *num_plocs - 1 )
1039 : {
1040 20421 : delta_tmp = ( plocs[m + 1] - plocs[m] - 1 ) / 2;
1041 20421 : if ( delta_tmp < DELTA_CORR )
1042 : {
1043 17235 : delta_corr_up = delta_tmp;
1044 : }
1045 : }
1046 :
1047 : /* Input Xph */
1048 75546 : while ( i < plocs[m] - delta_corr_dn )
1049 : {
1050 54561 : *seed = own_random( seed );
1051 :
1052 54561 : if ( *seed & 0x40 )
1053 : {
1054 27516 : sin_F = sincos[*seed >> 8];
1055 : }
1056 : else
1057 : {
1058 27045 : sin_F = -sincos[*seed >> 8];
1059 : }
1060 :
1061 54561 : if ( *seed & 0x80 )
1062 : {
1063 27582 : cos_F = sincos[-( *seed >> 8 )];
1064 : }
1065 : else
1066 : {
1067 26979 : cos_F = -sincos[-( *seed >> 8 )];
1068 : }
1069 :
1070 54561 : 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 54561 : tmp = one_peak_flag_mask * ( X[i] * cos_F - X[im_ind] * sin_F );
1078 54561 : X[im_ind] = one_peak_flag_mask * ( X[i] * sin_F + X[im_ind] * cos_F );
1079 : }
1080 :
1081 54561 : if ( alpha[k] < 1.0f )
1082 : {
1083 19230 : *seed = rand_phase( *seed, &sin_F, &cos_F );
1084 19230 : X[i] = alpha[k] * tmp + beta[k] * Xavg[k] * cos_F;
1085 19230 : X[im_ind] = alpha[k] * X[im_ind] + beta[k] * Xavg[k] * sin_F;
1086 : }
1087 : else
1088 : {
1089 35331 : X[i] = mag_chg[k] * tmp;
1090 35331 : X[im_ind] *= mag_chg[k];
1091 : }
1092 54561 : i++;
1093 54561 : im_ind--;
1094 54561 : if ( i >= ivas_gwlpr[k + 1] )
1095 : {
1096 504 : k++;
1097 : }
1098 : }
1099 :
1100 20985 : e = plocs[m] + delta_corr_up;
1101 20985 : if ( e > Lprot / 2 - 1 )
1102 : {
1103 0 : e = Lprot / 2 - 1;
1104 : }
1105 :
1106 20985 : Xph = corr_phase[m];
1107 : /* extract fractional phase integer index in the range [0...1023] */
1108 20985 : Xph_short = (int16_t) ( 0x000003ff & (int32_t) ( ( Xph * 512 ) / EVS_PI ) );
1109 :
1110 :
1111 20985 : if ( Xph_short >= 512 )
1112 : {
1113 10533 : sin_F = -sincos_t_ext[Xph_short - 512];
1114 10533 : if ( Xph_short < 768 )
1115 : {
1116 5289 : cos_F = -sincos_t_ext[Xph_short - 512 + 256];
1117 : }
1118 : else
1119 : {
1120 5244 : cos_F = sincos_t_ext[-Xph_short + 1024 + 256];
1121 : }
1122 : }
1123 : else
1124 : {
1125 10452 : sin_F = sincos_t_ext[Xph_short];
1126 10452 : if ( Xph_short < 256 )
1127 : {
1128 5505 : cos_F = sincos_t_ext[Xph_short + 256];
1129 : }
1130 : else
1131 : {
1132 4947 : cos_F = -sincos_t_ext[-Xph_short + 256 + 512];
1133 : }
1134 : }
1135 :
1136 149367 : while ( i <= e )
1137 : {
1138 128382 : mag_chg_local = mag_chg[k];
1139 :
1140 128382 : if ( ph_dith != 0.0f )
1141 : {
1142 : /* Call phase randomization only when needed */
1143 74196 : Xph = corr_phase[m];
1144 74196 : *seed = own_random( seed );
1145 74196 : 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 74196 : 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 74196 : mag_chg_local *= 0.5f + ( 1.0f - ( 1.0f / PHASE_DITH ) * ph_dith ) * 0.5f;
1152 : }
1153 :
1154 74196 : Xph_short = (int16_t) ( ( (int32_t) ( ( Xph * 512 ) / EVS_PI ) ) & 0x000003ff );
1155 :
1156 74196 : if ( Xph_short >= 512 )
1157 : {
1158 37062 : sin_F = -sincos_t_ext[Xph_short - 512];
1159 37062 : if ( Xph_short < 768 )
1160 : {
1161 18222 : cos_F = -sincos_t_ext[Xph_short - 512 + 256];
1162 : }
1163 : else
1164 : {
1165 18840 : cos_F = sincos_t_ext[-Xph_short + 1024 + 256];
1166 : }
1167 : }
1168 : else
1169 : {
1170 37134 : sin_F = sincos_t_ext[Xph_short];
1171 37134 : if ( Xph_short < 256 )
1172 : {
1173 18501 : cos_F = sincos_t_ext[Xph_short + 256];
1174 : }
1175 : else
1176 : {
1177 18633 : cos_F = -sincos_t_ext[-Xph_short + 256 + 512];
1178 : }
1179 : }
1180 : }
1181 :
1182 128382 : tmp = ( X[i] * cos_F - X[im_ind] * sin_F );
1183 128382 : X[im_ind] = ( X[i] * sin_F + X[im_ind] * cos_F );
1184 128382 : if ( alpha[k] < 1.0f )
1185 : {
1186 76338 : alpha_local = mag_chg_local;
1187 76338 : beta_local = (float) ( beta_mute * sqrt( 1.0f - SQR( alpha_local ) ) );
1188 76338 : if ( k >= LGW32k - 1 )
1189 : {
1190 29637 : beta_local *= 0.1f;
1191 : }
1192 46701 : else if ( k >= LGW16k - 1 )
1193 : {
1194 23199 : beta_local *= 0.5f;
1195 : }
1196 :
1197 76338 : *seed = rand_phase( *seed, &sin_F, &cos_F );
1198 76338 : X[i] = alpha_local * tmp + beta_local * Xavg[k] * cos_F;
1199 76338 : X[im_ind] = alpha_local * X[im_ind] + beta_local * Xavg[k] * sin_F;
1200 : }
1201 : else
1202 : {
1203 52044 : X[i] = mag_chg_local * tmp;
1204 52044 : X[im_ind] *= mag_chg_local;
1205 : }
1206 :
1207 128382 : i++;
1208 128382 : im_ind--;
1209 128382 : if ( i >= ivas_gwlpr[k + 1] )
1210 : {
1211 2637 : k++;
1212 : }
1213 : }
1214 : }
1215 :
1216 254049 : while ( i < Lprot / 2 )
1217 : {
1218 253425 : *seed = own_random( seed );
1219 :
1220 253425 : if ( *seed & 0x40 )
1221 : {
1222 126597 : sin_F = sincos[*seed >> 8];
1223 : }
1224 : else
1225 : {
1226 126828 : sin_F = -sincos[*seed >> 8];
1227 : }
1228 :
1229 253425 : if ( *seed & 0x80 )
1230 : {
1231 126654 : cos_F = sincos[-( *seed >> 8 )];
1232 : }
1233 : else
1234 : {
1235 126771 : cos_F = -sincos[-( *seed >> 8 )];
1236 : }
1237 :
1238 :
1239 253425 : 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 253425 : tmp = one_peak_flag_mask * ( X[i] * cos_F - X[im_ind] * sin_F );
1247 253425 : X[im_ind] = one_peak_flag_mask * ( X[i] * sin_F + X[im_ind] * cos_F );
1248 : }
1249 :
1250 253425 : if ( alpha[k] < 1.0f )
1251 : {
1252 75582 : *seed = rand_phase( *seed, &sin_F, &cos_F );
1253 75582 : X[i] = alpha[k] * tmp + beta[k] * Xavg[k] * cos_F;
1254 75582 : X[im_ind] = alpha[k] * X[im_ind] + beta[k] * Xavg[k] * sin_F;
1255 75582 : im_ind--;
1256 : }
1257 : else
1258 : {
1259 177843 : X[i] = mag_chg[k] * tmp;
1260 177843 : X[im_ind--] *= mag_chg[k];
1261 : }
1262 253425 : i++;
1263 :
1264 253425 : if ( i >= ivas_gwlpr[k + 1] )
1265 : {
1266 1686 : k++;
1267 : }
1268 : }
1269 :
1270 624 : return;
1271 : }
1272 :
1273 : /*--------------------------------------------------------------------------
1274 : * rec_wtda()
1275 : *
1276 : * Windowing and TDA of reconstructed frame
1277 : *--------------------------------------------------------------------------*/
1278 :
1279 624 : 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 624 : 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 624 : 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 624 : if ( output_frame == L_FRAME48k )
1312 : {
1313 459 : w_hamm = w_hamm_sana48k_2;
1314 459 : hamm_len2 = L_PROT_HAMM_LEN2_48k;
1315 : }
1316 165 : else if ( output_frame == L_FRAME32k )
1317 : {
1318 165 : w_hamm = w_hamm_sana32k_2;
1319 165 : 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 624 : if ( element_mode != EVS_MONO && *num_p > 0 && plocs[0] > 3 )
1328 : {
1329 : /* Perform inverse windowing of hammrect */
1330 264 : pX_start = X;
1331 264 : pX_end = X + Lprot - 1;
1332 68520 : for ( i = 0; i < hamm_len2; i++ )
1333 : {
1334 68256 : tmp = 1.0f / *w_hamm;
1335 68256 : *pX_start *= tmp;
1336 68256 : *pX_end *= tmp;
1337 68256 : pX_start++;
1338 68256 : pX_end--;
1339 68256 : w_hamm++;
1340 : }
1341 : }
1342 :
1343 : /* extract reconstructed frame with aldo window */
1344 624 : timesh = NS2SA( output_frame * FRAMES_PER_SEC, N_ZERO_MDCT_NS ) - ( 2 * output_frame - Lprot ) / 2;
1345 :
1346 624 : set_f( xsubst_, 0.0f, 2 * output_frame - Lprot + timesh );
1347 624 : 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 624 : if ( element_mode != EVS_MONO )
1351 : {
1352 624 : mvr2r( old_dec, xsubst_ + NS2SA( output_frame * FRAMES_PER_SEC, N_ZERO_MDCT_NS ), copy_len ); /* also need to scale to Q0 ?? */
1353 624 : pOld = old_dec + copy_len;
1354 624 : pNew = xsubst_ + copy_len + NS2SA( output_frame * FRAMES_PER_SEC, N_ZERO_MDCT_NS );
1355 624 : sinq( EVS_PI / ( ola_len * 2 ), 0.0f, ola_len, xfwin );
1356 624 : v_mult( xfwin, xfwin, xfwin, ola_len ); /* xfwin = sin^2 of 0..pi/4 */
1357 624 : pOldW = xfwin + ola_len - 1;
1358 624 : pNewW = xfwin;
1359 45006 : for ( i = 0; i < ola_len; i++ )
1360 : {
1361 44382 : *pNew = *pOld * *pOldW + *pNew * *pNewW;
1362 44382 : pOld += 1;
1363 44382 : pNew += 1;
1364 44382 : pOldW -= 1;
1365 44382 : 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 624 : wtda( xsubst_ + output_frame, ecu_rec, NULL, ALDO_WINDOW, ALDO_WINDOW, output_frame );
1384 :
1385 624 : return;
1386 : }
1387 :
1388 :
1389 : /*--------------------------------------------------------------------------
1390 : * rec_frame()
1391 : *
1392 : * Frame reconstruction
1393 : *--------------------------------------------------------------------------*/
1394 :
1395 624 : 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 624 : int16_t Lprot, LprotLog2 = 0;
1406 :
1407 624 : Lprot = 2 * output_frame * L_PROT32k / 1280;
1408 :
1409 624 : if ( output_frame == L_FRAME48k )
1410 : {
1411 459 : LprotLog2 = 9;
1412 : }
1413 165 : else if ( output_frame == L_FRAME32k )
1414 : {
1415 165 : LprotLog2 = 10;
1416 : }
1417 : else
1418 : {
1419 0 : LprotLog2 = 9;
1420 : }
1421 :
1422 : /* extend spectrum and IDFT */
1423 624 : if ( output_frame == L_FRAME48k )
1424 : {
1425 459 : ifft3( X, X, Lprot );
1426 : }
1427 : else
1428 : {
1429 165 : ifft_rel( X, Lprot, LprotLog2 );
1430 : }
1431 :
1432 624 : rec_wtda( X, ecu_rec, output_frame, Lprot, old_dec, element_mode, num_p, plocs );
1433 :
1434 624 : return;
1435 : }
1436 :
1437 :
1438 : /*--------------------------------------------------------------------------
1439 : * fir_dwn()
1440 : *
1441 : * FIR downsampling filter
1442 : *--------------------------------------------------------------------------*/
1443 :
1444 627 : 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 627 : k = 0;
1457 :
1458 : /* do the filtering */
1459 198132 : for ( i = K / 2; i < L; i += decimation )
1460 : {
1461 197505 : s = x[i] * h[0];
1462 :
1463 197505 : if ( i < K )
1464 : {
1465 3135 : mmax = i;
1466 : }
1467 : else
1468 : {
1469 194370 : mmax = K;
1470 : }
1471 :
1472 10938015 : for ( j = 1; j <= mmax; j++ )
1473 : {
1474 10740510 : s += h[j] * x[i - j];
1475 : }
1476 :
1477 197505 : y[k] = s;
1478 197505 : k++;
1479 : }
1480 :
1481 3762 : for ( ; i < L + K / 2; i += decimation )
1482 : {
1483 3135 : s = 0;
1484 :
1485 140175 : for ( j = i - L + 1; j <= K; j++ )
1486 : {
1487 137040 : s += h[j] * x[i - j];
1488 : }
1489 :
1490 3135 : y[k] = s;
1491 3135 : k++;
1492 : }
1493 :
1494 627 : return;
1495 : }
1496 :
1497 :
1498 : /*--------------------------------------------------------------------------
1499 : * fec_ecu_pitch()
1500 : *
1501 : * Pitch/correlation analysis and adaptive analysis frame length calculation
1502 : *--------------------------------------------------------------------------*/
1503 :
1504 627 : 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 627 : switch ( L )
1520 : {
1521 459 : case L_FRAME48k:
1522 459 : *decimatefator = 6;
1523 459 : filt_size = 60;
1524 459 : mvr2r( Asr_LP48, Asr_LP, filt_size + 1 );
1525 459 : break;
1526 168 : case L_FRAME32k:
1527 168 : *decimatefator = 4;
1528 168 : filt_size = 40;
1529 168 : mvr2r( Asr_LP32, Asr_LP, filt_size + 1 );
1530 168 : 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 627 : fir_dwn( prevsynth, Asr_LP, prevsynth_LP, 2 * L, filt_size, *decimatefator ); /* resampling without delay */
1547 :
1548 627 : Lon20 = (int16_t) ( ( L / 20 ) / *decimatefator );
1549 :
1550 : /* Correlation analysis */
1551 627 : *min_corr = 0;
1552 627 : accC = 0;
1553 30723 : for ( k = 0; k < 6 * Lon20; k++ )
1554 : {
1555 30096 : accC += prevsynth_LP[34 * Lon20 + k] * prevsynth_LP[34 * Lon20 + k];
1556 : }
1557 :
1558 627 : if ( HqVoicing == 1 )
1559 : {
1560 0 : cb_start = 0;
1561 0 : cb_end = 33 * Lon20;
1562 : }
1563 : else
1564 : {
1565 627 : cb_start = 0;
1566 627 : cb_end = 28 * Lon20;
1567 : }
1568 :
1569 627 : tmp_short = 34 * Lon20;
1570 627 : accB = 0;
1571 627 : delay_ind = cb_start;
1572 112716 : for ( i = cb_start; i < cb_end; i++ ) /* cb_end = 35 let 6 ms min of loop size */
1573 : {
1574 112245 : accA = 0;
1575 112245 : if ( i == cb_start )
1576 : {
1577 627 : accB = 0;
1578 30723 : for ( k = 0; k < 6 * Lon20; k++ )
1579 : {
1580 30096 : accA += prevsynth_LP[i + k] * prevsynth_LP[tmp_short + k];
1581 30096 : accB += prevsynth_LP[i + k] * prevsynth_LP[i + k];
1582 : }
1583 : }
1584 : else
1585 : {
1586 111618 : accB = accB - prevsynth_LP[i - 1] * prevsynth_LP[i - 1] + prevsynth_LP[i + 6 * Lon20 - 1] * prevsynth_LP[i + 6 * Lon20 - 1];
1587 5469282 : for ( k = 0; k < 6 * Lon20; k++ )
1588 : {
1589 5357664 : accA += prevsynth_LP[i + k] * prevsynth_LP[tmp_short + k];
1590 : }
1591 : }
1592 :
1593 : /* tests to avoid division by zero */
1594 112245 : if ( accB <= 0 )
1595 : {
1596 0 : accB = 1.0f;
1597 : }
1598 :
1599 112245 : if ( accC <= 0 )
1600 : {
1601 0 : accC = 1.0f;
1602 : }
1603 :
1604 112245 : if ( accB * accC <= 0 )
1605 : {
1606 0 : Ryy = accA;
1607 : }
1608 : else
1609 : {
1610 112245 : Ryy = accA / (float) sqrt( ( accB * accC ) );
1611 : }
1612 :
1613 112245 : if ( Ryy > *min_corr )
1614 : {
1615 4296 : *min_corr = Ryy;
1616 4296 : delay_ind = i;
1617 : }
1618 :
1619 112245 : if ( HqVoicing == 0 && *min_corr > 0.95f )
1620 : {
1621 156 : break;
1622 : }
1623 : }
1624 :
1625 627 : *N = 40 * Lon20 - delay_ind - 6 * Lon20;
1626 :
1627 627 : 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 3 : 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 3 : Lon20 = (int16_t) 160 / 20;
1654 3 : if ( element_mode == EVS_MONO )
1655 : {
1656 0 : alignment_point = 2 * 160 - 3 * Lon20;
1657 : }
1658 : else
1659 : {
1660 3 : alignment_point = 2 * 160;
1661 : }
1662 :
1663 :
1664 165 : for ( i = 0; i < N; i++ )
1665 : {
1666 162 : target[i] = prevsynth_LP[alignment_point - N + i];
1667 : }
1668 :
1669 : /* DFT */
1670 3 : *sum_Tf_abs = 0;
1671 3 : *Nfft = (int16_t) pow( 2, (int16_t) ceil( log( N ) / log( 2 ) ) );
1672 3 : tmp = ( (float) N - 1.0f ) / ( (float) *Nfft - 1.0f );
1673 :
1674 3 : set_f( Tfr, 0.0f, *Nfft );
1675 3 : set_f( Tfi, 0.0f, *Nfft );
1676 3 : Tfr[0] = target[0];
1677 3 : Tfr[*Nfft - 1] = target[N - 1];
1678 189 : for ( i = 1; i < *Nfft - 1; i++ ) /* interpolation for FFT */
1679 : {
1680 186 : tmp_short = (int16_t) floor( i * tmp );
1681 186 : Tfr[i] = target[tmp_short] + ( (float) i * tmp - ( (float) tmp_short ) ) * ( target[tmp_short + 1] - target[tmp_short] );
1682 : }
1683 :
1684 3 : DoRTFTn( Tfr, Tfi, *Nfft );
1685 3 : N_LP = (int16_t) ceil( *Nfft / 2 );
1686 :
1687 99 : for ( k = 0; k < N_LP; k++ )
1688 : {
1689 96 : Tf_abs[k] = (float) sqrt( Tfr[k] * Tfr[k] + Tfi[k] * Tfi[k] );
1690 96 : *sum_Tf_abs += Tf_abs[k];
1691 : }
1692 :
1693 3 : return;
1694 : }
1695 :
1696 : /*--------------------------------------------------------------------------*
1697 : * singenerator()
1698 : *
1699 : * fast cosinus generator Amp*cos(2*pi*freq+phi)
1700 : *--------------------------------------------------------------------------*/
1701 :
1702 :
1703 24 : 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 24 : L_C0 = a_re;
1718 24 : S0 = a_im;
1719 :
1720 24 : ptr = xx;
1721 :
1722 24 : *ptr = *ptr + L_C0;
1723 24 : ptr++;
1724 :
1725 15360 : for ( i = 0; i < L / 2 - 1; i++ )
1726 : {
1727 15336 : C0 = L_C0;
1728 15336 : L_C1 = C0 * cosfreq;
1729 15336 : L_C1 = L_C1 - S0 * sinfreq;
1730 15336 : L_S1 = C0 * sinfreq;
1731 15336 : S1 = L_S1 + S0 * cosfreq;
1732 15336 : *ptr = *ptr + L_C1;
1733 15336 : ptr++;
1734 :
1735 15336 : C1 = L_C1;
1736 15336 : L_C0 = C1 * cosfreq;
1737 15336 : L_C0 = L_C0 - S1 * sinfreq;
1738 15336 : L_S0 = C1 * sinfreq;
1739 15336 : S0 = L_S0 + S1 * cosfreq;
1740 15336 : *ptr = *ptr + L_C0;
1741 15336 : ptr++;
1742 : }
1743 :
1744 24 : C0 = L_C0;
1745 24 : L_C1 = C0 * cosfreq;
1746 24 : L_C1 = L_C1 - S0 * sinfreq;
1747 24 : *ptr = *ptr + L_C1;
1748 24 : ptr++;
1749 :
1750 24 : return;
1751 : }
1752 :
1753 :
1754 : /*--------------------------------------------------------------------------
1755 : * sinusoidal_synthesis()
1756 : *
1757 : * ECU frame sinusoid generation
1758 : *--------------------------------------------------------------------------*/
1759 :
1760 3 : 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 3 : 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 3 : 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 3 : p_pulses = pulses;
1787 3 : nb_pulses = 0;
1788 3 : new_s = Tf_abs[1];
1789 3 : glued = 1;
1790 3 : cpt = 0;
1791 3 : old = 0;
1792 :
1793 3 : PL = 0;
1794 3 : if ( N > Lon20 * 10 || HqVoicing )
1795 : {
1796 0 : PL = 1;
1797 : }
1798 78 : while ( cpt <= N / 2 - 1 - 2 )
1799 : {
1800 75 : if ( Tf_abs[cpt] > old && Tf_abs[cpt] > new_s )
1801 : {
1802 33 : glued = cpt;
1803 :
1804 66 : for ( i = glued; i < cpt + PL + 1; i++ )
1805 : {
1806 33 : *p_pulses++ = i;
1807 33 : nb_pulses++;
1808 : }
1809 33 : old = Tf_abs[cpt + PL];
1810 33 : new_s = Tf_abs[cpt + 2 + PL];
1811 33 : cpt = cpt + PL + 1;
1812 33 : glued = 1;
1813 : }
1814 : else
1815 : {
1816 42 : old = Tf_abs[cpt];
1817 42 : new_s = Tf_abs[cpt + 2];
1818 42 : cpt++;
1819 42 : glued = 0;
1820 : }
1821 : }
1822 :
1823 :
1824 3 : nb_pulses_final = 0;
1825 :
1826 : /* peak selection : keep the more energetics (max 20) */
1827 3 : tmp = 1.0f / (float) ( Nfft / 2 );
1828 3 : cumsum = 0;
1829 24 : for ( i = 0; i < min( FEC_NB_PULSE_MAX, nb_pulses ); i++ )
1830 : {
1831 24 : mmax = 0;
1832 288 : for ( k = 0; k < nb_pulses; k++ )
1833 : {
1834 264 : if ( Tf_abs[pulses[k]] > mmax )
1835 : {
1836 72 : mmax = Tf_abs[pulses[k]];
1837 72 : indmax = pulses[k];
1838 : }
1839 : }
1840 24 : cumsum += Tf_abs[indmax];
1841 :
1842 24 : if ( HqVoicing || cumsum < sum_Tf_abs * 0.7f )
1843 : {
1844 21 : a_re[i] = Tfr[indmax] * tmp;
1845 21 : a_im[i] = Tfi[indmax] * tmp;
1846 21 : freq[i] = (float) indmax * 2 / ( (float) N );
1847 21 : nb_pulses_final++;
1848 21 : Tf_abs[indmax] = -1;
1849 : }
1850 : else
1851 : {
1852 3 : a_re[i] = Tfr[indmax] * tmp;
1853 3 : a_im[i] = Tfi[indmax] * tmp;
1854 3 : freq[i] = (float) indmax * 2 / ( (float) N );
1855 3 : nb_pulses_final++;
1856 3 : break;
1857 : }
1858 : }
1859 :
1860 3 : nb_pulses = nb_pulses_final;
1861 :
1862 : /* sinusoidal synthesis */
1863 3 : set_f( synthesis, 0.0f, 40 * L / 20 );
1864 :
1865 3 : if ( HqVoicing )
1866 : {
1867 0 : k = 40 * L / 20;
1868 : }
1869 : else
1870 : {
1871 3 : k = 40 * L / 20;
1872 : }
1873 :
1874 27 : for ( i = 0; i < nb_pulses; i++ )
1875 : {
1876 24 : cosfreq = (float) cos( EVS_PI * freq[i] / (float) decimate_factor );
1877 24 : sinfreq = (float) sin( EVS_PI * freq[i] / (float) decimate_factor );
1878 24 : singenerator( k, cosfreq, sinfreq, a_re[i], a_im[i], synthesis );
1879 : }
1880 :
1881 3 : 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 3 : 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 3 : if ( element_mode == EVS_MONO )
1915 : {
1916 0 : alignment_point = 2 * L - 3 * L / 20;
1917 : }
1918 : else
1919 : {
1920 3 : alignment_point = 2 * L;
1921 : }
1922 3 : mvr2r( prevsynth + alignment_point - N, noisevect, N );
1923 :
1924 :
1925 : /* Noise addition on full band */
1926 : /* residual */
1927 3 : if ( N < L )
1928 : {
1929 3 : N_noise = ( N / 2 );
1930 651 : for ( k = 0; k < N; k++ )
1931 : {
1932 648 : 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 3 : if ( HqVoicing )
1945 : {
1946 0 : for ( i = 0; i < N; i++ )
1947 : {
1948 0 : noisevect[i] *= 0.25f;
1949 : }
1950 : }
1951 :
1952 3 : kk = 0;
1953 3 : k = 0;
1954 3 : Rnd_N_noise = N_noise;
1955 :
1956 63 : while ( k < 2 * L )
1957 : {
1958 60 : if ( kk == 0 )
1959 : {
1960 30 : tmp = ( ( own_random( ni_seed_forfec ) / (PCM16_TO_FLT_FAC) *0.2f ) + 0.5f );
1961 30 : kk = 1;
1962 : }
1963 : else
1964 : {
1965 30 : tmp = ( ( own_random( ni_seed_forfec ) / (PCM16_TO_FLT_FAC) *0.3f ) + 0.7f );
1966 30 : kk = 0;
1967 : }
1968 :
1969 60 : Rnd_N_noise = (int16_t) ( (float) N_noise * tmp );
1970 :
1971 60 : 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 3915 : for ( i = 0; i < Rnd_N_noise; i++ )
1974 : {
1975 3855 : if ( k < 2 * L )
1976 : {
1977 3840 : synthesis[k] += ( noisevect[N_noise - Rnd_N_noise + i] * SS[i] + noisevect[N_noise + i] * SS[Rnd_N_noise - i - 1] ); /* *noisefact; */
1978 : }
1979 3855 : k++;
1980 : }
1981 : }
1982 :
1983 3 : 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 3 : kk = NS2SA( L * FRAMES_PER_SEC, N_ZERO_MDCT_NS );
1991 3 : p_mdct_ola = old_out + kk;
1992 : }
1993 :
1994 : /* overlappadd with the ms of valid mdct of the last frame */
1995 3 : sinq( EVS_PI / ( 6.0f * L / 20.0f ), EVS_PI / ( 12.0f * L / 20.0f ), 3 * L / 20, SS );
1996 :
1997 291 : for ( k = 0; k < 3 * L / 20; k++ )
1998 : {
1999 288 : tmp = SS[k] * SS[k];
2000 288 : synthesis[k] = p_mdct_ola[k] * ( 1 - tmp ) + synthesis[k] * tmp;
2001 : }
2002 :
2003 3 : mvr2r( synthesis, synthesis + kk, 2 * L - kk );
2004 3 : mvr2r( synthesis + L, gapsynth, L );
2005 :
2006 3 : mvr2r( prevsynth + alignment_point - kk, synthesis, kk );
2007 :
2008 3 : return;
2009 : }
2010 :
2011 :
2012 : /*--------------------------------------------------------------------------
2013 : * fec_alg()
2014 : *
2015 : * Pitch based error-concealment algorithm with adaptive analysis frame
2016 : * length
2017 : *--------------------------------------------------------------------------*/
2018 :
2019 3 : 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 3 : fec_ecu_dft( prevsynth_LP, N, Tfr, Tfi, &sum_Tf_abs, Tf_abs, &Nfft, element_mode );
2042 :
2043 3 : sinusoidal_synthesis( Tfr, Tfi, Tf_abs, N, output_frame, decimatefactor, Nfft, sum_Tf_abs, synthesis, HqVoicing );
2044 :
2045 3 : fec_noise_filling( prevsynth, synthesis, output_frame, N * decimatefactor, HqVoicing, gapsynth, ni_seed_forfec, element_mode, old_out );
2046 :
2047 3 : n = (int16_t) ( (float) output_frame * N_ZERO_MDCT_NS / FRAME_SIZE_NS );
2048 3 : wtda( synthesis + ( output_frame - n ), ecu_rec, NULL, ALDO_WINDOW, ALDO_WINDOW, output_frame );
2049 :
2050 3 : return;
2051 : }
2052 :
2053 :
2054 : /*--------------------------------------------------------------------------
2055 : * hq_phase_ecu()
2056 : *
2057 : * Main routine for HQ phase ECU
2058 : *--------------------------------------------------------------------------*/
2059 :
2060 624 : 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 624 : noise_fac = 1.0f;
2090 :
2091 624 : 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 624 : ph_ecu_lookahead = 0;
2098 : }
2099 :
2100 :
2101 624 : Lprot = ( 2 * output_frame * 4 ) / 5;
2102 :
2103 624 : if ( !prev_bfi || ( prev_bfi && *last_fec && ( *time_offs == output_frame ) ) )
2104 : {
2105 360 : if ( !( prev_bfi && *last_fec && ( element_mode == EVS_MONO ) ) )
2106 : {
2107 360 : *time_offs = 0; /* IVAS reset of offset time counter, timeoffset variable later also used to calculate burst length */
2108 : }
2109 :
2110 360 : 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 360 : 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 360 : 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 264 : *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 264 : if ( *time_offs <= 0 )
2133 : {
2134 : /* detect wrap around of st->time_offs */
2135 30 : *time_offs = (int16_t) INT16_MAX; /* high value --> continued muting will ensure that the now saturated seed is not creating tones */
2136 : }
2137 :
2138 264 : 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 624 : mvr2r( X_sav, X, Lprot );
2142 :
2143 : /* seed for own_rand2 */
2144 624 : seed = *time_offs;
2145 624 : if ( *num_p > 0 )
2146 : {
2147 600 : 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 624 : 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 624 : old_dec = prevsynth + 2 * output_frame - NS2SA( output_frame * FRAMES_PER_SEC, N_ZERO_MDCT_NS );
2154 624 : rec_frame( X, ecu_rec, output_frame, old_dec, element_mode, num_p, plocs );
2155 :
2156 624 : return;
2157 : }
2158 :
2159 :
2160 : /*--------------------------------------------------------------------------
2161 : * hq_ecu()
2162 : *
2163 : * Main routine for HQ ECU
2164 : *--------------------------------------------------------------------------*/
2165 :
2166 627 : 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 627 : 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 627 : hHQ_core = st->hHQ_core;
2198 627 : 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 627 : fec_alg_input = prevsynth - NS2SA( output_frame * FRAMES_PER_SEC, PH_ECU_LOOKAHEAD_NS );
2205 : }
2206 :
2207 : /* find pitch and R value */
2208 627 : if ( !( output_frame < L_FRAME16k ) )
2209 : {
2210 627 : 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 3 : evs_mode_selection = ( st->total_brate >= 48000 && ( output_frame >= L_FRAME16k && !prev_bfi && ( !old_is_transient[0] || old_is_transient[1] ) &&
2220 1257 : ( ph_ecu_HqVoicing || ( ( ( hHQ_core->env_stab_plc > 0.5 ) && ( corr < 0.6 ) ) || ( hHQ_core->env_stab_plc < 0.5 && ( corr > 0.85 ) ) ) ) ) ) ||
2221 627 : ( st->total_brate < 48000 && ( ( ph_ecu_HqVoicing || corr > 0.85 ) && !prev_bfi && ( !old_is_transient[0] || old_is_transient[1] ) ) );
2222 :
2223 627 : ivas_mode_selection = ( N < PH_ECU_N_LIMIT ) || ( corr < PH_ECU_CORR_LIMIT );
2224 :
2225 627 : if ( ( ( st->element_mode == EVS_MONO ) && evs_mode_selection ) ||
2226 627 : ( ( st->element_mode != EVS_MONO ) && evs_mode_selection && ivas_mode_selection ) )
2227 :
2228 : {
2229 3 : 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 3 : *last_fec = 1;
2231 3 : *ph_ecu_active = 0;
2232 3 : *time_offs = output_frame;
2233 : }
2234 : else
2235 : {
2236 624 : 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 624 : *last_fec = 0;
2239 624 : *ph_ecu_active = 1;
2240 : }
2241 :
2242 627 : return;
2243 : }
|
