Line data Source code
1 : /******************************************************************************
2 : * ETSI TS 103 634 V1.5.1 *
3 : * Low Complexity Communication Codec Plus (LC3plus) *
4 : * *
5 : * Copyright licence is solely granted through ETSI Intellectual Property *
6 : * Rights Policy, 3rd April 2019. No patent licence is granted by implication, *
7 : * estoppel or otherwise. *
8 : ******************************************************************************/
9 :
10 : #include "options.h"
11 : #include "wmc_auto.h"
12 : #include "defines.h"
13 : #include "functions.h"
14 : #include "constants.h"
15 :
16 : static LC3_INT32 own_rand(LC3_INT32 seed);
17 : static Complex valley_magnitude_adj(Complex X_i_in, LC3_INT32 uni_seed, LC3_FLOAT cos_F);
18 : static LC3_INT32 rand_phase(LC3_INT32 seed_in, LC3_FLOAT* cos_F);
19 :
20 : #define ONE_SIDED_SINE_WIDTH (4) /* expected pure sine main lobe maximum width (4+1+4) bins *62.5 hz/bin => approx 560 Hz total width */
21 :
22 : static LC3_INT16 plc_phEcu_nonpure_tone_ana(const LC3_INT32* plocs, const LC3_INT32 n_plocs, const Complex* X, const LC3_FLOAT* Xavg, const LC3_INT32 Lprot);
23 :
24 0 : void plc_phEcu_subst_spec(LC3_INT32* plocs, LC3_INT32 n_plocs, LC3_FLOAT* f0est, LC3_INT32 time_offs, Complex* X, LC3_INT32 X_len,
25 : LC3_FLOAT* mag_chg_gr, LC3_INT32 *seed, LC3_FLOAT* alpha, LC3_FLOAT* beta, LC3_FLOAT* Xavg,
26 : LC3_INT32 t_adv_in, LC3_INT32 Lprot, LC3_INT32 delta_corr,
27 : LC3_INT16 fadeout, /* need for DC muting */
28 : LC3_INT16* nonpure_tone_flag_ptr,
29 : LC3_FLOAT *corr_phase_dbg,
30 : LC3_FLOAT *X_i_new_re_dbg, LC3_FLOAT *X_i_new_im_dbg) {
31 :
32 : LC3_INT32 i, i2, lprotBy2Minus1, one_peak_flag_mask, noise_mag_scale;
33 : LC3_INT32 t_adv;
34 0 : LC3_FLOAT corr_phase[MAX_PLC_NPLOCS] = {0};
35 : LC3_FLOAT cos_F, mag_chg_local, alpha_local, beta_local, tmp;
36 : Complex X_i, X_i_new;
37 : LC3_INT32 segmentLen, e;
38 : LC3_FLOAT Xph;
39 : LC3_FLOAT seed_local;
40 0 : LC3_INT32 binCounter = 1, subInd = 0;
41 : LC3_INT16 fs_idx;
42 :
43 : UNUSED(corr_phase_dbg);
44 : UNUSED(X_i_new_re_dbg);
45 : UNUSED(X_i_new_im_dbg);
46 :
47 0 : seed_local = (LC3_FLOAT) *seed;
48 :
49 0 : lprotBy2Minus1 = imin(320, Lprot/2 - 1); /* limit to 20 KHz */
50 :
51 0 : t_adv = t_adv_in + time_offs;
52 :
53 0 : for (i = 0; i < n_plocs; i++) {
54 0 : corr_phase[i] = (LC3_FLOAT)2.0 * (LC3_FLOAT)M_PI_LC3PLUS * (f0est[i]/Lprot)*(LC3_FLOAT)t_adv;
55 : }
56 :
57 : // EVOLVE PHASE -----------------
58 :
59 0 : one_peak_flag_mask = -1;
60 0 : fs_idx = (LC3_INT16)LC3_FLOOR((LC3_FLOAT)Lprot / 160.0); /* aquire, fs_idx for 10 ms frame sizes */
61 0 : if (n_plocs < 3 && n_plocs > 0)
62 : {
63 0 : one_peak_flag_mask = 0; /* initial crude single tone detection, only using n_plocs as a result from peak_locator() dynamics as input */
64 :
65 0 : if ( (*nonpure_tone_flag_ptr < 0 )
66 0 : && ( (fs_idx == 2) /*SemiSWB 24 kHz */ || (fs_idx >= 4) /* FB 48 kHz */ )
67 : )
68 : {
69 : /* in the first lost frame analyze spectra to possibly reverse initial pure sine assumption */
70 0 : *nonpure_tone_flag_ptr = plc_phEcu_nonpure_tone_ana(plocs, n_plocs, X, Xavg, Lprot );
71 : }
72 :
73 0 : if ( *nonpure_tone_flag_ptr > 0 ) {
74 0 : one_peak_flag_mask = -1; /* actually revert single pure tone detection */ /* 0-> mute all surrounding valley bins in evolution , 0xff -> generate noise in all valleys */
75 : }
76 :
77 : }
78 :
79 0 : noise_mag_scale = 0;
80 0 : if (n_plocs == 0 || time_offs != 0) {
81 0 : noise_mag_scale = 1;
82 : }
83 :
84 0 : if (n_plocs == 0) {
85 0 : X[0] = realtoc(0);
86 0 : X[X_len-1] = realtoc(0);
87 : }
88 :
89 : /* binary selection of fadeout scheme */
90 : assert(PLC2_FADEOUT_LONG_IN_MS >= PLC2_FADEOUT_IN_MS_MIN && PLC2_FADEOUT_IN_MS >= PLC2_FADEOUT_IN_MS_MIN);
91 : assert(PLC2_FADEOUT_LONG_IN_MS <= PLC2_FADEOUT_IN_MS_MAX && PLC2_FADEOUT_IN_MS <= PLC2_FADEOUT_IN_MS_MAX);
92 0 : i = (PLC2_FADEOUT_IN_MS - PLC2_FADEOUT_IN_MS_MIN) / PLC2_FADEOUT_RES;
93 :
94 0 : if (fadeout != 0)
95 : {
96 0 : i = (PLC2_FADEOUT_LONG_IN_MS - PLC2_FADEOUT_IN_MS_MIN) / PLC2_FADEOUT_RES;
97 : }
98 :
99 : /* calculate local burst_len for securing DC and fs/2 muting */
100 0 : i2 = (time_offs / ((Lprot * 100) / 160)) + 1; /* burst_len */
101 :
102 0 : if (i2 > (fade_scheme_tab[i][1] + 1))
103 : {
104 : /* start DC scaling attenuation */
105 0 : X[0].r = alpha[0] * X[0].r;
106 :
107 : /* start fs/by2 attenuation */
108 0 : X[X_len - 1].r = alpha[(xavg_N_grp[fs_idx] - 1)] * X[X_len - 1].r;
109 : }
110 :
111 0 : if (n_plocs != 0) {
112 0 : for (i = 0; i < n_plocs; i++) {
113 0 : LC3_INT32 delta_corr_dn = delta_corr;
114 0 : LC3_INT32 delta_corr_up = delta_corr;
115 :
116 0 : if (i > 0) {
117 0 : delta_corr_dn = imin( ((plocs[i] - plocs[i - 1] - 1) / 2), delta_corr_dn);
118 : }
119 :
120 0 : if (i < n_plocs - 1) {
121 0 : delta_corr_up = imin( ((plocs[i + 1] - plocs[i] - 1) / 2), delta_corr_up);
122 : }
123 :
124 0 : segmentLen = (plocs[i] - delta_corr_dn) - binCounter;
125 :
126 0 : for (i2 = 0; i2 < segmentLen; i2++) {
127 0 : seed_local = (LC3_FLOAT)rand_phase((LC3_INT32)seed_local, &cos_F);
128 :
129 0 : X_i = X[binCounter];
130 0 : X_i_new = cmul(X_i, cexpi((LC3_FLOAT)M_PI_LC3PLUS*seed_local / (LC3_FLOAT)32768.0));
131 :
132 :
133 0 : seed_local = (LC3_FLOAT)own_rand((LC3_INT32)seed_local);
134 :
135 0 : if (noise_mag_scale != 0) {
136 0 : X_i = valley_magnitude_adj(X_i_new,(LC3_INT32) seed_local, cos_F);
137 0 : X_i_new = X_i;
138 : }
139 :
140 0 : mag_chg_local = mag_chg_gr[subInd];
141 0 : alpha_local = alpha[subInd];
142 :
143 0 : if (beta[subInd] != 0) {
144 0 : tmp = beta[subInd] * Xavg[subInd];
145 0 : if (one_peak_flag_mask == 0) {
146 0 : tmp = 0;
147 0 : X_i_new = realtoc(0);
148 : }
149 0 : X[binCounter] = cadd(cmul(realtoc(alpha_local), X_i_new), cmul(realtoc(tmp), cexpi((LC3_FLOAT)M_PI_LC3PLUS*seed_local / (LC3_FLOAT)32768.0)));
150 : }
151 : else {
152 0 : if (one_peak_flag_mask == 0) {
153 0 : X_i_new = realtoc(0);
154 : }
155 :
156 0 : X[binCounter] = cmul(realtoc(mag_chg_local), X_i_new);
157 : }
158 :
159 0 : binCounter++;
160 :
161 0 : if (binCounter >= gwlpr[subInd + 1]) {
162 0 : subInd++;
163 : }
164 : }
165 :
166 0 : e = plocs[i] + delta_corr_up;
167 0 : if (e > lprotBy2Minus1) {
168 0 : e = lprotBy2Minus1;
169 : }
170 :
171 0 : Xph = corr_phase[i];
172 :
173 0 : segmentLen = e - (binCounter - 1);
174 :
175 0 : for (i2 = 0; i2 < segmentLen; i2++)
176 : {
177 0 : seed_local = (LC3_FLOAT)own_rand((LC3_INT32)seed_local);
178 0 : X_i = X[binCounter];
179 :
180 : {
181 0 : LC3_INT32 nrep =(LC3_INT32) LC3_FLOOR(Xph / (2.0f*(LC3_FLOAT)M_PI_LC3PLUS));
182 :
183 0 : X_i_new = cmul(X_i, cexpi(Xph - (2.0f*(LC3_FLOAT)M_PI_LC3PLUS*(LC3_FLOAT)nrep)));
184 : }
185 :
186 :
187 0 : seed_local = (LC3_FLOAT)own_rand((LC3_INT32)seed_local);
188 :
189 0 : mag_chg_local = mag_chg_gr[subInd];
190 0 : alpha_local = alpha[subInd];
191 0 : beta_local = beta[subInd];
192 0 : if (beta_local != 0) {
193 :
194 0 : assert(alpha_local == mag_chg_local);
195 0 : tmp = beta_local * Xavg[subInd];
196 :
197 0 : X[binCounter] = cadd(cmul(realtoc(alpha_local), X_i_new), cmul(realtoc(tmp), cexpi((LC3_FLOAT)M_PI_LC3PLUS*seed_local / (LC3_FLOAT)32768.0)));
198 : }
199 : else
200 : {
201 0 : X[binCounter] = cmul(realtoc(mag_chg_local), X_i_new);
202 : }
203 :
204 0 : binCounter++;
205 :
206 0 : if (binCounter >= gwlpr[subInd + 1]) {
207 0 : subInd++;
208 : }
209 : }
210 : }
211 : }
212 :
213 0 : segmentLen = lprotBy2Minus1 - (binCounter - 1);
214 :
215 0 : for (i = 0; i < segmentLen; i++) {
216 0 : seed_local = (LC3_FLOAT)rand_phase((LC3_INT32)seed_local, &cos_F);
217 :
218 0 : X_i = X[binCounter];
219 0 : X_i_new = cmul(X_i, cexpi((LC3_FLOAT)M_PI_LC3PLUS*seed_local/(LC3_FLOAT)32768.0));
220 :
221 0 : seed_local = (LC3_FLOAT)own_rand((LC3_INT32)seed_local);
222 :
223 0 : if (noise_mag_scale != 0) {
224 0 : X_i = valley_magnitude_adj(X_i_new, (LC3_INT32)seed_local, cos_F);
225 0 : X_i_new = X_i;
226 : }
227 :
228 0 : if (one_peak_flag_mask == 0) {
229 0 : X_i_new = realtoc(0);
230 : }
231 :
232 0 : alpha_local = alpha[subInd];
233 0 : mag_chg_local = mag_chg_gr[subInd];
234 :
235 0 : if (beta[subInd] != 0) {
236 0 : assert(alpha_local == mag_chg_local);
237 0 : tmp = beta[subInd]*Xavg[subInd];
238 :
239 0 : if (one_peak_flag_mask == 0) {
240 0 : tmp = 0;
241 : }
242 :
243 0 : X[binCounter] = cadd(cmul(realtoc(alpha_local), X_i_new), cmul(realtoc(tmp), cexpi((LC3_FLOAT)M_PI_LC3PLUS*seed_local/(LC3_FLOAT)32768.0)));
244 : }
245 : else
246 : {
247 0 : X[binCounter] = cmul(realtoc(mag_chg_local), X_i_new);
248 : }
249 :
250 0 : binCounter++;
251 :
252 0 : if (binCounter >= gwlpr[subInd + 1]) {
253 0 : subInd++;
254 : }
255 : }
256 :
257 :
258 0 : *seed = (LC3_INT32)seed_local;
259 0 : }
260 :
261 0 : static LC3_INT32 own_rand(LC3_INT32 seed) {
262 : LC3_INT32 retSeed;
263 0 : assert(seed <= 32767 && seed >= -32768);
264 0 : retSeed = (13849 + (seed + 32768) * 31821) & 65535;
265 0 : retSeed -= 32768;
266 0 : assert(retSeed <= 32767 && retSeed >= -32768);
267 0 : return retSeed;
268 : }
269 :
270 0 : static Complex valley_magnitude_adj(Complex X_i_in, LC3_INT32 uni_seed, LC3_FLOAT cos_F) {
271 0 : LC3_FLOAT scale = ((LC3_FLOAT)0.5*(LC3_FLOAT)uni_seed/(LC3_FLOAT)32768.0) + (LC3_FLOAT)0.5*cos_F;
272 0 : scale = (LC3_FLOAT)1.0 + (LC3_FLOAT)0.25*scale;
273 :
274 0 : assert(scale <= (LC3_FLOAT)1.25);
275 0 : assert(scale >= (LC3_FLOAT)0.75);
276 :
277 0 : return cmul(X_i_in, realtoc(scale));
278 : }
279 :
280 0 : static LC3_INT32 rand_phase(LC3_INT32 seed_in, LC3_FLOAT* cos_F) {
281 0 : LC3_FLOAT seed = (LC3_FLOAT)own_rand(seed_in);
282 0 : *cos_F = LC3_COS((LC3_FLOAT)M_PI_LC3PLUS*seed/(LC3_FLOAT)32768.0);
283 0 : return (LC3_INT32) seed;
284 : }
285 :
286 0 : static LC3_INT16 plc_phEcu_nonpure_tone_ana(const LC3_INT32* plocs, const LC3_INT32 n_plocs, const Complex* X, const LC3_FLOAT* Xavg, const LC3_INT32 Lprot)
287 : {
288 :
289 : LC3_INT16 nonpure_tone_detect;
290 : LC3_INT16 n_ind, tone_ind, low_ind, high_ind;
291 : LC3_FLOAT peak_amp, peak_amp2, valley_amp, x_abs[(1 + 2 * ONE_SIDED_SINE_WIDTH + 2 * 1)];
292 : LC3_INT16 sineband_ind_low, sineband_ind_high;
293 : LC3_INT16 i, fs_idx, N_grp;
294 : LC3_FLOAT tmp, tmp_dB, tot_inc_HF, tot_inc_LF;
295 :
296 :
297 :
298 : /* use compressed hearing sensitivity curve to allow more deviation in highest and lowest bands */
299 : /* ROM table LC3_FLOAT scATHFx[MAX_LGW - 1] */
300 :
301 : /* init */
302 0 : nonpure_tone_detect = 0;
303 0 : tot_inc_HF = 0.0;
304 0 : tot_inc_LF = 0.0;
305 :
306 : /* limit single sine optimization to when 2 peaks are close enough to represent a single sinusoid */
307 0 : if (n_plocs == 2 && (plocs[1] - plocs[0]) >= ONE_SIDED_SINE_WIDTH) /* NB, plocs is an ordered vector */
308 : {
309 0 : nonpure_tone_detect |= 0x1;
310 : }
311 :
312 : /* local bin wise dynamics analysis, if 2 peaks, we do the analysis based on the location of the largest peak */
313 : {
314 0 : tone_ind = 0;
315 0 : plc_phEcu_fft_spec2_sqrt_approx(&(X[plocs[0]]), 1, &peak_amp); /* get 1st peak amplitude = approx_sqrt(Re^2+Im^2) */
316 :
317 :
318 0 : if ((n_plocs - 2) == 0)
319 : {
320 0 : plc_phEcu_fft_spec2_sqrt_approx(&(X[plocs[1]]), 1, &peak_amp2); /* get 2nd peak amplitude */
321 0 : if (peak_amp2 > peak_amp)
322 : {
323 0 : tone_ind = 1;
324 0 : peak_amp = peak_amp2;
325 : }
326 : }
327 :
328 0 : low_ind = MAX(1, plocs[tone_ind] - (ONE_SIDED_SINE_WIDTH + 1)); /* DC is not allowed as valley */
329 0 : high_ind = MIN((Lprot >> 1) - 2, plocs[tone_ind] + (ONE_SIDED_SINE_WIDTH + 1)); /* Fs/2 is not allowed as valley */
330 :
331 0 : n_ind = high_ind - low_ind + 1;
332 :
333 : /* find lowest amplitudes around the assumed main lobe center location */
334 0 : plc_phEcu_fft_spec2_sqrt_approx(&(X[low_ind]), n_ind, x_abs);
335 0 : valley_amp = peak_amp;
336 0 : for (i = 0; i < n_ind; i++) {
337 0 : valley_amp = MIN(x_abs[i], valley_amp);
338 : }
339 :
340 : /* at least a localized amplitude ratio of 16 (24 dB) required to declare a pure sinusoid */
341 0 : if (peak_amp < 16 * valley_amp) /* 1/16 easily implemented in BASOP */
342 : {
343 :
344 0 : nonpure_tone_detect |= 0x2;/* not a pure tone due to too low local SNR */
345 :
346 : }
347 : }
348 :
349 : /* analyze LF/ HF bands energy dynamics vs the assumed single tone band ( one or two peaks found) */
350 : {
351 0 : fs_idx = (LC3_INT16)floor(Lprot / 160); /* fs_idx */
352 0 : assert(fs_idx < 5);
353 :
354 : /* Xavg , is a vector of rather rough MDCT based band energy estimates in perceptually motivated bands. from approx the last 26 ms of synthesis */
355 :
356 : /* eval amplitude relations for assumed tonal band vs lower and higher bands */
357 0 : N_grp = xavg_N_grp[fs_idx]; /* { 4 NB , 5 WB , 6 SSWB , 7 SWB, 8 FB }; */
358 :
359 : /* establish band(s) with assumed sinusoid tone */
360 : /* if tone freq location is below first MDCT-band definition, use first band as location anyway */
361 0 : i = 0; /* band 0 , 1 , 2 , 3 , ...*/
362 0 : while (plocs[tone_ind] >= gwlpr[i + 1]) { /* gwplr= [ 1, 12(750Hz), 20(1250Hz) , 36 , .. */
363 : /* dct-inds "0"...11, 12...19, 20...35, 36 ... */
364 0 : i++;
365 : }
366 0 : sineband_ind_low = i;
367 0 : sineband_ind_high = i; /* typically in the same band as low */
368 :
369 : /* a single tone may end up on a band border
370 : , handle case when assumed tone is more or less right in between two perceptual bands +/-4 62.5 Hz */
371 0 : if ((sineband_ind_high > 0) &&
372 0 : (plocs[tone_ind] - ONE_SIDED_SINE_WIDTH) >= gwlpr[sineband_ind_high + 1]
373 : ) {
374 0 : sineband_ind_low = sineband_ind_high - 1;
375 : }
376 :
377 0 : if ( (sineband_ind_low < (N_grp - 1)) &&
378 0 : (plocs[tone_ind] + ONE_SIDED_SINE_WIDTH) >= gwlpr[sineband_ind_low + 1]
379 : ) {
380 0 : sineband_ind_high = sineband_ind_low + 1;
381 : }
382 : }
383 :
384 :
385 :
386 : /* intraframe(26 ms), weighted LB and HB envelope dynamics/variation analysis */
387 : /* envelope analysis ,
388 : require at least two HF or two LF bands in the envelope taper/roll-off analysis , otherwise skip this condition */
389 :
390 :
391 0 : if (nonpure_tone_detect == 0 &&
392 0 : (((sineband_ind_high + 2) < N_grp) ||
393 : ((sineband_ind_low - 2) >= 1)
394 : )
395 : )
396 : {
397 : /* delta taper-off analysis solution, less sensitive to input bandwidth limitation and levels */
398 :
399 : /* verify that an assumed clean sine does not have any odd HF content indications by thresholding the accumulated delta rise in LF/HF side lobes */
400 0 : for (i = (sineband_ind_high + 1); i < (N_grp - 1); i++) {
401 0 : tmp = (Xavg[i + 1] + LC3_EPS) / (Xavg[i] + LC3_EPS);
402 0 : tmp_dB = 20.0*LC3_LOGTEN(tmp);
403 0 : if ((Xavg[i] + LC3_EPS) > (Xavg[i + 1] + LC3_EPS)) {
404 0 : tmp_dB = 0;
405 : }
406 0 : tot_inc_HF += scATHFx[i] * tmp_dB; /* i is ATH factor between band i, i+1 based on Hearing sensitivity */
407 : }
408 :
409 : /* verify that an assumed clean sine does not have any odd LF content by thresholding the accumulated LF reverse up tilt */
410 0 : for (i = MAX(0, (sineband_ind_low - 1)); i > 0; i--) {
411 0 : tmp = (Xavg[i - 1] + LC3_EPS) / (Xavg[i] + LC3_EPS);
412 0 : tmp_dB = 20.0*LC3_LOGTEN(tmp); /* switch to log2() to simplify BASOP */
413 :
414 0 : if ((Xavg[i - 1] + LC3_EPS) < (Xavg[i] + LC3_EPS)) {
415 0 : tmp_dB = 0;
416 : }
417 0 : tot_inc_LF += scATHFx[i - 1] * tmp_dB; /* "psycho" scale using i-1 is ATH factor between band i-1, i , based on Hearing sensitivity */
418 : }
419 :
420 0 : if (tot_inc_HF > 4.5){ /* 4.5 dB in log2 is 0.7474 */
421 0 : nonpure_tone_detect |= 0x10; /* still not a pure tone, too great HF side increase */
422 : }
423 :
424 0 : if (tot_inc_LF > 4.5) { /* 4.5 dB limit in 4.5 = 20log10(x) corresponds to limit value 0.7474 in log2(x) */
425 0 : nonpure_tone_detect |= 0x20; /* still not a pure tone, too great accumulated LF side increase */
426 : }
427 :
428 : /* verify that an assumed clean sine does not have any odd LF+HF content by thresholding the accumulated LF+HF unexpected tilt */
429 0 : if ((tot_inc_LF + tot_inc_HF) > 6.0) { /* 6 dB limit in 20log10(x) corresponds to limit value 1.0 in log2(x) */
430 0 : nonpure_tone_detect |= 0x40; /* still not a pure tone, to great LF+HF side variation/increase */
431 : }
432 : } /* bands available*/
433 :
434 0 : return nonpure_tone_detect;
435 : }
|
