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 : #include <stdint.h>
34 : #include "options.h"
35 : #include "cnst.h"
36 : #include "prot.h"
37 : #include "ivas_cnst.h"
38 : #include "ivas_prot.h"
39 : #include "math.h"
40 : #ifdef DEBUGGING
41 : #include "debug.h"
42 : #endif
43 : #include "wmc_auto.h"
44 :
45 : /*---------------------------------------------------------------
46 : * Local constants
47 : * ---------------------------------------------------------------*/
48 :
49 : #define ZP8k 15 /* zero padding in 8kHz DFT analysis */
50 : #define OFFSET8k 55 /* offset in 8 kHz */
51 : #define STEREO_DFT_PLC_STEP21 ( L_FRAME8k - OFFSET8k ) /* Step from subframe 2 in frame n to subframe 1 in frame n+1 */
52 : #define STEREO_DFT_PLC_PH_C ( 1.0f / 3.0329f ) /* Phase estimation constant, for estimating phase of fractional frequency */
53 :
54 :
55 : /*---------------------------------------------------------------
56 : * stereo_dft_res_ecu()
57 : *
58 : * Error concealment of DFT Stereo residual, including memory
59 : * updates of DFT analysis memory and IMDCT OLA
60 : * ---------------------------------------------------------------*/
61 :
62 59088 : void stereo_dft_res_ecu(
63 : STEREO_DFT_DEC_DATA_HANDLE hStereoDft, /* i/o: Decoder DFT stereo handle */
64 : float *pDFT_RES, /* i/o: residual signal */
65 : float *const DFT_PRED_RES, /* i/o: residual prediction signal */
66 : const int16_t k, /* i : Subframe index */
67 : const int16_t output_frame, /* i : Output frame length */
68 : const int16_t prev_bfi, /* i : Previous BFI */
69 : const float dmx_nrg, /* i : Down-mix energy */
70 : int16_t *num_plocs, /* i/o: Number of peak locations */
71 : int16_t *plocs, /* i/o: Peak locations (bin) */
72 : float *plocsi, /* i/o: Peak locations (fractional) */
73 : float *input_mem /* o : Residual DFT buffer input mem */
74 : )
75 : {
76 : float res_buf[L_FRAME8k];
77 : int16_t i;
78 : int16_t L_res;
79 : float step;
80 : float fac;
81 : float trigo_dec[STEREO_DFT32MS_N_8k / 2 + 1];
82 : int16_t trigo_step;
83 : int16_t time_offs;
84 :
85 59088 : set_zero( pDFT_RES, L_FRAME8k );
86 :
87 59088 : L_res = hStereoDft->band_limits[hStereoDft->res_cod_band_max];
88 :
89 59088 : stereo_dft_res_subst_spec( hStereoDft, pDFT_RES, DFT_PRED_RES, hStereoDft->time_offs, L_res, L_FRAME8k, k, num_plocs, plocs, plocsi, k == 0 );
90 :
91 59088 : fac = (float) ( L_FRAME8k ) / (float) ( hStereoDft->NFFT );
92 :
93 59088 : if ( hStereoDft->core_hist[0] == ACELP_CORE )
94 : {
95 31470 : fac *= 0.25f;
96 : }
97 :
98 59088 : trigo_step = STEREO_DFT_TRIGO_SRATE_8k_STEP * STEREO_DFT_TRIGO_DEC_STEP;
99 2422608 : for ( i = 0; i < STEREO_DFT32MS_N_8k / 4; i++ )
100 : {
101 2363520 : trigo_dec[i] = hStereoDft->dft_trigo_8k[i * trigo_step];
102 2363520 : trigo_dec[STEREO_DFT32MS_N_8k / 2 - i] = hStereoDft->dft_trigo_8k[i * trigo_step];
103 : }
104 59088 : trigo_dec[STEREO_DFT32MS_N_8k / 4] = hStereoDft->dft_trigo_8k[STEREO_DFT32MS_N_8k / 4 * trigo_step];
105 :
106 : /* estimation of res_cod_mem (ola part in imdct residual signal) and input_mem (memory for buffer in DFT analysis)*/
107 59088 : if ( k == 0 )
108 : {
109 29544 : mvr2r( pDFT_RES, res_buf, L_FRAME8k );
110 29544 : time_offs = min( MAX16B, hStereoDft->time_offs + output_frame );
111 29544 : stereo_dft_res_subst_spec( hStereoDft, res_buf, DFT_PRED_RES, time_offs, L_res, L_FRAME8k, k, num_plocs, plocs, plocsi, FALSE );
112 :
113 29544 : rfft( res_buf, trigo_dec, L_FRAME8k, +1 );
114 :
115 29544 : v_multc( res_buf, fac, res_buf, L_FRAME8k );
116 29544 : mvr2r( res_buf + ( OFFSET8k - ZP8k ), hStereoDft->res_cod_mem, STEREO_DFT_OVL_8k );
117 :
118 29544 : mvr2r( res_buf + ZP8k, input_mem, NS2SA( 8000, STEREO_DFT32MS_OVL_NS ) ); /* Store memory for cross-fade to next frame, in case of good frame */
119 : }
120 : else
121 : {
122 29544 : mvr2r( pDFT_RES, res_buf, L_FRAME8k );
123 :
124 29544 : rfft( res_buf, trigo_dec, L_FRAME8k, +1 );
125 :
126 29544 : v_multc( res_buf, fac, res_buf, L_FRAME8k );
127 :
128 : /* Cross-fade memory */
129 29544 : fac = 0;
130 29544 : step = 1.0f / NS2SA( 8000, STEREO_DFT32MS_OVL_NS );
131 768144 : for ( i = 0; i < NS2SA( 8000, STEREO_DFT32MS_OVL_NS ); i++ )
132 : {
133 738600 : input_mem[i] = ( 1 - fac ) * res_buf[i + L_FRAME8k - NS2SA( 8000, STEREO_DFT32MS_OVL_NS ) - ZP8k] + fac * input_mem[i];
134 738600 : fac += step;
135 : }
136 :
137 : /*in case of burst error*/
138 29544 : hStereoDft->time_offs = min( MAX16B, hStereoDft->time_offs + L_FRAME8k );
139 : }
140 :
141 59088 : set_zero( DFT_PRED_RES, 2 * L_res );
142 :
143 59088 : if ( prev_bfi )
144 : {
145 18936 : stereo_dft_res_ecu_burst_att( hStereoDft, pDFT_RES, dmx_nrg, L_res, L_FRAME8k );
146 : }
147 :
148 59088 : return;
149 : }
150 :
151 :
152 : /*---------------------------------------------------------------
153 : * stereo_dft_res_subst_spec()
154 : *
155 : * Generate error concealment frame in DFT domain
156 : * ---------------------------------------------------------------*/
157 :
158 88632 : void stereo_dft_res_subst_spec(
159 : STEREO_DFT_DEC_DATA_HANDLE hStereoDft, /* i/o: Decoder DFT stereo handle */
160 : float *pDFT_RES, /* i/o: residual signal */
161 : const float *const DFT_PRED_RES, /* i : residual prediction signal */
162 : const int16_t time_offs, /* i : Time offset for phase adjustment*/
163 : const int16_t L_res, /* i : bandwidth of residual signal */
164 : const int16_t L_ana, /* i : Length of FFT analysis */
165 : const int16_t k, /* i : Subframe index */
166 : int16_t *num_plocs, /* i/o: Number of peak locations */
167 : int16_t *plocs, /* i/o: Peak locations (bin) */
168 : float *plocsi, /* i/o: Peak locations (fractional) */
169 : const int16_t analysis_flag /* i : Flag for running peak analysis */
170 : )
171 : {
172 : int16_t i, idx;
173 : float fac;
174 : float s1, s2, abs1, abs2, abs3, abs4;
175 : float abs_res[( STEREO_DFT_RES_BW_MAX ) / 2];
176 : float Xmax, Xmin;
177 : float sel;
178 : float corr_phase;
179 : float *p_mem;
180 : float f_frac;
181 : float peak_phase;
182 : float phase_tmp;
183 : float phase;
184 : float conj_sign;
185 : int16_t Np;
186 : float cos_F, sin_F;
187 :
188 : /* initialization */
189 88632 : mvr2r( DFT_PRED_RES, pDFT_RES, 2 * L_res );
190 88632 : p_mem = hStereoDft->res_mem;
191 88632 : Np = 1;
192 :
193 88632 : if ( analysis_flag )
194 : {
195 : /* Perform spectral analysis on 2nd subframe of last good frame */
196 29544 : abs_res[0] = 0.5f * ( p_mem[0] * p_mem[0] ); /* DC */
197 620424 : for ( i = 1; i < L_res; i++ )
198 : {
199 590880 : abs_res[i] = ( p_mem[2 * i] * p_mem[2 * i] + p_mem[2 * i + 1] * p_mem[2 * i + 1] );
200 : }
201 :
202 : /* Find maxima */
203 29544 : maximum( abs_res, L_res, &Xmax );
204 29544 : minimum( abs_res, L_res, &Xmin );
205 29544 : sel = ( Xmax - Xmin ) * ( 1.0f - 0.97f );
206 :
207 29544 : peakfinder( abs_res, L_res, plocs, num_plocs, sel, FALSE );
208 : /* Refine peaks */
209 113950 : for ( i = 0; i < *num_plocs; i++ )
210 : {
211 84406 : if ( plocs[i] == 0 )
212 : {
213 0 : plocsi[i] = plocs[i] + imax_pos( &abs_res[plocs[i]] );
214 : }
215 84406 : else if ( plocs[i] == L_res )
216 : {
217 0 : plocsi[i] = plocs[i] - 2 + imax_pos( &abs_res[plocs[i] - 2] );
218 : }
219 : else
220 : {
221 84406 : plocsi[i] = plocs[i] - 1 + imax_pos( &abs_res[plocs[i] - 1] );
222 : }
223 : }
224 : }
225 :
226 : /* Apply phase of stereo filling on noise spectrum */
227 1861272 : for ( i = 1; i < L_res; i++ )
228 : {
229 1772640 : s1 = sign( pDFT_RES[2 * i] );
230 1772640 : s2 = sign( pDFT_RES[2 * i + 1] );
231 1772640 : abs1 = fabsf( pDFT_RES[2 * i] );
232 1772640 : abs2 = fabsf( pDFT_RES[2 * i + 1] );
233 1772640 : abs3 = fabsf( p_mem[2 * i] );
234 1772640 : abs4 = fabsf( p_mem[2 * i + 1] );
235 :
236 1772640 : fac = 1.0f;
237 :
238 : /* Low-complex phase matching that brings the angle within pi/4 of the target angle */
239 1772640 : if ( ( ( abs1 > abs2 ) && ( abs3 < abs4 ) ) || ( ( abs1 <= abs2 ) && ( abs3 >= abs4 ) ) )
240 : {
241 899430 : pDFT_RES[2 * i] = fac * s1 * abs4;
242 899430 : pDFT_RES[2 * i + 1] = fac * s2 * abs3;
243 : }
244 : else
245 : {
246 873210 : pDFT_RES[2 * i] = fac * s1 * abs3;
247 873210 : pDFT_RES[2 * i + 1] = fac * s2 * abs4;
248 : }
249 : }
250 :
251 : /* Apply phase adjustment of identified peaks, including Np=1 peak neighbors on each side */
252 341850 : for ( i = *num_plocs - 1; i >= 0; i-- )
253 : {
254 : #ifdef NONBE_FIX_NONBE_BETWEEN_OPTIMIZATION_LEVELS_2
255 : volatile float corr_phase_tmp;
256 : #endif
257 253218 : if ( k == 0 )
258 : {
259 : /* For 1st subframe, apply reversed time ECU to get correct analysis window */
260 168812 : f_frac = plocsi[i] - plocs[i];
261 168812 : peak_phase = atan2f( p_mem[2 * plocs[i] + 1], p_mem[2 * plocs[i]] );
262 168812 : phase_tmp = peak_phase - f_frac * STEREO_DFT_PLC_PH_C;
263 168812 : phase = phase_tmp - f_frac * EVS_PI;
264 168812 : corr_phase = -2 * phase - PI2 * ( STEREO_DFT_PLC_STEP21 + L_ana + time_offs ) * ( plocsi[i] / L_ana );
265 168812 : conj_sign = -1.0f;
266 : }
267 : else
268 : {
269 : /* For 2nd subframe, do regular phase shift */
270 84406 : corr_phase = PI2 * ( L_ana + time_offs ) * ( plocsi[i] / L_ana );
271 84406 : conj_sign = 1.0f;
272 : }
273 :
274 : #ifdef NONBE_FIX_NONBE_BETWEEN_OPTIMIZATION_LEVELS_2
275 253218 : corr_phase_tmp = corr_phase;
276 253218 : cos_F = cosf( corr_phase_tmp );
277 253218 : sin_F = sinf( corr_phase_tmp );
278 : #else
279 : cos_F = cosf( corr_phase );
280 : sin_F = sinf( corr_phase );
281 : #endif
282 :
283 253218 : idx = max( 0, plocs[i] - Np ); /* Iterate over plocs[i]-1:plocs[i]+1, considering the edges of the spectrum */
284 1012872 : while ( ( idx < plocs[i] + Np + 1 ) && ( idx < L_res ) )
285 : {
286 759654 : pDFT_RES[2 * idx] = p_mem[2 * idx] * cos_F - p_mem[2 * idx + 1] * sin_F;
287 759654 : pDFT_RES[2 * idx + 1] = conj_sign * ( p_mem[2 * idx] * sin_F + p_mem[2 * idx + 1] * cos_F );
288 759654 : idx++;
289 : }
290 : }
291 :
292 88632 : return;
293 : }
294 :
295 :
296 : /*---------------------------------------------------------------
297 : * stereo_dft_res_ecu_burst_att()
298 : *
299 : * scaling residual PLC in burst error, considering DMX PLC attenuation
300 : * ---------------------------------------------------------------*/
301 :
302 18936 : void stereo_dft_res_ecu_burst_att(
303 : STEREO_DFT_DEC_DATA_HANDLE hStereoDft, /* i/o: Decoder DFT stereo handle */
304 : float *pDFT_RES, /* i/o: residual signal /att. residual */
305 : const float dmx_nrg, /* i : dmx energy of current frame */
306 : const int16_t L_res, /* i : Bandwidth of residual */
307 : const int16_t L_ana /* i : Length of FFT analysis */
308 : )
309 : {
310 : float fac;
311 :
312 : /* attenuation of residual; follow attenuation of DMX */
313 18936 : if ( hStereoDft->core_hist[0] == ACELP_CORE )
314 : {
315 8608 : fac = 0.1f * sqrtf( dmx_nrg / hStereoDft->past_dmx_nrg );
316 : }
317 : else
318 : {
319 10328 : fac = (int16_t) ( 1 - ( hStereoDft->time_offs - L_ana ) / ( hStereoDft->time_offs + L_ana ) );
320 : }
321 :
322 18936 : v_multc( pDFT_RES, fac, pDFT_RES, 2 * L_res );
323 :
324 18936 : return;
325 : }
326 :
327 :
328 : /*---------------------------------------------------------------
329 : * stereo_dft_dmx_swb_nrg()
330 : *
331 : * Calculate DMX energy
332 : * ---------------------------------------------------------------*/
333 :
334 : /*! r: total energy of downmix with maximum swb bandwidth max */
335 136164 : float stereo_dft_dmx_swb_nrg(
336 : const float *dmx_k0, /* i : first subframe spectrum */
337 : const float *dmx_k1, /* i : second subframe spectrum */
338 : const int16_t frame_length /* i : frame lanegth */
339 : )
340 : {
341 : int16_t i;
342 : float dmx_nrg;
343 :
344 136164 : dmx_nrg = EPSILON;
345 36590884 : for ( i = 0; i < frame_length / 2; i++ )
346 : {
347 36454720 : dmx_nrg += 0.5f * ( dmx_k0[2 * i] * dmx_k0[2 * i] + dmx_k0[2 * i + 1] * dmx_k0[2 * i + 1] +
348 36454720 : dmx_k1[2 * i] * dmx_k1[2 * i] + dmx_k1[2 * i + 1] * dmx_k1[2 * i + 1] );
349 : }
350 :
351 136164 : return dmx_nrg;
352 : }
353 :
354 :
355 : /*---------------------------------------------------------------
356 : * stereo_dft_sg_recovery()
357 : *
358 : * estimates panning measure
359 : * updates recovery flag that might enbale recovery of side gain
360 : * ---------------------------------------------------------------*/
361 :
362 2000493 : int16_t stereo_dft_sg_recovery(
363 : STEREO_DFT_DEC_DATA_HANDLE hStereoDft /* i/o: Decoder DFT stereo handle */
364 : )
365 : {
366 : int16_t b;
367 : float *pSideGain;
368 : float sg_m;
369 : float beta;
370 :
371 2000493 : if ( !hStereoDft->sg_mem_corrupt )
372 : {
373 1949574 : pSideGain = hStereoDft->side_gain + 2 * STEREO_DFT_BAND_MAX;
374 1949574 : beta = 0.425f;
375 :
376 1949574 : sg_m = EPSILON;
377 21367001 : for ( b = 0; b < hStereoDft->nbands; b++ )
378 : {
379 19417427 : sg_m += pSideGain[b];
380 : }
381 1949574 : sg_m /= hStereoDft->nbands;
382 :
383 1949574 : if ( sg_m < 0.6f && sg_m > -0.6f )
384 : {
385 1769372 : hStereoDft->sg_mean = 0.0f;
386 : }
387 : else
388 : {
389 180202 : hStereoDft->sg_mean = beta * sg_m + ( 1 - beta ) * hStereoDft->sg_mean; /* LP filter delta_sg to obtain side gain stability measure */
390 : }
391 : }
392 50919 : else if ( hStereoDft->sg_mean > 0.6f || hStereoDft->sg_mean < -0.6f )
393 : {
394 6732 : return 1;
395 : }
396 :
397 1993761 : return 0;
398 : }
|
