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 <string.h>
34 : #include <stdio.h>
35 : #include <stdlib.h>
36 : #include <assert.h>
37 : #include <math.h>
38 : #include "options.h"
39 : #include "cnst.h"
40 : #include "rom_enc.h"
41 : #include "rom_com.h"
42 : #include "prot.h"
43 : #include "ivas_prot.h"
44 : #include "ivas_stat_dec.h"
45 : #include "ivas_cnst.h"
46 : #include "ivas_rom_com.h"
47 : #include "ivas_rom_dec.h"
48 : #ifdef DEBUGGING
49 : #include "debug.h"
50 : #endif
51 : #include "wmc_auto.h"
52 : #include "rom_dec.h"
53 :
54 : /*-----------------------------------------------------------------------*
55 : * Local constants
56 : *-----------------------------------------------------------------------*/
57 :
58 : #define SQRT_EPSILON 3.16227755e-08 /* square root of EPSILON */
59 :
60 : /*-------------------------------------------------------------------*
61 : * ivas_dirac_dec_output_synthesis_cov_open()
62 : *
63 : * Sets up the state and parameters for the Covariance Synthesis
64 : *-------------------------------------------------------------------*/
65 :
66 1023 : ivas_error ivas_dirac_dec_output_synthesis_cov_open(
67 : DIRAC_OUTPUT_SYNTHESIS_PARAMS *h_dirac_output_synthesis_params, /* i/o: handle for the covariance synthesis parameters */
68 : DIRAC_OUTPUT_SYNTHESIS_COV_STATE *h_dirac_output_synthesis_state, /* i/o: hanlde for the covariance synthesis state */
69 : const int16_t max_band_decorr, /* i : uppermost frequency band where decorrelation is applied */
70 : const int16_t interp_length, /* i : length for interpolating the mixing matrices in time slots */
71 : const int16_t num_param_bands, /* i : number of parameter bands */
72 : const int16_t num_param_bands_residual, /* i : number of parameter bands with a residual mixing matrix (i.e. decorrelation */
73 : const int16_t nchan_in, /* i : number of input (transport) channels */
74 : const int16_t nchan_out, /* i : number of output channels */
75 : const float *proto_matrix /* i : the prototype (upmix) matrix (only used if mode == 1) */
76 : )
77 : {
78 : int16_t idx;
79 :
80 1023 : h_dirac_output_synthesis_params->max_band_decorr = max_band_decorr;
81 :
82 : /*-----------------------------------------------------------------*
83 : * memory allocation
84 : *-----------------------------------------------------------------*/
85 :
86 : /* buffer length and interpolator */
87 1023 : h_dirac_output_synthesis_params->alpha_synthesis = NULL;
88 1023 : if ( ( h_dirac_output_synthesis_params->proto_matrix = (float *) malloc( nchan_out * nchan_in * sizeof( float ) ) ) == NULL )
89 : {
90 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for DirAC synthesis covariance\n" ) );
91 : }
92 :
93 : /* cov buffers */
94 13335 : for ( idx = 0; idx < num_param_bands; idx++ )
95 : {
96 :
97 12312 : if ( ( h_dirac_output_synthesis_state->cx_old[idx] = (float *) malloc( nchan_in * nchan_in * sizeof( float ) ) ) == NULL )
98 : {
99 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for DirAC synthesis covariance\n" ) );
100 : }
101 12312 : if ( ( h_dirac_output_synthesis_state->cy_old[idx] = (float *) malloc( nchan_out * nchan_out * sizeof( float ) ) ) == NULL )
102 : {
103 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for DirAC synthesis covariance\n" ) );
104 : }
105 12312 : if ( ( h_dirac_output_synthesis_state->mixing_matrix_old[idx] = (float *) malloc( nchan_out * nchan_in * sizeof( float ) ) ) == NULL )
106 : {
107 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for DirAC synthesis covariance\n" ) );
108 : }
109 12312 : set_zero( h_dirac_output_synthesis_state->cx_old[idx], nchan_in * nchan_in );
110 12312 : set_zero( h_dirac_output_synthesis_state->cy_old[idx], nchan_out * nchan_out );
111 12312 : set_zero( h_dirac_output_synthesis_state->mixing_matrix_old[idx], nchan_out * nchan_in );
112 :
113 12312 : if ( ( h_dirac_output_synthesis_state->mixing_matrix[idx] = (float *) malloc( nchan_out * nchan_in * sizeof( float ) ) ) == NULL )
114 : {
115 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for DirAC synthesis matrix\n" ) );
116 : }
117 12312 : set_zero( h_dirac_output_synthesis_state->mixing_matrix[idx], nchan_out * nchan_in );
118 : }
119 50091 : for ( ; idx < CLDFB_NO_CHANNELS_MAX; idx++ )
120 : {
121 49068 : h_dirac_output_synthesis_state->cx_old[idx] = NULL;
122 49068 : h_dirac_output_synthesis_state->cy_old[idx] = NULL;
123 49068 : h_dirac_output_synthesis_state->mixing_matrix_old[idx] = NULL;
124 49068 : h_dirac_output_synthesis_state->mixing_matrix[idx] = NULL;
125 : }
126 :
127 10689 : for ( idx = 0; idx < num_param_bands_residual; idx++ )
128 : {
129 9666 : if ( ( h_dirac_output_synthesis_state->mixing_matrix_res_old[idx] = (float *) malloc( nchan_out * nchan_out * sizeof( float ) ) ) == NULL )
130 : {
131 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for DirAC synthesis covariance\n" ) );
132 : }
133 9666 : set_zero( h_dirac_output_synthesis_state->mixing_matrix_res_old[idx], nchan_out * nchan_out );
134 :
135 9666 : if ( ( h_dirac_output_synthesis_state->mixing_matrix_res[idx] = (float *) malloc( nchan_out * nchan_out * sizeof( float ) ) ) == NULL )
136 : {
137 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for DirAC synthesis matrix\n" ) );
138 : }
139 9666 : set_zero( h_dirac_output_synthesis_state->mixing_matrix_res[idx], nchan_out * nchan_out );
140 : }
141 52737 : for ( ; idx < CLDFB_NO_CHANNELS_MAX; idx++ )
142 : {
143 51714 : h_dirac_output_synthesis_state->mixing_matrix_res_old[idx] = NULL;
144 51714 : h_dirac_output_synthesis_state->mixing_matrix_res[idx] = NULL;
145 : }
146 :
147 : /*-----------------------------------------------------------------*
148 : * prepare processing parameters
149 : *-----------------------------------------------------------------*/
150 :
151 : /* compute interpolator */
152 1023 : if ( ( h_dirac_output_synthesis_params->interpolator = (float *) malloc( interp_length * sizeof( float ) ) ) == NULL )
153 : {
154 0 : return ( IVAS_ERROR( IVAS_ERR_FAILED_ALLOC, "Can not allocate memory for DirAC synthesis covariance\n" ) );
155 : }
156 :
157 33759 : for ( idx = 1; idx <= interp_length; ++idx )
158 : {
159 32736 : h_dirac_output_synthesis_params->interpolator[idx - 1] = (float) idx / (float) interp_length;
160 : }
161 :
162 1023 : mvr2r( proto_matrix, h_dirac_output_synthesis_params->proto_matrix, nchan_in * nchan_out );
163 :
164 1023 : return IVAS_ERR_OK;
165 : }
166 :
167 :
168 : /*-------------------------------------------------------------------*
169 : * ivas_dirac_dec_output_synthesis_get_interpolator()
170 : *
171 : *
172 : *-------------------------------------------------------------------*/
173 :
174 983961 : void ivas_dirac_dec_output_synthesis_get_interpolator(
175 : DIRAC_OUTPUT_SYNTHESIS_PARAMS *h_dirac_output_synthesis_params, /* i/o: handle for the covariance synthesis parameters */
176 : const uint16_t interp_length /* i : interpolator length */
177 : )
178 : {
179 : int16_t idx;
180 :
181 4893804 : for ( idx = 1; idx <= interp_length; ++idx )
182 : {
183 3909843 : h_dirac_output_synthesis_params->interpolator[idx - 1] = (float) idx / (float) interp_length;
184 : }
185 :
186 983961 : return;
187 : }
188 :
189 :
190 : /*-------------------------------------------------------------------*
191 : * ivas_dirac_dec_output_synthesis_cov_init()
192 : *
193 : * initialize the states for the covariance synthesis
194 : *-------------------------------------------------------------------*/
195 :
196 1023 : void ivas_dirac_dec_output_synthesis_cov_init(
197 : DIRAC_OUTPUT_SYNTHESIS_COV_STATE *h_dirac_output_synthesis_state, /* i/o: pointer to the state of the covariance synthesis */
198 : const int16_t nchan_in, /* i : number of input (tranport) channels */
199 : const int16_t nchan_out, /* i : number of output channels */
200 : const int16_t n_param_bands, /* i : number of total parameter bands */
201 : const int16_t n_param_bands_res /* i : number of parameter bands with a residual mixing matrix (i.e. decorrelation */
202 : )
203 : {
204 :
205 : int16_t idx;
206 :
207 : /* initialize buffers */
208 13335 : for ( idx = 0; idx < n_param_bands; idx++ )
209 : {
210 12312 : set_zero( h_dirac_output_synthesis_state->cx_old[idx], nchan_in * nchan_in );
211 12312 : set_zero( h_dirac_output_synthesis_state->cy_old[idx], nchan_out * nchan_out );
212 12312 : set_zero( h_dirac_output_synthesis_state->mixing_matrix_old[idx], nchan_out * nchan_in );
213 12312 : set_zero( h_dirac_output_synthesis_state->mixing_matrix[idx], nchan_out * nchan_in );
214 : }
215 :
216 10689 : for ( idx = 0; idx < n_param_bands_res; idx++ )
217 : {
218 9666 : set_zero( h_dirac_output_synthesis_state->mixing_matrix_res_old[idx], nchan_out * nchan_out );
219 9666 : set_zero( h_dirac_output_synthesis_state->mixing_matrix_res[idx], nchan_out * nchan_out );
220 : }
221 :
222 1023 : return;
223 : }
224 :
225 :
226 : /*-------------------------------------------------------------------*
227 : * ivas_dirac_dec_output_synthesis_cov_close()
228 : *
229 : * deallocate dynamic memory in the covariance synthesis state
230 : *-------------------------------------------------------------------*/
231 :
232 1023 : void ivas_dirac_dec_output_synthesis_cov_close(
233 : DIRAC_OUTPUT_SYNTHESIS_PARAMS *h_dirac_output_synthesis_params, /* i : handle for the covariance synthesis parameters */
234 : DIRAC_OUTPUT_SYNTHESIS_COV_STATE *h_dirac_output_synthesis_state /* i/o: handle for the covariance synthesis state */
235 : )
236 : {
237 : int16_t idx;
238 :
239 : /*-----------------------------------------------------------------*
240 : * memory deallocation
241 : *-----------------------------------------------------------------*/
242 :
243 : /* free interpolator */
244 1023 : if ( h_dirac_output_synthesis_params->interpolator != NULL )
245 : {
246 1023 : free( h_dirac_output_synthesis_params->interpolator );
247 1023 : h_dirac_output_synthesis_params->interpolator = NULL;
248 : }
249 :
250 : /* free alpha */
251 1023 : if ( h_dirac_output_synthesis_params->alpha_synthesis != NULL )
252 : {
253 0 : free( h_dirac_output_synthesis_params->alpha_synthesis );
254 0 : h_dirac_output_synthesis_params->alpha_synthesis = NULL;
255 : }
256 :
257 : /* free proto_matrix */
258 1023 : if ( h_dirac_output_synthesis_params->proto_matrix != NULL )
259 : {
260 1023 : free( h_dirac_output_synthesis_params->proto_matrix );
261 1023 : h_dirac_output_synthesis_params->proto_matrix = NULL;
262 : }
263 :
264 : /* free cov buffers */
265 62403 : for ( idx = 0; idx < CLDFB_NO_CHANNELS_MAX; idx++ )
266 : {
267 61380 : if ( h_dirac_output_synthesis_state->cx_old[idx] != NULL )
268 : {
269 12312 : free( h_dirac_output_synthesis_state->cx_old[idx] );
270 12312 : h_dirac_output_synthesis_state->cx_old[idx] = NULL;
271 : }
272 :
273 61380 : if ( h_dirac_output_synthesis_state->cy_old[idx] != NULL )
274 : {
275 12312 : free( h_dirac_output_synthesis_state->cy_old[idx] );
276 12312 : h_dirac_output_synthesis_state->cy_old[idx] = NULL;
277 : }
278 :
279 61380 : if ( h_dirac_output_synthesis_state->mixing_matrix_old[idx] != NULL )
280 : {
281 12312 : free( h_dirac_output_synthesis_state->mixing_matrix_old[idx] );
282 12312 : h_dirac_output_synthesis_state->mixing_matrix_old[idx] = NULL;
283 : }
284 :
285 61380 : if ( h_dirac_output_synthesis_state->mixing_matrix_res_old[idx] != NULL )
286 : {
287 9666 : free( h_dirac_output_synthesis_state->mixing_matrix_res_old[idx] );
288 9666 : h_dirac_output_synthesis_state->mixing_matrix_res_old[idx] = NULL;
289 : }
290 :
291 61380 : if ( h_dirac_output_synthesis_state->mixing_matrix[idx] != NULL )
292 : {
293 12312 : free( h_dirac_output_synthesis_state->mixing_matrix[idx] );
294 12312 : h_dirac_output_synthesis_state->mixing_matrix[idx] = NULL;
295 : }
296 :
297 61380 : if ( h_dirac_output_synthesis_state->mixing_matrix_res[idx] != NULL )
298 : {
299 9666 : free( h_dirac_output_synthesis_state->mixing_matrix_res[idx] );
300 9666 : h_dirac_output_synthesis_state->mixing_matrix_res[idx] = NULL;
301 : }
302 : }
303 :
304 1023 : return;
305 : }
306 :
307 :
308 : /*-------------------------------------------------------------------*
309 : * ivas_dirac_dec_output_synthesis_cov_param_mc_collect_slot()
310 : *
311 : * collect the multi channel input covariance for one filter bank time slot
312 : *-------------------------------------------------------------------*/
313 :
314 6255285 : void ivas_dirac_dec_output_synthesis_cov_param_mc_collect_slot(
315 : float *RealBuffer, /* i : input channel filter bank samples (real part) */
316 : float *ImagBuffer, /* i : input channel filter bank samples (imaginary part */
317 : float cx[PARAM_MC_MAX_TRANSPORT_CHANS * PARAM_MC_MAX_TRANSPORT_CHANS], /* o : accumulated input covariance (real part) */
318 : float cx_imag[PARAM_MC_MAX_TRANSPORT_CHANS * PARAM_MC_MAX_TRANSPORT_CHANS], /* o : accumulated input covariance (imaginary part) */
319 : const int16_t param_band, /* i : parameter band */
320 : PARAM_MC_DEC_HANDLE hParamMC, /* i : handle to Parametric MC state */
321 : const int16_t nchan_in /* i : number of input channels */
322 : )
323 : {
324 : int16_t band_idx, ch_idx;
325 : int16_t brange[2];
326 : float real_in_buffer[PARAM_MC_MAX_BANDS_IN_PARAMETER_BAND * MAX_TRANSPORT_CHANNELS];
327 : float imag_in_buffer[PARAM_MC_MAX_BANDS_IN_PARAMETER_BAND * MAX_TRANSPORT_CHANNELS];
328 : float real_buffer[PARAM_MC_MAX_TRANSPORT_CHANS * PARAM_MC_MAX_TRANSPORT_CHANS];
329 : float imag_buffer[PARAM_MC_MAX_TRANSPORT_CHANS * PARAM_MC_MAX_TRANSPORT_CHANS];
330 :
331 : /* estimate input covariance */
332 : /* Already stack here instead of in the process_subframe */
333 :
334 : /* collect input frame */
335 : int16_t num_bands;
336 6255285 : brange[0] = hParamMC->band_grouping[param_band];
337 6255285 : brange[1] = hParamMC->band_grouping[param_band + 1];
338 6255285 : num_bands = brange[1] - brange[0];
339 :
340 32753805 : for ( band_idx = 0; band_idx < num_bands; band_idx++ )
341 : {
342 26498520 : int16_t band = brange[0] + band_idx;
343 79701840 : for ( ch_idx = 0; ch_idx < nchan_in; ch_idx++ )
344 : {
345 53203320 : real_in_buffer[band_idx + num_bands * ch_idx] = RealBuffer[ch_idx * hParamMC->num_freq_bands + band];
346 53203320 : imag_in_buffer[band_idx + num_bands * ch_idx] = ImagBuffer[ch_idx * hParamMC->num_freq_bands + band];
347 : }
348 : }
349 :
350 6255285 : cmplx_matrix_square( real_in_buffer, imag_in_buffer, num_bands, nchan_in, real_buffer, imag_buffer );
351 :
352 6255285 : v_add( cx, real_buffer, cx, nchan_in * nchan_in );
353 :
354 6255285 : v_add( cx_imag, imag_buffer, cx_imag, nchan_in * nchan_in );
355 :
356 :
357 6255285 : return;
358 : }
359 :
360 :
361 : /*-------------------------------------------------------------------*
362 : * ivas_dirac_dec_output_synthesis_cov_param_mc_synthesise_slot()
363 : *
364 : * synthesize one filter bank slot of multi channel output filter bank
365 : * samples with the covariance synthesis
366 : *-------------------------------------------------------------------*/
367 :
368 491193 : void ivas_dirac_dec_output_synthesis_cov_param_mc_synthesise_slot(
369 : float *Cldfb_RealBuffer_in, /* i : input channel filter bank samples (real part) */
370 : float *Cldfb_ImagBuffer_in, /* i : input channel filter bank samples (imaginary part) */
371 : float Cldfb_RealBuffer[][MAX_PARAM_SPATIAL_SUBFRAMES][CLDFB_NO_CHANNELS_MAX], /* o : output channel filter bank samples (real part) */
372 : float Cldfb_ImagBuffer[][MAX_PARAM_SPATIAL_SUBFRAMES][CLDFB_NO_CHANNELS_MAX], /* o : output channel filter bank samples (imaginary part) */
373 : float *mixing_matrix[], /* i : parameter band wise mixing matrices (direct part) */
374 : float *mixing_matrix_res[], /* i : parameter band wise mixing matrices (residual part) */
375 : const uint16_t slot_idx_sfr, /* i : time slot index for the current slot within the current subframe */
376 : const uint16_t slot_idx_tot, /* i : time slot index for the current slot within the frame */
377 : const int16_t nX, /* i : number of input channels */
378 : const int16_t nY, /* i : number of output channels */
379 : PARAM_MC_DEC_HANDLE hParamMC /* i : handle to the Parametric MC decoder state */
380 : )
381 : {
382 : int16_t param_band_idx, band, ch_idx;
383 : float mixing_matrix_smooth[MAX_CICP_CHANNELS * PARAM_MC_MAX_TRANSPORT_CHANS];
384 : float mixing_matrix_res_smooth[MAX_CICP_CHANNELS * MAX_CICP_CHANNELS];
385 : float mixing_matrix_buffer[MAX_CICP_CHANNELS * MAX_CICP_CHANNELS];
386 : int16_t have_residual;
387 : float input_f_real[PARAM_MC_MAX_TRANSPORT_CHANS];
388 : float input_f_imag[PARAM_MC_MAX_TRANSPORT_CHANS];
389 : float output_f_real[MAX_CICP_CHANNELS];
390 : float output_f_imag[MAX_CICP_CHANNELS];
391 : float diff_f_real[MAX_CICP_CHANNELS];
392 : float diff_f_imag[MAX_CICP_CHANNELS];
393 : int16_t brange[2];
394 491193 : DIRAC_OUTPUT_SYNTHESIS_COV_STATE h_synthesis_state = hParamMC->h_output_synthesis_cov_state;
395 :
396 491193 : set_zero( input_f_real, PARAM_MC_MAX_TRANSPORT_CHANS );
397 491193 : set_zero( input_f_imag, PARAM_MC_MAX_TRANSPORT_CHANS );
398 491193 : set_zero( output_f_real, MAX_CICP_CHANNELS );
399 491193 : set_zero( output_f_imag, MAX_CICP_CHANNELS );
400 491193 : set_zero( diff_f_real, MAX_CICP_CHANNELS );
401 491193 : set_zero( diff_f_imag, MAX_CICP_CHANNELS );
402 :
403 6867210 : for ( param_band_idx = 0; param_band_idx < hParamMC->num_param_bands_synth; param_band_idx++ )
404 : {
405 : /* final mixing */
406 6376017 : have_residual = 0;
407 6376017 : brange[0] = hParamMC->band_grouping[param_band_idx];
408 6376017 : brange[1] = hParamMC->band_grouping[param_band_idx + 1];
409 :
410 6376017 : if ( brange[0] < hParamMC->h_output_synthesis_params.max_band_decorr )
411 : {
412 5130792 : have_residual = 1;
413 : }
414 :
415 : /* mixing matrix interpolation*/
416 6376017 : v_multc( mixing_matrix[param_band_idx], hParamMC->h_output_synthesis_params.interpolator[slot_idx_tot], mixing_matrix_smooth, nY * nX );
417 6376017 : v_multc( h_synthesis_state.mixing_matrix_old[param_band_idx], 1.0f - hParamMC->h_output_synthesis_params.interpolator[slot_idx_tot], mixing_matrix_buffer, nY * nX );
418 6376017 : v_add( mixing_matrix_smooth, mixing_matrix_buffer, mixing_matrix_smooth, nY * nX );
419 :
420 6376017 : if ( have_residual )
421 : {
422 : /* residual mixing matrix interpolation*/
423 5130792 : v_multc( mixing_matrix_res[param_band_idx], hParamMC->h_output_synthesis_params.interpolator[slot_idx_tot], mixing_matrix_res_smooth, nY * nY );
424 5130792 : v_multc( h_synthesis_state.mixing_matrix_res_old[param_band_idx], 1.0f - hParamMC->h_output_synthesis_params.interpolator[slot_idx_tot], mixing_matrix_buffer, nY * nY );
425 5130792 : v_add( mixing_matrix_res_smooth, mixing_matrix_buffer, mixing_matrix_res_smooth, nY * nY );
426 : }
427 :
428 33366777 : for ( band = brange[0]; band < brange[1]; band++ )
429 : {
430 26990760 : assert( band >= 0 );
431 :
432 26990760 : if ( have_residual )
433 : {
434 : /* collect diffuse prototypes */
435 9823860 : assert( band < hParamMC->h_output_synthesis_params.max_band_decorr );
436 69823380 : for ( ch_idx = 0; ch_idx < nY; ch_idx++ )
437 : {
438 59999520 : diff_f_real[ch_idx] = Cldfb_RealBuffer[ch_idx][slot_idx_sfr][band];
439 59999520 : diff_f_imag[ch_idx] = Cldfb_ImagBuffer[ch_idx][slot_idx_sfr][band];
440 : }
441 :
442 : /* apply residual mixing */
443 9823860 : matrix_product( mixing_matrix_res_smooth, nY, nY, 0, diff_f_real, nY, 1, 0, output_f_real );
444 9823860 : matrix_product( mixing_matrix_res_smooth, nY, nY, 0, diff_f_imag, nY, 1, 0, output_f_imag );
445 :
446 69823380 : for ( ch_idx = 0; ch_idx < nY; ch_idx++ )
447 : {
448 59999520 : Cldfb_RealBuffer[ch_idx][slot_idx_sfr][band] = output_f_real[ch_idx];
449 59999520 : Cldfb_ImagBuffer[ch_idx][slot_idx_sfr][band] = output_f_imag[ch_idx];
450 : }
451 : }
452 : else
453 : {
454 121858260 : for ( ch_idx = 0; ch_idx < nY; ch_idx++ )
455 : {
456 104691360 : Cldfb_RealBuffer[ch_idx][slot_idx_sfr][band] = 0.0f;
457 104691360 : Cldfb_ImagBuffer[ch_idx][slot_idx_sfr][band] = 0.0f;
458 : }
459 : }
460 :
461 : /* collect input signals, still in cldfb buffers */
462 81202680 : for ( ch_idx = 0; ch_idx < nX; ch_idx++ )
463 : {
464 54211920 : input_f_real[ch_idx] = Cldfb_RealBuffer_in[ch_idx * hParamMC->num_freq_bands + band];
465 54211920 : input_f_imag[ch_idx] = Cldfb_ImagBuffer_in[ch_idx * hParamMC->num_freq_bands + band];
466 : }
467 :
468 : /* apply mixing matrix */
469 26990760 : matrix_product( mixing_matrix_smooth, nY, nX, 0, input_f_real, nX, 1, 0, output_f_real );
470 26990760 : matrix_product( mixing_matrix_smooth, nY, nX, 0, input_f_imag, nX, 1, 0, output_f_imag );
471 :
472 : /* collect output */
473 191681640 : for ( ch_idx = 0; ch_idx < nY; ch_idx++ )
474 : {
475 164690880 : Cldfb_RealBuffer[ch_idx][slot_idx_sfr][band] += output_f_real[ch_idx];
476 164690880 : Cldfb_ImagBuffer[ch_idx][slot_idx_sfr][band] += output_f_imag[ch_idx];
477 : }
478 : }
479 : }
480 :
481 491193 : return;
482 : }
483 :
484 :
485 : /*-------------------------------------------------------------------*
486 : * computeMixingMatrices()
487 : *
488 : * compute a mixing matrix using the convariance synthesis approach
489 : *-------------------------------------------------------------------*/
490 :
491 398496 : int16_t computeMixingMatrices(
492 : const int16_t num_inputs, /* i : number of input channels */
493 : const int16_t num_outputs, /* i : number of output channels */
494 : const float *Cx, /* i : input channel covariance matrix */
495 : const float *Cy, /* i : target covariance matrix */
496 : const float *Q, /* i : prototype matrix (usually a upmix matrix) */
497 : const int16_t energy_compensation_flag, /* i : flag indicating that the energy compensation should be performed (i.e. no residual mixing matrix will follow) */
498 : const float reg_Sx, /* i : regularization factor for the input channel singular values */
499 : const float reg_ghat, /* i : regularization factor for the normalization matrix */
500 : float *mixing_matrix, /* o : resulting mixing matrix */
501 : float *Cr /* o : residual covariance matrix */
502 : )
503 : {
504 : int16_t i, j;
505 398496 : int16_t out = EXIT_SUCCESS;
506 : int16_t nL, nC;
507 398496 : int16_t lengthCx = num_inputs;
508 398496 : int16_t lengthCy = num_outputs;
509 : float *Cr_p, *Cy_tilde_p;
510 : const float *Cy_p;
511 : float *adj;
512 : float limit;
513 : float svd_in_buffer[MAX_OUTPUT_CHANNELS][MAX_OUTPUT_CHANNELS];
514 : float svd_u_buffer[MAX_OUTPUT_CHANNELS][MAX_OUTPUT_CHANNELS];
515 : float svd_s_buffer[MAX_OUTPUT_CHANNELS];
516 : float svd_v_buffer[MAX_OUTPUT_CHANNELS][MAX_OUTPUT_CHANNELS];
517 : float Kx[MAX_TRANSPORT_CHANNELS * MAX_TRANSPORT_CHANNELS];
518 : float Ky[MAX_OUTPUT_CHANNELS * MAX_OUTPUT_CHANNELS];
519 : float Kx_reg_inv[MAX_TRANSPORT_CHANNELS * MAX_TRANSPORT_CHANNELS];
520 : float Q_Cx[MAX_OUTPUT_CHANNELS * MAX_OUTPUT_CHANNELS];
521 : float Cy_hat_diag[MAX_OUTPUT_CHANNELS];
522 : float G_hat[MAX_OUTPUT_CHANNELS];
523 : float mat_mult_buffer1[MAX_OUTPUT_CHANNELS * MAX_OUTPUT_CHANNELS];
524 : float mat_mult_buffer2[MAX_OUTPUT_CHANNELS * MAX_OUTPUT_CHANNELS];
525 : float mat_mult_buffer3[MAX_OUTPUT_CHANNELS * MAX_OUTPUT_CHANNELS];
526 :
527 398496 : push_wmops( "dirac_cov_mix_mat" );
528 :
529 398496 : set_zero( svd_s_buffer, MAX_OUTPUT_CHANNELS );
530 6774432 : for ( i = 0; i < MAX_OUTPUT_CHANNELS; i++ )
531 : {
532 6375936 : set_zero( svd_in_buffer[i], MAX_OUTPUT_CHANNELS );
533 6375936 : set_zero( svd_u_buffer[i], MAX_OUTPUT_CHANNELS );
534 6375936 : set_zero( svd_v_buffer[i], MAX_OUTPUT_CHANNELS );
535 : }
536 :
537 : /*-----------------------------------------------------------------*
538 : * Decomposition of Cy
539 : *-----------------------------------------------------------------*/
540 :
541 : /* Processing the SVD */
542 398496 : mat2svdMat( Cy, svd_in_buffer, lengthCy, lengthCy, 0 );
543 :
544 398496 : svd( svd_in_buffer, svd_u_buffer, svd_s_buffer, svd_v_buffer, lengthCy, lengthCy );
545 :
546 : /* Computing Ky */
547 2438223 : for ( i = 0; i < lengthCy; ++i )
548 : {
549 12633024 : for ( j = 0; j < lengthCy; ++j )
550 : {
551 10593297 : Ky[i + j * lengthCy] = svd_u_buffer[i][j] * sqrtf( svd_s_buffer[j] );
552 : }
553 : }
554 :
555 : /*-----------------------------------------------------------------*
556 : * Decomposition of Cx
557 : *-----------------------------------------------------------------*/
558 :
559 : /* Processing the SVD */
560 398496 : mat2svdMat( Cx, svd_in_buffer, lengthCx, lengthCx, 0 );
561 :
562 398496 : svd( svd_in_buffer, svd_u_buffer, svd_s_buffer, svd_v_buffer, lengthCx, lengthCx );
563 :
564 : /* Computing Kx */
565 1199568 : for ( i = 0; i < lengthCx; ++i )
566 : {
567 2415456 : for ( j = 0; j < lengthCx; ++j )
568 : {
569 1614384 : Kx[i + j * lengthCx] = svd_u_buffer[i][j] * sqrtf( svd_s_buffer[j] );
570 : }
571 : }
572 :
573 : /* Sx = sqrt(Sx) */
574 1199568 : for ( i = 0; i < lengthCx; ++i )
575 : {
576 801072 : svd_s_buffer[i] = sqrtf( svd_s_buffer[i] );
577 : }
578 :
579 : /*-----------------------------------------------------------------*
580 : * Regularization of Sx
581 : *-----------------------------------------------------------------*/
582 :
583 398496 : maximum( svd_s_buffer, lengthCx, &limit );
584 398496 : limit = (float) max( limit * reg_Sx, SQRT_EPSILON );
585 1199568 : for ( i = 0; i < lengthCx; ++i )
586 : {
587 801072 : svd_s_buffer[i] = ( ( svd_s_buffer[i] > limit ) ? svd_s_buffer[i] : limit );
588 : }
589 :
590 398496 : limit = 0.0f;
591 :
592 : /*-----------------------------------------------------------------*
593 : * regularized Kx-1
594 : *-----------------------------------------------------------------*/
595 :
596 1199568 : for ( i = 0; i < lengthCx; ++i )
597 : {
598 801072 : float reg_fac = ( 1.0f / svd_s_buffer[i] );
599 :
600 2415456 : for ( j = 0; j < lengthCx; ++j )
601 : {
602 1614384 : Kx_reg_inv[i + j * lengthCx] = reg_fac * svd_u_buffer[j][i];
603 : }
604 : }
605 :
606 : /*-----------------------------------------------------------------*
607 : * normalization matrix G hat
608 : *-----------------------------------------------------------------*/
609 :
610 : /* Computing Q*Cx*Q' */
611 398496 : matrix_product( Q, lengthCy, lengthCx, 0, Cx, lengthCx, lengthCx, 0, Q_Cx );
612 :
613 398496 : matrix_product_diag( Q_Cx, lengthCy, lengthCx, 0, Q, lengthCy, lengthCx, 1, Cy_hat_diag );
614 :
615 : /* Computing Cy_hat_diag */
616 2438223 : for ( i = 0; i < lengthCy; ++i )
617 : {
618 2039727 : if ( Cy_hat_diag[i] > limit )
619 : {
620 1016814 : limit = Cy_hat_diag[i];
621 : }
622 : }
623 :
624 398496 : limit = limit * reg_ghat + EPSILON;
625 :
626 : /* Computing G_hat */
627 2438223 : for ( i = 0; i < lengthCy; ++i )
628 : {
629 2039727 : if ( limit > Cy_hat_diag[i] ) /* Computing Cy_hat_diag = max(Cy_hat_diag,limit) */
630 : {
631 54 : Cy_hat_diag[i] = limit;
632 : }
633 :
634 2039727 : G_hat[i] = sqrtf( Cy[i + i * lengthCy] / Cy_hat_diag[i] );
635 : }
636 :
637 : /*-----------------------------------------------------------------*
638 : * Formulate optimal P
639 : *-----------------------------------------------------------------*/
640 :
641 : /* Computing the input matrix Kx'*Q'*G_hat'*Ky */
642 398496 : matrix_product( Kx, lengthCx, lengthCx, 1, Q, lengthCy, lengthCx, 1, mat_mult_buffer1 );
643 398496 : matrix_diag_product( mat_mult_buffer1, lengthCx, lengthCy, 0, G_hat, lengthCy, mat_mult_buffer2 );
644 398496 : matrix_product( mat_mult_buffer2, lengthCx, lengthCy, 0, Ky, lengthCy, lengthCy, 0, mat_mult_buffer1 );
645 :
646 398496 : if ( lengthCx < lengthCy )
647 : {
648 398496 : mat2svdMat( mat_mult_buffer1, svd_in_buffer, lengthCx, lengthCy, 1 );
649 :
650 398496 : nL = lengthCy;
651 398496 : nC = lengthCx;
652 :
653 398496 : svd( svd_in_buffer, svd_v_buffer, svd_s_buffer, svd_u_buffer, nL, nC );
654 : }
655 : else
656 : {
657 0 : mat2svdMat( mat_mult_buffer1, svd_in_buffer, lengthCx, lengthCy, 0 );
658 :
659 0 : nL = lengthCx;
660 0 : nC = lengthCy;
661 :
662 0 : svd( svd_in_buffer, svd_u_buffer, svd_s_buffer, svd_v_buffer, nL, nC );
663 : }
664 :
665 : /* Actually Processing P */
666 :
667 : /* can be skipped: lambda is always column-truncated identity matrix, so this operation just
668 : truncates V to num_input_channel columns */
669 :
670 398496 : svdMat2mat( svd_v_buffer, mat_mult_buffer1, lengthCy, lengthCx );
671 398496 : svdMat2mat( svd_u_buffer, mat_mult_buffer2, lengthCx, lengthCx );
672 :
673 398496 : matrix_product( mat_mult_buffer1, lengthCy, lengthCx, 0,
674 : mat_mult_buffer2, lengthCx, lengthCx, 1,
675 : mat_mult_buffer3 );
676 :
677 : /************************ Formulate M **********************/
678 :
679 398496 : matrix_product( Ky, lengthCy, lengthCy, 0, mat_mult_buffer3, lengthCy, lengthCx, 0, mat_mult_buffer1 );
680 398496 : matrix_product( mat_mult_buffer1, lengthCy, lengthCx, 0, Kx_reg_inv, lengthCx, lengthCx, 0, mixing_matrix );
681 :
682 : /*-----------------------------------------------------------------*
683 : * Formulate Cr
684 : *-----------------------------------------------------------------*/
685 :
686 : /* Compute Cy_tilde = M*Cx*M' */
687 398496 : matrix_product( mixing_matrix, lengthCy, lengthCx, 0, Cx, lengthCx, lengthCx, 0, mat_mult_buffer1 );
688 398496 : matrix_product( mat_mult_buffer1, lengthCy, lengthCx, 0, mixing_matrix, lengthCy, lengthCx, 1, mat_mult_buffer2 );
689 :
690 : /* Compute Cr = Cy - Cy_tilde */
691 398496 : Cr_p = Cr;
692 398496 : Cy_p = Cy;
693 398496 : Cy_tilde_p = mat_mult_buffer2;
694 2438223 : for ( i = 0; i < lengthCy; ++i )
695 : {
696 12633024 : for ( j = 0; j < lengthCy; ++j )
697 : {
698 10593297 : *( Cr_p++ ) = *( Cy_p++ ) - *( Cy_tilde_p++ );
699 : }
700 :
701 : /* Avoid Meaningless negative main diagonal elements */
702 2039727 : if ( Cr[i + i * lengthCy] < 0.0f )
703 : {
704 232755 : Cr[i + i * lengthCy] = 0.0f;
705 : }
706 : }
707 :
708 : /*-----------------------------------------------------------------*
709 : * Energy Compensation
710 : *-----------------------------------------------------------------*/
711 :
712 398496 : if ( energy_compensation_flag == 1 )
713 : {
714 77826 : adj = svd_s_buffer;
715 77826 : Cy_tilde_p = mat_mult_buffer2;
716 475116 : for ( i = 0; i < lengthCy; ++i )
717 : {
718 : /* Avoid correction for very small energies,
719 : main diagonal elements of Cy_tilde_p may be negative */
720 397290 : if ( Cy_tilde_p[i + i * lengthCy] < 0.0f )
721 : {
722 0 : adj[i] = 1.0f;
723 : }
724 : else
725 : {
726 397290 : adj[i] = sqrtf( Cy[i + lengthCy * i] / ( Cy_tilde_p[i + i * lengthCy] + EPSILON ) );
727 : }
728 :
729 397290 : if ( adj[i] > 4.0f )
730 : {
731 15 : adj[i] = 4.0f;
732 : }
733 : }
734 :
735 77826 : diag_matrix_product( adj, lengthCy, mixing_matrix, lengthCy, lengthCx, 0, mat_mult_buffer3 );
736 :
737 77826 : mvr2r( mat_mult_buffer3, mixing_matrix, lengthCy * lengthCx );
738 : }
739 398496 : pop_wmops();
740 :
741 398496 : return out;
742 : }
743 :
744 :
745 : /*-------------------------------------------------------------------*
746 : * computeMixingMatricesResidual()
747 : *
748 : * compute a residual mixing matrix using the covariance synthesis approach
749 : *-------------------------------------------------------------------*/
750 :
751 320670 : int16_t computeMixingMatricesResidual(
752 : const int16_t num_outputs, /* i : number of output channels */
753 : const float *Cx, /* i : vector containing the diagonal diffuse prototype covariance */
754 : const float *Cy, /* i : matrix containing the missing cov (Cr from computeMixingMatrices()) */
755 : const float reg_Sx, /* i : regularization factor for the input channel singular values */
756 : const float reg_ghat, /* i : regularization factor for the normalization matrix */
757 : float *mixing_matrix /* o : resulting residual mixing matrix */
758 : )
759 : {
760 : int16_t i, j;
761 320670 : int16_t out = EXIT_SUCCESS;
762 320670 : int16_t lengthCx = num_outputs;
763 320670 : int16_t lengthCy = num_outputs;
764 : float *adj;
765 : float limit;
766 : float svd_in_buffer[MAX_OUTPUT_CHANNELS][MAX_OUTPUT_CHANNELS];
767 : float svd_u_buffer[MAX_OUTPUT_CHANNELS][MAX_OUTPUT_CHANNELS];
768 : float svd_s_buffer[MAX_OUTPUT_CHANNELS];
769 : float svd_v_buffer[MAX_OUTPUT_CHANNELS][MAX_OUTPUT_CHANNELS];
770 : float Kx[MAX_OUTPUT_CHANNELS];
771 : float Ky[MAX_OUTPUT_CHANNELS * MAX_OUTPUT_CHANNELS];
772 : float Kx_reg_inv[MAX_OUTPUT_CHANNELS];
773 : float Cy_hat_diag[MAX_OUTPUT_CHANNELS];
774 : float G_hat[MAX_OUTPUT_CHANNELS];
775 : float Cy_tilde[MAX_OUTPUT_CHANNELS];
776 : float mat_mult_buffer1[MAX_OUTPUT_CHANNELS * MAX_OUTPUT_CHANNELS];
777 : float mat_mult_buffer2[MAX_OUTPUT_CHANNELS * MAX_OUTPUT_CHANNELS];
778 : float mat_mult_buffer3[MAX_OUTPUT_CHANNELS * MAX_OUTPUT_CHANNELS];
779 :
780 320670 : push_wmops( "dirac_cov_mix_mat_r" );
781 :
782 : /*-----------------------------------------------------------------*
783 : * Decomposition of Cy
784 : *-----------------------------------------------------------------*/
785 :
786 : /* Processing the SVD */
787 :
788 : /* linear array to svd buffer */
789 320670 : mat2svdMat( Cy, svd_in_buffer, lengthCy, lengthCy, 0 );
790 :
791 320670 : svd( svd_in_buffer, svd_u_buffer, svd_s_buffer, svd_v_buffer, lengthCy, lengthCy );
792 :
793 : /* Computing Ky */
794 1963107 : for ( i = 0; i < lengthCy; ++i )
795 : {
796 10179204 : for ( j = 0; j < lengthCy; ++j )
797 : {
798 8536767 : Ky[i + j * lengthCy] = svd_u_buffer[i][j] * sqrtf( svd_s_buffer[j] );
799 : }
800 : }
801 :
802 : /*-----------------------------------------------------------------*
803 : * Decomposition of Cx
804 : *-----------------------------------------------------------------*/
805 :
806 : /* Processing the SVD of Cx*/
807 : /* Cx is a diagonal matrix, so SVD would lead to the sorted diagonal as S and u
808 : * would be just indicating the sorting index, so go straight to Kx as the
809 : * square root of the diagonal of Cx */
810 :
811 : /* Computing Kx */
812 1963107 : for ( i = 0; i < lengthCx; ++i )
813 : {
814 1642437 : Kx[i] = sqrtf( Cx[i] );
815 : }
816 :
817 : /*-----------------------------------------------------------------*
818 : * Regularization of Sx
819 : *-----------------------------------------------------------------*/
820 :
821 320670 : maximum( Kx, lengthCx, &limit );
822 320670 : limit = limit * reg_Sx + EPSILON;
823 :
824 1963107 : for ( i = 0; i < lengthCx; ++i )
825 : {
826 1642437 : float reg_fac = ( 1.0f / ( ( Kx[i] > limit ) ? Kx[i] : limit ) );
827 1642437 : Kx_reg_inv[i] = reg_fac;
828 : }
829 :
830 320670 : limit = 0.0f;
831 :
832 : /*-----------------------------------------------------------------*
833 : * regularized Kx-1
834 : *-----------------------------------------------------------------*/
835 :
836 : /*-----------------------------------------------------------------*
837 : * normalization matrix G hat
838 : *-----------------------------------------------------------------*/
839 :
840 : /* Computing Cy_hat_diag */
841 320670 : mvr2r( Cx, Cy_hat_diag, num_outputs );
842 :
843 1963107 : for ( i = 0; i < lengthCy; ++i )
844 : {
845 1642437 : if ( Cy_hat_diag[i] > limit )
846 : {
847 818316 : limit = Cy_hat_diag[i];
848 : }
849 : }
850 :
851 320670 : limit = limit * reg_ghat + EPSILON;
852 :
853 : /* Computing G_hat */
854 1963107 : for ( i = 0; i < lengthCy; ++i )
855 : {
856 1642437 : if ( limit > Cy_hat_diag[i] ) /* Computing Cy_hat_diag = max(Cy_hat_diag,limit) */
857 : {
858 45 : Cy_hat_diag[i] = limit;
859 : }
860 1642437 : G_hat[i] = sqrtf( Cy[i + i * lengthCy] / Cy_hat_diag[i] );
861 : }
862 :
863 : /*-----------------------------------------------------------------*
864 : * Formulate optimal P
865 : *-----------------------------------------------------------------*/
866 :
867 1963107 : for ( i = 0; i < num_outputs; i++ )
868 : {
869 1642437 : Kx[i] *= G_hat[i];
870 : }
871 :
872 1963107 : for ( i = 0; i < num_outputs; i++ )
873 : {
874 1642437 : float fac = Kx[i];
875 :
876 10179204 : for ( j = 0; j < num_outputs; j++ )
877 : {
878 8536767 : mat_mult_buffer1[i + j * num_outputs] = Ky[i + j * num_outputs] * fac;
879 : }
880 : }
881 :
882 320670 : mat2svdMat( mat_mult_buffer1, svd_in_buffer, lengthCx, lengthCy, 0 );
883 :
884 320670 : svd( svd_in_buffer, svd_u_buffer, svd_s_buffer, svd_v_buffer, lengthCx, lengthCy );
885 :
886 : /* Actually Processing P */
887 320670 : svdMat2mat( svd_v_buffer, mat_mult_buffer1, lengthCy, lengthCx );
888 320670 : svdMat2mat( svd_u_buffer, mat_mult_buffer2, lengthCx, lengthCx );
889 :
890 320670 : matrix_product( mat_mult_buffer1, lengthCy, lengthCx, 0,
891 : mat_mult_buffer2, lengthCx, lengthCx, 1,
892 : mat_mult_buffer3 );
893 :
894 : /*-----------------------------------------------------------------*
895 : * Formulate M
896 : *-----------------------------------------------------------------*/
897 :
898 320670 : matrix_product( Ky, lengthCy, lengthCy, 0, mat_mult_buffer3, lengthCy, lengthCx, 0, mat_mult_buffer1 );
899 :
900 1963107 : for ( i = 0; i < num_outputs; i++ )
901 : {
902 1642437 : float fac = Kx_reg_inv[i];
903 :
904 10179204 : for ( j = 0; j < num_outputs; j++ )
905 : {
906 8536767 : mixing_matrix[j + i * num_outputs] = mat_mult_buffer1[j + i * num_outputs] * fac;
907 : }
908 : }
909 :
910 : /*-----------------------------------------------------------------*
911 : * Formulate Cr
912 : *-----------------------------------------------------------------*/
913 :
914 : /* Compute Cy_tilde = M*Cx*M' */
915 320670 : matrix_diag_product( mixing_matrix, lengthCy, lengthCx, 0, Cx, lengthCx, mat_mult_buffer1 );
916 :
917 320670 : matrix_product_diag( mat_mult_buffer1, lengthCy, lengthCx, 0, mixing_matrix, lengthCy, lengthCx, 1, Cy_tilde );
918 :
919 : /*-----------------------------------------------------------------*
920 : * Energy Compensation
921 : *-----------------------------------------------------------------*/
922 :
923 320670 : adj = svd_s_buffer;
924 :
925 1963107 : for ( i = 0; i < lengthCy; ++i )
926 : {
927 1642437 : adj[i] = sqrtf( Cy[i + lengthCy * i] / ( Cy_tilde[i] + EPSILON ) );
928 1642437 : if ( adj[i] > 4.0f )
929 : {
930 9 : adj[i] = 4.0f;
931 : }
932 : }
933 :
934 320670 : diag_matrix_product( adj, lengthCy, mixing_matrix, lengthCy, lengthCx, 0, mat_mult_buffer3 );
935 :
936 320670 : mvr2r( mat_mult_buffer3, mixing_matrix, lengthCy * lengthCx );
937 :
938 320670 : pop_wmops();
939 :
940 320670 : return out;
941 : }
942 :
943 :
944 : /*-------------------------------------------------------------------*
945 : * computeMixingMatricesISM()
946 : *
947 : *
948 : *-------------------------------------------------------------------*/
949 :
950 1216440 : int16_t computeMixingMatricesISM(
951 : const int16_t num_inputs,
952 : const int16_t num_responses,
953 : const int16_t num_outputs,
954 : const float *responses,
955 : const float *ener,
956 : const float *Cx_diag,
957 : const float *Cy_diag,
958 : const float *Q,
959 : const int16_t energy_compensation_flag,
960 : const float reg_Sx,
961 : const float reg_ghat,
962 : float *mixing_matrix )
963 : {
964 : int16_t i, out;
965 : int16_t lengthCx, lengthCy;
966 : float *Cy_tilde_p;
967 : float *adj;
968 : float limit;
969 : float svd_in_buffer[MAX_OUTPUT_CHANNELS][MAX_OUTPUT_CHANNELS];
970 : float svd_u_buffer[MAX_OUTPUT_CHANNELS][MAX_OUTPUT_CHANNELS];
971 : float svd_s_buffer[MAX_OUTPUT_CHANNELS];
972 : float svd_v_buffer[MAX_OUTPUT_CHANNELS][MAX_OUTPUT_CHANNELS];
973 : float Kx[MAX_TRANSPORT_CHANNELS];
974 : float Ky[MAX_OUTPUT_CHANNELS * MAX_OUTPUT_CHANNELS];
975 : float Kx_reg_inv[MAX_TRANSPORT_CHANNELS];
976 : float Q_Cx[MAX_OUTPUT_CHANNELS * MAX_OUTPUT_CHANNELS];
977 : float Cy_hat_diag[MAX_OUTPUT_CHANNELS];
978 : float G_hat[MAX_OUTPUT_CHANNELS];
979 : float mat_mult_buffer1[MAX_OUTPUT_CHANNELS * MAX_OUTPUT_CHANNELS];
980 : float mat_mult_buffer2[MAX_OUTPUT_CHANNELS * MAX_OUTPUT_CHANNELS];
981 : float mat_mult_buffer3[MAX_OUTPUT_CHANNELS * MAX_OUTPUT_CHANNELS];
982 : int16_t nL, nC;
983 :
984 1216440 : push_wmops( "dirac_cov_mix_mat" );
985 :
986 1216440 : out = EXIT_SUCCESS;
987 1216440 : lengthCx = num_inputs;
988 1216440 : lengthCy = num_outputs;
989 :
990 1216440 : set_zero( svd_s_buffer, MAX_OUTPUT_CHANNELS );
991 20679480 : for ( i = 0; i < MAX_OUTPUT_CHANNELS; i++ )
992 : {
993 19463040 : set_zero( svd_in_buffer[i], MAX_OUTPUT_CHANNELS );
994 19463040 : set_zero( svd_u_buffer[i], MAX_OUTPUT_CHANNELS );
995 19463040 : set_zero( svd_v_buffer[i], MAX_OUTPUT_CHANNELS );
996 : }
997 :
998 : /* Decomposition of Cy = Ky*Ky' */
999 : /* Ky = responses*diag(ener) */
1000 1216440 : matrix_diag_product( responses, lengthCy, num_responses, 0, ener, num_responses, Ky );
1001 :
1002 : /* Decomposition of Cx -> Computing Kx */
1003 1216440 : v_sqrt( Cx_diag, Kx, lengthCx );
1004 :
1005 : /* Regularization of Sx */
1006 1216440 : maximum( Kx, lengthCx, &limit );
1007 1216440 : limit = limit * reg_Sx + EPSILON;
1008 :
1009 3649320 : for ( i = 0; i < lengthCx; ++i )
1010 : {
1011 2432880 : svd_s_buffer[i] = ( ( Kx[i] > limit ) ? Kx[i] : limit );
1012 : }
1013 :
1014 1216440 : limit = 0.0f;
1015 :
1016 : /* regularized Kx-1 */
1017 :
1018 3649320 : for ( i = 0; i < lengthCx; ++i )
1019 : {
1020 2432880 : float reg_fac = ( 1.0f / svd_s_buffer[i] );
1021 2432880 : Kx_reg_inv[i] = reg_fac;
1022 : }
1023 :
1024 : /************************ normalization matrix G hat **********************/
1025 :
1026 : /* Computing Q*Cx*Q' */
1027 1216440 : matrix_diag_product( Q, lengthCy, lengthCx, 0, Cx_diag, lengthCx, Q_Cx );
1028 1216440 : matrix_product_diag( Q_Cx, lengthCy, lengthCx, 0, Q, lengthCy, lengthCx, 1, Cy_hat_diag );
1029 :
1030 : /* Computing Cy_hat_diag */
1031 12833100 : for ( i = 0; i < lengthCy; ++i )
1032 : {
1033 11616660 : if ( Cy_hat_diag[i] > limit )
1034 : {
1035 1820064 : limit = Cy_hat_diag[i];
1036 : }
1037 : }
1038 :
1039 :
1040 1216440 : limit = limit * reg_ghat + EPSILON;
1041 :
1042 : /* Computing G_hat */
1043 12833100 : for ( i = 0; i < lengthCy; ++i )
1044 : {
1045 11616660 : if ( limit > Cy_hat_diag[i] ) /* Computing Cy_hat_diag = max(Cy_hat_diag,limit) */
1046 : {
1047 130254 : Cy_hat_diag[i] = limit;
1048 : }
1049 11616660 : G_hat[i] = sqrtf( Cy_diag[i] / Cy_hat_diag[i] );
1050 : }
1051 :
1052 : /************************ Formulate optimal P **********************/
1053 :
1054 : /* Computing the input matrix Kx'*Q'*G_hat'*Ky */
1055 1216440 : diag_matrix_product( Kx, lengthCx, Q, lengthCy, lengthCx, 1, mat_mult_buffer1 );
1056 1216440 : matrix_diag_product( mat_mult_buffer1, lengthCx, lengthCy, 0, G_hat, lengthCy, mat_mult_buffer2 );
1057 1216440 : matrix_product( mat_mult_buffer2, lengthCx, lengthCy, 0, Ky, lengthCy, num_responses, 0, mat_mult_buffer1 );
1058 :
1059 1216440 : if ( lengthCx < num_responses )
1060 : {
1061 10800 : mat2svdMat( mat_mult_buffer1, svd_in_buffer, lengthCx, num_responses, 1 );
1062 10800 : nL = num_responses;
1063 10800 : nC = lengthCx;
1064 10800 : svd( svd_in_buffer, svd_v_buffer, svd_s_buffer, svd_u_buffer, nL, nC );
1065 : }
1066 : else
1067 : {
1068 1205640 : mat2svdMat( mat_mult_buffer1, svd_in_buffer, lengthCx, num_responses, 0 );
1069 1205640 : nL = lengthCx;
1070 1205640 : nC = num_responses;
1071 1205640 : svd( svd_in_buffer, svd_u_buffer, svd_s_buffer, svd_v_buffer, nL, nC );
1072 : }
1073 :
1074 : /* Actually Processing P */
1075 :
1076 : /* can be skipped: lambda is always column-truncated identity matrix, so this operation just truncates V to num_input_channel columns */
1077 1216440 : svdMat2mat( svd_v_buffer, mat_mult_buffer1, num_responses, lengthCx );
1078 1216440 : svdMat2mat( svd_u_buffer, mat_mult_buffer2, lengthCx, lengthCx );
1079 :
1080 1216440 : matrix_product( mat_mult_buffer1, num_responses, lengthCx, 0, mat_mult_buffer2, lengthCx, lengthCx, 1, mat_mult_buffer3 );
1081 :
1082 : /************************ Formulate M **********************/
1083 :
1084 1216440 : matrix_product( Ky, lengthCy, num_responses, 0, mat_mult_buffer3, num_responses, lengthCx, 0, mat_mult_buffer1 );
1085 :
1086 1216440 : matrix_diag_product( mat_mult_buffer1, lengthCy, lengthCx, 0, Kx_reg_inv, lengthCx, mixing_matrix );
1087 :
1088 : /*********************** Energy Compensation ****************/
1089 :
1090 : /* Compute Cy_tilde = M*Cx*M' */
1091 1216440 : matrix_diag_product( mixing_matrix, lengthCy, lengthCx, 0, Cx_diag, lengthCx, mat_mult_buffer1 );
1092 1216440 : matrix_product( mat_mult_buffer1, lengthCy, lengthCx, 0, mixing_matrix, lengthCy, lengthCx, 1, mat_mult_buffer2 );
1093 :
1094 1216440 : if ( energy_compensation_flag == 1 )
1095 : {
1096 1216440 : adj = svd_s_buffer;
1097 1216440 : Cy_tilde_p = mat_mult_buffer2;
1098 12833100 : for ( i = 0; i < lengthCy; ++i )
1099 : {
1100 : /* Avoid correction for very small energies, main diagonal elements of Cy_tilde_p may be negative */
1101 11616660 : if ( Cy_tilde_p[i + i * lengthCy] < 0.0f )
1102 : {
1103 0 : adj[i] = 1.0f;
1104 : }
1105 : else
1106 : {
1107 11616660 : adj[i] = sqrtf( Cy_diag[i] / ( Cy_tilde_p[i + i * lengthCy] + EPSILON ) );
1108 : }
1109 :
1110 11616660 : if ( adj[i] > 4.0f )
1111 : {
1112 39564 : adj[i] = 4.0f;
1113 : }
1114 : }
1115 :
1116 1216440 : diag_matrix_product( adj, lengthCy, mixing_matrix, lengthCy, lengthCx, 0, mat_mult_buffer3 );
1117 :
1118 1216440 : mvr2r( mat_mult_buffer3, mixing_matrix, lengthCy * lengthCx );
1119 : }
1120 :
1121 1216440 : pop_wmops();
1122 :
1123 1216440 : return out;
1124 : }
|
