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 "prot.h"
36 : #include "ivas_prot_rend.h"
37 : #ifdef DEBUGGING
38 : #include "debug.h"
39 : #endif
40 : #include <assert.h>
41 : #include <math.h>
42 : #include <stdint.h>
43 : #include "wmc_auto.h"
44 :
45 :
46 : /*------------------------------------------------------------------------------------------*
47 : * Local constants
48 : *------------------------------------------------------------------------------------------*/
49 :
50 : #define STEP_LIMIT_PIVOT_FREQ ( 1000.0f ) /* Pivot (initial reference) frequency for response gradient limit */
51 : #define RESPONSE_STEP_LIMIT_LF ( 4.0f ) /* Maximum step in dB per bin at low frequencies (< pivot frequency) */
52 : #define RESPONSE_STEP_LIMIT_HF ( 1.5f ) /* Maximum step in dB per bin at high frequencies (> pivot frequency) */
53 : #define EPS ( 1e-30f )
54 :
55 :
56 : /*------------------------------------------------------------------------------------------*
57 : * Function definitions
58 : *------------------------------------------------------------------------------------------*/
59 :
60 : /*-------------------------------------------------------------------*
61 : * calc_min_phase()
62 : *
63 : * Compute the minimum phase spectrum that can be derived
64 : * from the amplitude spectrum of the input.
65 : *-------------------------------------------------------------------*/
66 :
67 19132 : static void calc_min_phase(
68 : rv_fftwf_type_complex *pSpectrum,
69 : const int16_t fft_size,
70 : float *pMinPhase )
71 : {
72 : int16_t idx;
73 19132 : const int16_t half_fft_size = fft_size >> 1;
74 19132 : const int16_t log2_fft_size = int_log2( fft_size );
75 19132 : const int16_t spectrum_size = half_fft_size + 1;
76 19132 : const int16_t cepstrum_smoothing_extent = fft_size >> 3;
77 19132 : const float scale_factor = 1.0f / (float) fft_size;
78 :
79 : float pCepstrum[RV_FILTER_MAX_FFT_SIZE];
80 : float pFolded_cepstrum_re[RV_FILTER_MAX_FFT_SIZE];
81 : float pFolded_cepstrum_im[RV_FILTER_MAX_FFT_SIZE];
82 : float pFolded_cepstrum_smoothing_win[RV_FILTER_MAX_FFT_SIZE];
83 : float angle_increment;
84 : float initial_angle;
85 :
86 :
87 : /* Prepare the smoothing window. Portion of the pCepstrum that will be smoothed out by a window to prevent aliasing. */
88 :
89 : /* The first portion of the pCepstrum is kept unchanged */
90 3034684 : for ( idx = 0; idx < half_fft_size - cepstrum_smoothing_extent; idx++ )
91 : {
92 3015552 : pFolded_cepstrum_smoothing_win[idx] = 1.0f;
93 : }
94 :
95 : /* Adding Hann half-window for smooth transition */
96 19132 : angle_increment = EVS_PI / (float) ( 2 * cepstrum_smoothing_extent - 1 );
97 :
98 : /* Initial angle set to match Hann window alignment in Matlab */
99 19132 : initial_angle = 0.5f * ( EVS_PI + angle_increment );
100 :
101 1024316 : for ( idx = 0; idx < cepstrum_smoothing_extent; idx++ )
102 : {
103 1005184 : const float sine_value = sinf( initial_angle + angle_increment * idx );
104 :
105 1005184 : pFolded_cepstrum_smoothing_win[idx + half_fft_size - cepstrum_smoothing_extent] = sine_value * sine_value;
106 : }
107 :
108 : /* Padding the rest with zeros */
109 4039868 : for ( idx = half_fft_size; idx < fft_size; idx++ )
110 : {
111 4020736 : pFolded_cepstrum_smoothing_win[idx] = 0.0f;
112 : }
113 :
114 : /* Compute the log amplitude spectrum */
115 4059000 : for ( idx = 0; idx < spectrum_size; idx++ )
116 : {
117 4039868 : pCepstrum[idx] = 0.5f * logf( ( pSpectrum[idx][0] * pSpectrum[idx][0] ) + ( pSpectrum[idx][1] * pSpectrum[idx][1] ) + EPS );
118 : }
119 :
120 : /* Extending one-sided spectrum to double-sided one */
121 4020736 : for ( idx = spectrum_size; idx < fft_size; idx++ )
122 : {
123 4001604 : pCepstrum[idx] = pCepstrum[fft_size - idx];
124 : }
125 :
126 : /* Compute the real pCepstrum of the log amplitude spectrum */
127 19132 : fft_rel( pCepstrum, fft_size, log2_fft_size );
128 : /* Note: using real-input fft rather than cmplx ifft to ensure support for 1024 size */
129 : /* Extra scaling is though required on the real/img parts, and img part is later inverted */
130 :
131 8060604 : for ( idx = 0; idx < fft_size; idx++ )
132 : {
133 8041472 : pCepstrum[idx] *= scale_factor;
134 : }
135 :
136 : /* Fold the pCepstrum to ensure that zeros outside the unit circle move inside it, making it minimum phase. */
137 19132 : pFolded_cepstrum_re[0] = pCepstrum[0];
138 19132 : pFolded_cepstrum_im[0] = 0.0f;
139 :
140 4020736 : for ( idx = 1; idx < half_fft_size; idx++ )
141 : {
142 4001604 : pFolded_cepstrum_re[idx] = 2.0f * pCepstrum[idx];
143 4001604 : pFolded_cepstrum_im[idx] = -2.0f * pCepstrum[fft_size - idx];
144 : /* Note: sign inverted because of fft rather than ifft used before */
145 : }
146 19132 : pFolded_cepstrum_re[half_fft_size] = pCepstrum[half_fft_size];
147 19132 : pFolded_cepstrum_im[half_fft_size] = 0.0f;
148 :
149 4020736 : for ( idx = spectrum_size; idx < fft_size; idx++ )
150 : {
151 4001604 : pFolded_cepstrum_re[idx] = 0.0f;
152 4001604 : pFolded_cepstrum_im[idx] = 0.0f;
153 : }
154 :
155 : /* Smoothing out the folded pCepstrum, to remove discontinuity at Ns,
156 : hoping to significantly reduce the secondary response in the minimum-phase IR */
157 8060604 : for ( idx = 0; idx < fft_size; idx++ )
158 : {
159 8041472 : pFolded_cepstrum_re[idx] *= pFolded_cepstrum_smoothing_win[idx];
160 8041472 : pFolded_cepstrum_im[idx] *= pFolded_cepstrum_smoothing_win[idx];
161 : }
162 :
163 : /* Convert back and isolate the phase. */
164 19132 : if ( fft_size <= 512 ) /* for size <= 512 using complex-value FFT (more effecient, but available only up to 512 size) */
165 : {
166 19132 : DoRTFTn( pFolded_cepstrum_re, pFolded_cepstrum_im, fft_size );
167 :
168 : /* Copying the img part into the output */
169 4020736 : for ( idx = 1; idx < half_fft_size; idx++ )
170 : {
171 4001604 : pMinPhase[idx] = pFolded_cepstrum_im[idx];
172 : }
173 : }
174 : else /* for fft_size > 512 using real-values FFT twice (for real + imag parts) */
175 : {
176 : /* Real part */
177 0 : fft_rel( pFolded_cepstrum_re, fft_size, log2_fft_size );
178 : /* Imag part */
179 0 : fft_rel( pFolded_cepstrum_im, fft_size, log2_fft_size );
180 :
181 : /* Copying the img part into the output */
182 0 : for ( idx = 1; idx < half_fft_size; idx++ )
183 : {
184 0 : pMinPhase[idx] = pFolded_cepstrum_re[fft_size - idx] + pFolded_cepstrum_im[idx];
185 : }
186 : }
187 19132 : pMinPhase[0] = 0.0f;
188 19132 : pMinPhase[half_fft_size] = 0.0f;
189 :
190 19132 : return;
191 : }
192 :
193 :
194 : /*-------------------------------------------------------------------*
195 : * calc_min_phase_filter()
196 : *
197 : * Convert FFT-domain pFilter pH_flt into a minimum-phase pFilter.
198 : * This function expects only the positive frequency bins up to Nyquist/2
199 : *-------------------------------------------------------------------*/
200 :
201 19132 : static void calc_min_phase_filter(
202 : rv_fftwf_type_complex *pH_flt,
203 : const int16_t fft_size,
204 : float delay )
205 : {
206 19132 : const int16_t spectrum_size = ( fft_size >> 1 ) + 1;
207 : int16_t idx;
208 : float phase_increment;
209 : float pMin_phase[RV_FILTER_MAX_FFT_SIZE];
210 :
211 19132 : phase_increment = -PI2 * delay / (float) fft_size;
212 :
213 19132 : calc_min_phase( pH_flt, fft_size, pMin_phase );
214 :
215 4059000 : for ( idx = 0; idx < spectrum_size; idx++ )
216 : {
217 : /* Cancel out initial phase by computing amplitude */
218 : /* Note: slightly different (but mathematically equivalent) approach used for better efficiency */
219 4039868 : const float current_ampl = sqrtf( pH_flt[idx][0] * pH_flt[idx][0] + pH_flt[idx][1] * pH_flt[idx][1] );
220 :
221 : /* Using the phase computed by calc_min_phase() + additional delay */
222 4039868 : const float current_phase = pMin_phase[idx] + phase_increment * idx;
223 :
224 : /* Apply the computed phase */
225 4039868 : pH_flt[idx][0] = current_ampl * cosf( current_phase );
226 4039868 : pH_flt[idx][1] = current_ampl * sinf( current_phase );
227 : }
228 :
229 19132 : return;
230 : }
231 :
232 :
233 : /*-------------------------------------------------------------------*
234 : * apply_window_fft()
235 : *
236 : * Apply the smoothing (anti-aliasing) window in the time domain
237 : *-------------------------------------------------------------------*/
238 :
239 19132 : static void apply_window_fft(
240 : rv_fftwf_type_complex *pH_flt,
241 : const float *pWindow,
242 : const int16_t fft_size )
243 : {
244 : int16_t idx;
245 : float pFilter[RV_FILTER_MAX_FFT_SIZE];
246 19132 : const int16_t half_fft_size = fft_size >> 1;
247 19132 : const int16_t log2_fft_size = int_log2( fft_size );
248 :
249 : /* Converting pFilter to ifft function input format */
250 : /* Real parts */
251 4059000 : for ( idx = 0; idx <= half_fft_size; idx++ )
252 : {
253 4039868 : pFilter[idx] = pH_flt[idx][0];
254 : }
255 : /* Img parts */
256 4020736 : for ( idx = 1; idx < half_fft_size; idx++ )
257 : {
258 4001604 : pFilter[fft_size - idx] = pH_flt[idx][1];
259 : }
260 :
261 : /* Do inverse fft to go to the time domain */
262 19132 : ifft_rel( pFilter, fft_size, log2_fft_size );
263 :
264 : /* Apply the window in the time domain */
265 8060604 : for ( idx = 0; idx < fft_size; idx++ )
266 : {
267 8041472 : pFilter[idx] *= pWindow[idx];
268 : }
269 :
270 : /* Convert back to the frequency domain */
271 19132 : fft_rel( pFilter, fft_size, log2_fft_size );
272 :
273 : /* Copy data to the output with format conversion */
274 19132 : pH_flt[0][0] = pFilter[0];
275 19132 : pH_flt[half_fft_size][0] = pFilter[half_fft_size];
276 19132 : pH_flt[0][1] = 0.0f;
277 19132 : pH_flt[half_fft_size][1] = 0.0f;
278 4020736 : for ( idx = 1; idx < half_fft_size; idx++ )
279 : {
280 4001604 : pH_flt[idx][0] = pFilter[idx];
281 4001604 : pH_flt[idx][1] = pFilter[fft_size - idx];
282 : }
283 :
284 19132 : return;
285 : }
286 :
287 :
288 : /*-------------------------------------------------------------------*
289 : * response_step_limit()
290 : *
291 : * Limit the gain vs frequency slope to T db per bin
292 : *-------------------------------------------------------------------*/
293 :
294 9566 : static void response_step_limit(
295 : float *X,
296 : const int16_t dim_x,
297 : const float step_limit_lf_dB,
298 : const float step_limit_hf_dB,
299 : const int16_t pivot_bin_idx )
300 : {
301 : int16_t i;
302 9566 : const float positive_step_limit_lf = powf( 10.0f, 0.05f * step_limit_lf_dB );
303 9566 : const float negative_step_limit_lf = 1.0f / positive_step_limit_lf;
304 9566 : const float positive_step_limit_hf = powf( 10.0f, 0.05f * step_limit_hf_dB );
305 9566 : const float negative_step_limit_hf = 1.0f / positive_step_limit_hf;
306 :
307 : /* Go up from pivot frequency and limit the slope to the maximum given by T. */
308 1882078 : for ( i = pivot_bin_idx + 1; i < dim_x; i++ )
309 : {
310 1872512 : float desiredChange = X[i] / X[i - 1];
311 1872512 : float change = fabsf( desiredChange );
312 1872512 : if ( desiredChange >= 1 )
313 : {
314 276614 : if ( change > positive_step_limit_hf )
315 : {
316 0 : change = positive_step_limit_hf;
317 : }
318 : }
319 : else
320 : {
321 1595898 : if ( change < negative_step_limit_hf )
322 : {
323 124654 : change = negative_step_limit_hf;
324 : }
325 : }
326 1872512 : X[i] = X[i - 1] * change;
327 : }
328 :
329 : /* Go down from pivot frequency and limit the slope to the maximum given by T. */
330 147422 : for ( i = pivot_bin_idx - 1; i >= 0; i-- )
331 : {
332 137856 : float desiredChange = X[i] / X[i + 1];
333 137856 : float change = fabsf( desiredChange );
334 :
335 137856 : if ( desiredChange >= 1 )
336 : {
337 51354 : if ( change > positive_step_limit_lf )
338 : {
339 8 : change = positive_step_limit_lf;
340 : }
341 : }
342 : else
343 : {
344 86502 : if ( change < negative_step_limit_lf )
345 : {
346 25456 : change = negative_step_limit_lf;
347 : }
348 : }
349 137856 : X[i] = X[i + 1] * change;
350 : }
351 :
352 9566 : return;
353 : }
354 :
355 :
356 : /*-------------------------------------------------------------------*
357 : * ivas_reverb_define_window_fft()
358 : *
359 : * Compute a smoothing window used later to avoid aliasing in FFT filters
360 : *-------------------------------------------------------------------*/
361 :
362 4783 : void ivas_reverb_define_window_fft(
363 : float *pWindow,
364 : const int16_t transitionStart,
365 : const int16_t transitionLength,
366 : const int16_t spectrumLength )
367 : {
368 : int16_t idx, fftLength;
369 : float angle_increment;
370 : float initial_angle;
371 :
372 4783 : fftLength = ( spectrumLength - 1 ) * 2;
373 :
374 : /* The first portion of the sequence is kept unchanged, window == 1 */
375 508923 : for ( idx = 0; idx < transitionStart; idx++ )
376 : {
377 504140 : pWindow[idx] = 1.0f;
378 : }
379 :
380 : /* Adding Hann half-window for smooth transition */
381 4783 : angle_increment = EVS_PI / ( 2.0f * transitionLength - 1 );
382 :
383 : /* Initial angle set to match Hann window alignment in Matlab */
384 4783 : initial_angle = 0.5f * ( EVS_PI + angle_increment );
385 193791 : for ( idx = 0; idx < transitionLength; idx++ )
386 : {
387 189008 : const float sine_value = sinf( initial_angle + angle_increment * idx );
388 :
389 189008 : pWindow[idx + transitionStart] = sine_value * sine_value;
390 : }
391 :
392 : /* Padding the rest with zeros */
393 1322003 : for ( idx = transitionStart + transitionLength; idx < fftLength; idx++ )
394 : {
395 1317220 : pWindow[idx] = 0.0f;
396 : }
397 :
398 4783 : return;
399 : }
400 :
401 : /*-------------------------------------------------------------------*
402 : * apply_window_fft()
403 : *
404 : * Applies the smoothing (anti-aliasing) window in the time domain
405 : *-------------------------------------------------------------------*/
406 :
407 : /* Computes colorations filters for the target frequency responses */
408 4783 : int16_t ivas_reverb_calc_color_filters(
409 : const float *pTargetL,
410 : const float *pTargetR,
411 : const float *pWindow,
412 : const int16_t fft_size,
413 : const float delay,
414 : rv_fftwf_type_complex *pBeqL,
415 : rv_fftwf_type_complex *pBeqR )
416 : {
417 : int16_t idx, half_fft_size;
418 :
419 4783 : half_fft_size = ( fft_size >> 1 );
420 :
421 1014750 : for ( idx = 0; idx <= half_fft_size; idx++ )
422 : {
423 1009967 : pBeqL[idx][0] = pTargetL[idx];
424 1009967 : pBeqL[idx][1] = 0.0f;
425 1009967 : pBeqR[idx][0] = pTargetR[idx];
426 1009967 : pBeqR[idx][1] = 0.0f;
427 : }
428 :
429 : /* Make the response minimum phase. Does make the spectrum complex, but
430 : will avoid aliasing when applied in the FFT domain. */
431 4783 : calc_min_phase_filter( pBeqL, fft_size, delay );
432 4783 : calc_min_phase_filter( pBeqR, fft_size, delay );
433 4783 : apply_window_fft( pBeqL, pWindow, fft_size );
434 4783 : apply_window_fft( pBeqR, pWindow, fft_size );
435 :
436 1014750 : for ( idx = 0; idx <= half_fft_size; idx++ )
437 : {
438 1009967 : if ( ( pTargetL[idx] < 0.0f ) || ( pTargetR[idx] < 0.0f ) ) /* Shouldn't have negative gains. */
439 : {
440 0 : return 1;
441 : }
442 : }
443 :
444 4783 : return 0;
445 : }
446 :
447 :
448 : /*-------------------------------------------------------------------*
449 : * ivas_reverb_calc_correl_filters()
450 : *
451 : * Compute correlation filters for the target frequency response
452 : *-------------------------------------------------------------------*/
453 :
454 4783 : int16_t ivas_reverb_calc_correl_filters(
455 : const float *pTargetICC,
456 : const float *pWindow,
457 : const int16_t fft_size,
458 : const float delay,
459 : rv_fftwf_type_complex *pU,
460 : rv_fftwf_type_complex *pV )
461 : {
462 : int16_t idx, half_fft_size;
463 :
464 4783 : half_fft_size = fft_size >> 1;
465 :
466 1014750 : for ( idx = 0; idx <= half_fft_size; idx++ )
467 : {
468 1009967 : if ( ( 1.0f + pTargetICC[idx] < 0.0f ) || ( 1.0f - pTargetICC[idx] < 0.0f ) ) /* Shouldn't have negative gains. */
469 : {
470 0 : return 1;
471 : }
472 : }
473 :
474 1014750 : for ( idx = 0; idx <= half_fft_size; idx++ )
475 : {
476 1009967 : pU[idx][0] = sqrtf( ( 1.0f + pTargetICC[idx] ) * 0.5f );
477 1009967 : pU[idx][1] = 0.0f;
478 1009967 : pV[idx][0] = sqrtf( ( 1.0f - pTargetICC[idx] ) * 0.5f );
479 1009967 : pV[idx][1] = 0.0f;
480 : }
481 :
482 : /* Make the response minimum phase. Does make the spectrum complex, but
483 : will avoid aliasing when applied in the FFT domain. */
484 4783 : calc_min_phase_filter( pU, fft_size, delay );
485 4783 : calc_min_phase_filter( pV, fft_size, delay );
486 4783 : apply_window_fft( pU, pWindow, fft_size );
487 4783 : apply_window_fft( pV, pWindow, fft_size );
488 :
489 4783 : return 0;
490 : }
491 :
492 :
493 : /*-------------------------------------------------------------------*
494 : * ivas_reverb_calc_color_levels()
495 : *
496 : * Compute the target levels (gains) for the coloration filters
497 : *-------------------------------------------------------------------*/
498 :
499 4783 : void ivas_reverb_calc_color_levels(
500 : const int32_t output_Fs,
501 : const int16_t freq_count,
502 : const int16_t loop_count,
503 : const float *pFc,
504 : const float *pAcoustic_dsr,
505 : const float *pHrtf_avg_pwr_L,
506 : const float *pHrtf_avg_pwr_R,
507 : const int16_t *pLoop_delays,
508 : const float *pT60_filter_coeff,
509 : float *pTarget_color_L,
510 : float *pTarget_color_R )
511 : {
512 : int16_t i, freq_idx, idx_pivot, nrcoefs, loop_idx;
513 : float t60[RV_LENGTH_NR_FC];
514 :
515 : /* Pre-computing inverse values for optimisation (to avoid divisions in inner loops) */
516 4783 : const float fs_inverted = 1.0f / (float) output_Fs;
517 4783 : const float loop_count_inverted = 1.0f / (float) loop_count;
518 :
519 4783 : const float log__0_001 = logf( 0.001f );
520 : /* Pre-computed for better efficiency */
521 :
522 4783 : const float revEnergyFactor = 1.0f;
523 : /* revEnergyFactor == 1 for IVAS, as the pDsr values are already scaled previously for predelay delta */
524 :
525 : /* Computing minimum delays in a simpler way, as they are already provided in descending order */
526 4783 : const int16_t minDelay = pLoop_delays[loop_count - 1];
527 4783 : const int16_t minDelayDiff = pLoop_delays[loop_count - 2] - pLoop_delays[loop_count - 1];
528 4783 : const float freq_step = 0.5f * output_Fs / (float) ( freq_count - 1 );
529 :
530 1014750 : for ( freq_idx = 0; freq_idx < freq_count; freq_idx++ )
531 : {
532 1009967 : t60[freq_idx] = 0.0f;
533 : }
534 :
535 43047 : for ( loop_idx = 0; loop_idx < loop_count; loop_idx++ )
536 : {
537 : float coefA[2];
538 : float coefB[2];
539 :
540 38264 : nrcoefs = 2;
541 : /* for pFilter order == 1 */
542 :
543 : /* Obtaining T60 filters coefficients for the current loop */
544 114792 : for ( i = 0; i < nrcoefs; i++ )
545 : {
546 76528 : coefA[i] = pT60_filter_coeff[2 * nrcoefs * loop_idx + i + nrcoefs];
547 76528 : coefB[i] = pT60_filter_coeff[2 * nrcoefs * loop_idx + i];
548 : }
549 :
550 : /* Loop over frequencies */
551 8118000 : for ( freq_idx = 0; freq_idx < freq_count; freq_idx++ )
552 : {
553 : /* Compute T60 pFilter gain at the given frequency */
554 : float cos_w;
555 : float H_filter;
556 : float T60_est;
557 :
558 8079736 : cos_w = cosf( PI2 * pFc[freq_idx] * fs_inverted );
559 8079736 : H_filter = coefB[0] * coefB[0] + coefB[1] * coefB[1] + 2.0f * coefB[0] * coefB[1] * cos_w;
560 8079736 : H_filter /= 1.0f + coefA[1] * coefA[1] + 2.0f * coefA[1] * cos_w;
561 8079736 : H_filter = sqrtf( H_filter );
562 :
563 : /* Compute the T60 value on the basis of the pFilter gain */
564 8079736 : T60_est = -3.0f * pLoop_delays[loop_idx] / ( log10f( H_filter ) * output_Fs );
565 :
566 8079736 : t60[freq_idx] += T60_est;
567 : }
568 : }
569 :
570 : /* Dividing by the number of loops to compute the average T60 estimate */
571 1014750 : for ( freq_idx = 0; freq_idx < freq_count; freq_idx++ )
572 : {
573 1009967 : t60[freq_idx] *= loop_count_inverted;
574 : }
575 :
576 : /* Compute target gains on the basis of the estimated T60 values */
577 1014750 : for ( freq_idx = 0; freq_idx < freq_count; freq_idx++ )
578 : {
579 1009967 : const float A0_square_est = powf( 10.0f, -6.0f * minDelay / ( t60[freq_idx] * output_Fs ) );
580 1009967 : const float alpha = -log__0_001 / t60[freq_idx];
581 1009967 : const float revPredNormEnergy = ( A0_square_est / ( 2.0f * alpha ) ) * output_Fs / ( 0.8776f * minDelayDiff + 26.7741f );
582 :
583 1009967 : pTarget_color_L[freq_idx] = sqrtf( pAcoustic_dsr[freq_idx] * revEnergyFactor * pHrtf_avg_pwr_L[freq_idx] ) / max( sqrtf( revPredNormEnergy ), EPS );
584 :
585 1009967 : pTarget_color_R[freq_idx] = sqrtf( pAcoustic_dsr[freq_idx] * revEnergyFactor * pHrtf_avg_pwr_R[freq_idx] ) / max( sqrtf( revPredNormEnergy ), EPS );
586 : }
587 :
588 : /* Limiting the frequency response gradients
589 : Find frequency band closest to chosen pivot frequency. */
590 4783 : idx_pivot = (int16_t) roundf( STEP_LIMIT_PIVOT_FREQ / freq_step );
591 :
592 : /* Perform step limiting */
593 4783 : response_step_limit( pTarget_color_L, freq_count, RESPONSE_STEP_LIMIT_LF, RESPONSE_STEP_LIMIT_HF, idx_pivot );
594 4783 : response_step_limit( pTarget_color_R, freq_count, RESPONSE_STEP_LIMIT_LF, RESPONSE_STEP_LIMIT_HF, idx_pivot );
595 :
596 4783 : return;
597 : }
598 :
599 :
600 : /*-------------------------------------------------------------------*
601 : * ivas_reverb_interpolate_acoustic_data()
602 : *
603 : * Interpolates data from the input T60 and DSR tables to the FFT pFilter uniform grid
604 : * Note: the fc frequencies both for the input and the output must be in the ascending order
605 : *-------------------------------------------------------------------*/
606 :
607 23597 : void ivas_reverb_interpolate_acoustic_data(
608 : const int16_t input_table_size,
609 : const float *pInput_fc,
610 : const float *pInput_t60,
611 : const float *pInput_dsr,
612 : const int16_t output_table_size,
613 : const float *pOutput_fc,
614 : float *pOutput_t60,
615 : float *pOutput_dsr )
616 : {
617 : int16_t input_idx, output_idx;
618 : float rel_offset;
619 :
620 23597 : input_idx = 0;
621 :
622 2162404 : for ( output_idx = 0; output_idx < output_table_size; output_idx++ )
623 : {
624 : /* if the bin frequency is lower than the 1st frequency point in the input table, take this 1st point */
625 2138807 : if ( pOutput_fc[output_idx] < pInput_fc[0] )
626 : {
627 4783 : input_idx = 0;
628 4783 : rel_offset = 0.0f;
629 : }
630 : else
631 : {
632 : /* if the bin frequency is higher than the last frequency point in the input table, take this last point */
633 2134024 : if ( pOutput_fc[output_idx] > pInput_fc[input_table_size - 1] )
634 : {
635 365162 : input_idx = input_table_size - 2;
636 365162 : rel_offset = 1.0f;
637 : }
638 : /* otherwise use linear interpolation between 2 consecutive points in the input table */
639 : else
640 : {
641 4113867 : while ( pOutput_fc[output_idx] > pInput_fc[input_idx + 1] )
642 : {
643 2345005 : input_idx++;
644 : }
645 1768862 : rel_offset = ( pOutput_fc[output_idx] - pInput_fc[input_idx] ) / ( pInput_fc[input_idx + 1] - pInput_fc[input_idx] );
646 : }
647 : }
648 :
649 2138807 : pOutput_t60[output_idx] = pInput_t60[input_idx] + rel_offset * ( pInput_t60[input_idx + 1] - pInput_t60[input_idx] );
650 :
651 2138807 : pOutput_dsr[output_idx] = pInput_dsr[input_idx] + rel_offset * ( pInput_dsr[input_idx + 1] - pInput_dsr[input_idx] );
652 : }
653 :
654 23597 : return;
655 : }
|
