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 : #ifdef DEBUGGING
36 : #include "debug.h"
37 : #endif
38 : #include "prot.h"
39 : #include "ivas_prot.h"
40 : #include "ivas_cnst.h"
41 : #include <math.h>
42 : #include <assert.h>
43 : #include "wmc_auto.h"
44 :
45 :
46 : /*-----------------------------------------------------------------------*
47 : * Local constants
48 : *-----------------------------------------------------------------------*/
49 :
50 : #define IVAS_PCA_SM_FAC 0.0234375f
51 : #define IVAS_PCA_N_ITER_QR 10
52 :
53 :
54 : /*-----------------------------------------------------------------------*
55 : * Local function definitions
56 : *-----------------------------------------------------------------------*/
57 :
58 12 : static void ivas_bitstream_write_int32(
59 : BSTR_ENC_HANDLE hMetaData,
60 : const int32_t val,
61 : const int16_t bits )
62 : {
63 : /* MSBs */
64 12 : push_next_indice( hMetaData, (uint16_t) ( val >> 16 ), bits - 16 );
65 :
66 : /* LSBs */
67 12 : push_next_indice( hMetaData, (uint16_t) ( val & 0x0000FFFF ), 16 );
68 :
69 12 : return;
70 : }
71 :
72 :
73 4 : static void pca_enc_reset(
74 : PCA_ENC_STATE *hPCA )
75 : {
76 : int16_t i;
77 :
78 : /* reset states for interpolation and multiplexing */
79 4 : eye_matrix( hPCA->prev_eigVec, FOA_CHANNELS, 1.0f );
80 4 : set_zero( hPCA->prev_ql, IVAS_PCA_INTERP );
81 4 : hPCA->prev_ql[0] = 1.0f;
82 4 : set_zero( hPCA->prev_qr, IVAS_PCA_INTERP );
83 4 : hPCA->prev_qr[0] = 1.0f;
84 :
85 20 : for ( i = 0; i < FOA_CHANNELS; i++ )
86 : {
87 16 : hPCA->prev_D[i] = 0.25f;
88 : }
89 :
90 4 : eye_matrix( hPCA->mem_eigVec_interp, FOA_CHANNELS, 1.0f );
91 4 : set_zero( hPCA->old_r_sm, FOA_CHANNELS * FOA_CHANNELS );
92 :
93 4 : return;
94 : }
95 :
96 :
97 400 : static void pca_transform_sub(
98 : float *eigVec,
99 : float *transformed_data[8], /* i : input/transformed audio channels */
100 : const int16_t start,
101 : const int16_t len,
102 : const int16_t n_channels )
103 : {
104 : int16_t i, j, k;
105 : float temp;
106 : float buffer_data[FOA_CHANNELS];
107 10000 : for ( j = 0; j < len; j++ )
108 : {
109 48000 : for ( k = 0; k < n_channels; k++ )
110 : {
111 38400 : buffer_data[k] = transformed_data[k][j + start];
112 : }
113 48000 : for ( k = 0; k < n_channels; k++ )
114 : {
115 38400 : temp = 0.0f;
116 192000 : for ( i = 0; i < n_channels; i++ )
117 : {
118 153600 : temp += eigVec[i * IVAS_PCA_INTERP + k] * buffer_data[i];
119 : }
120 38400 : transformed_data[k][j + start] = temp;
121 : }
122 : }
123 :
124 400 : return;
125 : }
126 :
127 10 : static void pca_enc_transform(
128 : PCA_ENC_STATE *hPCA,
129 : float *ql,
130 : float *qr,
131 : float *transformed_data[8], /* i : input/transformed audio channels */
132 : const int16_t input_frame,
133 : const int16_t n_channels )
134 : {
135 : float eigVec_interp[FOA_CHANNELS * FOA_CHANNELS]; /* eigenvectors in current frame */
136 : float ql_interp[IVAS_PCA_LEN_INTERP_Q], qr_interp[IVAS_PCA_LEN_INTERP_Q];
137 : int16_t time_slot;
138 : int16_t slot_len;
139 :
140 10 : quat_shortestpath( hPCA->prev_ql, ql, hPCA->prev_qr, qr );
141 :
142 10 : pca_interp_preproc( hPCA->prev_ql, hPCA->prev_qr, ql, qr, IVAS_PCA_N_SLOTS, ql_interp, qr_interp );
143 :
144 10 : slot_len = (int16_t) ( input_frame / IVAS_PCA_N_SLOTS );
145 :
146 410 : for ( time_slot = 0; time_slot < IVAS_PCA_N_SLOTS; time_slot++ )
147 : {
148 : /* convert from double quaternion to 4D matrix */
149 400 : dquat2mat( &ql_interp[IVAS_PCA_INTERP * time_slot], &qr_interp[IVAS_PCA_INTERP * time_slot], eigVec_interp );
150 400 : pca_transform_sub( eigVec_interp, transformed_data, slot_len * time_slot, slot_len, n_channels );
151 : }
152 :
153 10 : return;
154 : }
155 :
156 :
157 1750 : static void pca_update_state(
158 : PCA_ENC_STATE *hPCA,
159 : float *ql,
160 : float *qr,
161 : float *eigVec,
162 : const int16_t n_channels )
163 : {
164 1750 : mvr2r( qr, hPCA->prev_qr, IVAS_PCA_INTERP );
165 1750 : mvr2r( ql, hPCA->prev_ql, IVAS_PCA_INTERP );
166 1750 : mvr2r( eigVec, hPCA->prev_eigVec, n_channels * n_channels );
167 :
168 1750 : return;
169 : }
170 :
171 :
172 1983 : static void swap_eigvec(
173 : float *eigVec,
174 : const int16_t i,
175 : const int16_t j )
176 : {
177 : float eigVec_tmp[FOA_CHANNELS];
178 : int16_t k;
179 :
180 9915 : for ( k = 0; k < FOA_CHANNELS; k++ )
181 : {
182 7932 : eigVec_tmp[k] = eigVec[k * FOA_CHANNELS + i];
183 : }
184 9915 : for ( k = 0; k < FOA_CHANNELS; k++ )
185 : {
186 7932 : eigVec[k * FOA_CHANNELS + i] = eigVec[k * FOA_CHANNELS + j];
187 : }
188 9915 : for ( k = 0; k < FOA_CHANNELS; k++ )
189 : {
190 7932 : eigVec[k * FOA_CHANNELS + j] = eigVec_tmp[k];
191 : }
192 :
193 1983 : return;
194 : }
195 :
196 :
197 1750 : static void sort4_D_eigVec(
198 : float *D,
199 : float *eigVec )
200 : {
201 : float tempr;
202 :
203 1750 : if ( D[0] < D[1] )
204 : {
205 0 : SWAP( D[0], D[1] );
206 0 : swap_eigvec( eigVec, 0, 1 );
207 : }
208 :
209 1750 : if ( D[2] < D[3] )
210 : {
211 60 : SWAP( D[2], D[3] );
212 60 : swap_eigvec( eigVec, 2, 3 );
213 : }
214 :
215 1750 : if ( D[0] < D[2] )
216 : {
217 0 : SWAP( D[0], D[2] );
218 0 : swap_eigvec( eigVec, 0, 2 );
219 : }
220 :
221 1750 : if ( D[1] < D[3] )
222 : {
223 0 : SWAP( D[1], D[3] );
224 0 : swap_eigvec( eigVec, 1, 3 );
225 : }
226 :
227 1750 : if ( D[1] < D[2] )
228 : {
229 26 : SWAP( D[1], D[2] );
230 26 : swap_eigvec( eigVec, 1, 2 );
231 : }
232 :
233 : /* swap last 2 values */
234 1750 : SWAP( D[2], D[3] );
235 1750 : swap_eigvec( eigVec, 2, 3 );
236 :
237 1750 : return;
238 : }
239 :
240 :
241 : /*-------------------------------------------------------------------------
242 : * ivas_pca_enc_init()
243 : *
244 : * Initialize PCA encoder
245 : *------------------------------------------------------------------------*/
246 :
247 4 : void ivas_pca_enc_init(
248 : PCA_ENC_STATE *hPCA /* i/o: PCA encoder structure */
249 : )
250 : {
251 4 : hPCA->prev_bypass_decision = PCA_MODE_INACTIVE;
252 4 : pca_enc_reset( hPCA );
253 :
254 4 : return;
255 : }
256 :
257 :
258 : /*-------------------------------------------------------------------------
259 : * ivas_pca_enc()
260 : *
261 : * PCA encoder
262 : *------------------------------------------------------------------------*/
263 :
264 1750 : void ivas_pca_enc(
265 : const ENCODER_CONFIG_HANDLE hEncoderConfig, /* i : configuration structure */
266 : PCA_ENC_STATE *hPCA, /* i : PCA encoder structure */
267 : BSTR_ENC_HANDLE hMetaData, /* i/o: MetaData handle */
268 : float *data_f[8], /* i : input/transformed audio channels */
269 : const int16_t input_frame, /* i : input frame length */
270 : const int16_t n_channels /* i : number of channels */
271 : )
272 : {
273 : int16_t i, j, k, l;
274 : int16_t len_subfr;
275 : float r[FOA_CHANNELS * FOA_CHANNELS]; /* covariance matrix */
276 : float eigVec[FOA_CHANNELS * FOA_CHANNELS]; /* eigenvectors in current frame */
277 : float eigVec_tmp[FOA_CHANNELS * FOA_CHANNELS]; /* transpose / tmp matrix for permutation */
278 : float cost_mtx[FOA_CHANNELS * FOA_CHANNELS]; /* cost matrix */
279 : int16_t path[FOA_CHANNELS];
280 : float D[FOA_CHANNELS]; /* eigenvalues in current frame */
281 : float ql[IVAS_PCA_INTERP], qr[IVAS_PCA_INTERP];
282 : float det;
283 : float temp;
284 : float *ptr_sig[FOA_CHANNELS];
285 : int32_t index[2];
286 : float fac_D, sum_D;
287 : float fac_cost;
288 : float r_sm[16];
289 : float alpha;
290 : float D_tmp[FOA_CHANNELS];
291 : int32_t ivas_total_brate;
292 : float min_dot, dotl, dotr;
293 : float min_dot2;
294 : int16_t bypass_decision;
295 : float dist_alt;
296 :
297 1750 : ivas_total_brate = hEncoderConfig->ivas_total_brate;
298 :
299 : /* if PCA is disabled, just pass-through */
300 1750 : if ( hEncoderConfig->Opt_PCA_ON == 0 )
301 : {
302 : /* write by-pass indicator */
303 0 : push_next_indice( hMetaData, PCA_MODE_INACTIVE, 1 );
304 :
305 0 : return;
306 : }
307 :
308 : /* handle bitrate switching */
309 1750 : if ( ivas_total_brate != PCA_BRATE )
310 : {
311 0 : pca_enc_reset( hPCA );
312 :
313 0 : if ( hEncoderConfig->last_ivas_total_brate != PCA_BRATE )
314 : {
315 0 : eye_matrix( hPCA->mem_eigVec_interp, FOA_CHANNELS, 1.0f );
316 : /* copy input data into output directly as previous frame was already in by-pass mode */
317 0 : for ( k = 0; k < n_channels; k++ )
318 : {
319 : // ToDo: TBV
320 : }
321 : }
322 : else
323 : {
324 : /* set PCA by-pass mode in current frame and interpolate transform as previous frame used PCA */
325 0 : pca_enc_transform( hPCA, ql, qr, data_f, input_frame, n_channels );
326 : }
327 :
328 0 : pca_update_state( hPCA, ql, qr, eigVec, n_channels );
329 :
330 0 : hPCA->prev_bypass_decision = PCA_MODE_INACTIVE;
331 :
332 0 : return;
333 : }
334 :
335 : #ifdef DEBUGGING
336 : assert( ivas_total_brate == PCA_BRATE ); /* the remaining code is defined at 256k where there are 4 dmx channel */
337 : #endif
338 : /*-----------------------------------------------------------------*
339 : * Covariance
340 : *-----------------------------------------------------------------*/
341 :
342 1750 : len_subfr = input_frame / 2;
343 :
344 5250 : for ( i = 0; i < input_frame; i += len_subfr )
345 : {
346 17500 : for ( k = 0; k < FOA_CHANNELS; k++ )
347 : {
348 14000 : ptr_sig[k] = &data_f[k][i];
349 : }
350 :
351 3500 : cov_subfr( ptr_sig, r, n_channels, input_frame / 2 );
352 :
353 3500 : alpha = IVAS_PCA_SM_FAC;
354 59500 : for ( k = 0; k < 16; k++ )
355 : {
356 56000 : r_sm[k] = alpha * r[k] + ( 1.f - alpha ) * hPCA->old_r_sm[k];
357 56000 : hPCA->old_r_sm[k] = r_sm[k];
358 : }
359 : }
360 :
361 : /* conditioning */
362 8750 : for ( k = 0; k < FOA_CHANNELS; k++ )
363 : {
364 7000 : temp = r_sm[k * FOA_CHANNELS + k];
365 :
366 :
367 7000 : if ( temp < IVAS_PCA_COV_THRES )
368 : {
369 0 : temp = IVAS_PCA_COV_THRES;
370 : }
371 7000 : r_sm[k * FOA_CHANNELS + k] = temp; /* pointer reuse */
372 : }
373 :
374 : /*-----------------------------------------------------------------*
375 : * Eigenvalue decomposition
376 : *-----------------------------------------------------------------*/
377 :
378 1750 : eig_qr( r_sm, IVAS_PCA_N_ITER_QR, eigVec, D, n_channels );
379 :
380 : /* force positive eigenvalues */
381 1750 : sum_D = 0.0;
382 8750 : for ( k = 0; k < n_channels; k++ )
383 : {
384 7000 : if ( D[k] < 0.f )
385 : {
386 0 : D[k] = -D[k];
387 :
388 0 : for ( l = 0; l < n_channels; l++ )
389 : {
390 0 : eigVec[l * n_channels + k] = -eigVec[l * n_channels + k];
391 : }
392 : }
393 7000 : sum_D += D[k];
394 : }
395 :
396 : /* normalize */
397 1750 : fac_D = 1 / sum_D;
398 8750 : for ( k = 0; k < n_channels; k++ )
399 : {
400 7000 : D[k] *= fac_D;
401 : }
402 :
403 : /* sorting (sanity check) with SPAR inversion of last two channels */
404 1750 : sort4_D_eigVec( D, eigVec );
405 :
406 : /* normalize and compute amount of decorrelation */
407 1750 : dist_alt = 0.0;
408 8750 : for ( k = 0; k < FOA_CHANNELS; k++ )
409 : {
410 7000 : dist_alt += logf( r[k * FOA_CHANNELS + k] * fac_D + 1e-8f );
411 7000 : dist_alt -= logf( D[k] + 1e-8f );
412 : }
413 :
414 : /*-----------------------------------------------------------------*
415 : * Eigenvector alignment
416 : *-----------------------------------------------------------------*/
417 :
418 1750 : if ( hPCA->prev_bypass_decision == PCA_MODE_ACTIVE )
419 : {
420 : /* compute absolute cost matrix */
421 30 : for ( k = 0; k < n_channels; k++ ) /* column */
422 : {
423 120 : for ( l = 0; l < n_channels; l++ ) /* row */
424 : {
425 96 : fac_cost = hPCA->prev_D[l] * D[k];
426 96 : temp = 0.0;
427 480 : for ( i = 0; i < n_channels; i++ ) /* row */
428 : {
429 384 : temp += hPCA->prev_eigVec[i * n_channels + l] * eigVec[i * n_channels + k];
430 : }
431 96 : if ( temp < 0 )
432 : {
433 38 : temp = -temp;
434 : }
435 96 : cost_mtx[k * n_channels + l] = temp * fac_cost;
436 : }
437 : }
438 :
439 : /* find optimal permutation */
440 6 : exhst_4x4( cost_mtx, path, 1 );
441 : }
442 : else
443 : {
444 : /* no alignment needed if previous PCA is inactive */
445 8720 : for ( i = 0; i < 4; i++ )
446 : {
447 6976 : path[i] = i;
448 : }
449 : }
450 :
451 : /* permute eigenvectors */
452 8750 : for ( k = 0; k < n_channels; k++ )
453 : {
454 7000 : j = path[k];
455 : /* copy j-th column to column k */
456 35000 : for ( l = 0; l < n_channels; l++ )
457 : {
458 28000 : eigVec_tmp[l * n_channels + k] = eigVec[l * n_channels + j];
459 : }
460 7000 : D_tmp[k] = D[j];
461 : }
462 :
463 29750 : for ( k = 0; k < n_channels * n_channels; k++ )
464 : {
465 28000 : eigVec[k] = eigVec_tmp[k];
466 : }
467 :
468 : /* check for sign inversions */
469 8750 : for ( k = 0; k < n_channels; k++ ) /* column */
470 : {
471 7000 : temp = 0.0;
472 35000 : for ( i = 0; i < n_channels; i++ ) /* row */
473 : {
474 28000 : temp += hPCA->prev_eigVec[i * n_channels + k] * eigVec[i * n_channels + k];
475 : }
476 :
477 7000 : if ( temp < 0 )
478 : {
479 9155 : for ( i = 0; i < n_channels; i++ )
480 : {
481 7324 : eigVec[i * n_channels + k] = -eigVec[i * n_channels + k];
482 : }
483 : }
484 : }
485 :
486 : /* force rotation matrix(det = +1) */
487 1750 : det = mat_det4( eigVec );
488 1750 : if ( det < 0 )
489 : {
490 147 : swap_eigvec( eigVec, 2, 3 );
491 147 : temp = D_tmp[3];
492 147 : D_tmp[3] = D_tmp[2];
493 147 : D_tmp[2] = temp;
494 : }
495 :
496 : /* update state */
497 8750 : for ( k = 0; k < n_channels; k++ )
498 : {
499 7000 : hPCA->prev_D[k] = D_tmp[k];
500 : }
501 :
502 : /*-----------------------------------------------------------------*
503 : * Rotation matrix parametrization and quantization
504 : *-----------------------------------------------------------------*/
505 :
506 : /* convert frrm rotation matrix to double quaternion */
507 1750 : mat2dquat( eigVec, ql, qr );
508 :
509 1750 : dotl = dotp( hPCA->prev_ql, ql, 4 );
510 1750 : dotr = dotp( hPCA->prev_qr, qr, 4 );
511 1750 : if ( dotl < dotr )
512 : {
513 743 : min_dot = dotl;
514 : }
515 : else
516 : {
517 1007 : min_dot = dotr;
518 : }
519 :
520 1750 : if ( ql[0] < qr[0] )
521 : {
522 744 : min_dot2 = ql[0];
523 : }
524 : else
525 : {
526 1006 : min_dot2 = qr[0];
527 : }
528 :
529 1750 : bypass_decision = PCA_MODE_ACTIVE;
530 1750 : if ( dist_alt < IVAS_PCA_THRES_DIST_ALT )
531 : {
532 1682 : bypass_decision = PCA_MODE_INACTIVE;
533 : }
534 1750 : if ( min_dot < IVAS_PCA_THRES_MIN_DOT )
535 : {
536 1099 : bypass_decision = PCA_MODE_INACTIVE;
537 : }
538 1750 : if ( min_dot2 < IVAS_PCA_THRES_MIN_DOT2 )
539 : {
540 139 : bypass_decision = PCA_MODE_INACTIVE;
541 : }
542 :
543 : /* if PCA is inactive */
544 1750 : if ( bypass_decision == PCA_MODE_INACTIVE )
545 : {
546 1744 : eye_matrix( eigVec, 4, 1. );
547 1744 : set_zero( ql, 4 );
548 1744 : set_zero( qr, 4 );
549 1744 : ql[0] = 1.;
550 1744 : qr[0] = 1.;
551 :
552 : /* write by-pass indicator */
553 1744 : push_next_indice( hMetaData, PCA_MODE_INACTIVE, 1 );
554 1744 : if ( hPCA->prev_bypass_decision == PCA_MODE_INACTIVE )
555 : {
556 1740 : eye_matrix( hPCA->mem_eigVec_interp, FOA_CHANNELS, 1.0f );
557 : }
558 : else
559 : {
560 : /* set PCA by-pass mode in current frame and interpolate transform as previous frame used PCA */
561 4 : pca_enc_transform( hPCA, ql, qr, data_f, input_frame, n_channels );
562 : }
563 :
564 1744 : pca_update_state( hPCA, ql, qr, eigVec, n_channels );
565 1744 : hPCA->prev_bypass_decision = PCA_MODE_INACTIVE;
566 :
567 1744 : return;
568 : }
569 :
570 : /* force ql to have first component with positive sign */
571 6 : if ( ql[0] < 0 )
572 : {
573 0 : for ( i = 0; i < 4; i++ )
574 : {
575 0 : ql[i] = -ql[i];
576 0 : qr[i] = -qr[i];
577 : }
578 : }
579 :
580 : /* quantize double quaternion */
581 6 : pca_enc_s3( ql, &index[0] );
582 6 : pca_enc_s3( qr, &index[1] );
583 :
584 : /* write bypass flag to bitstream to indicate active mode */
585 6 : push_next_indice( hMetaData, PCA_MODE_ACTIVE, 1 );
586 :
587 6 : ivas_bitstream_write_int32( hMetaData, index[0], IVAS_PCA_QBITS - 1 );
588 6 : ivas_bitstream_write_int32( hMetaData, index[1], IVAS_PCA_QBITS );
589 :
590 : /* transform */
591 6 : pca_enc_transform( hPCA, ql, qr, data_f, input_frame, n_channels );
592 :
593 : /*-----------------------------------------------------------------*
594 : * update state for next frame
595 : *-----------------------------------------------------------------*/
596 :
597 6 : hPCA->prev_bypass_decision = bypass_decision;
598 6 : pca_update_state( hPCA, ql, qr, eigVec, n_channels );
599 :
600 6 : return;
601 : }
|
