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 : /*====================================================================================
34 : EVS Codec 3GPP TS26.443 Nov 04, 2021. Version 12.14.0 / 13.10.0 / 14.6.0 / 15.4.0 / 16.3.0
35 : ====================================================================================*/
36 :
37 : #include <stdint.h>
38 : #include "options.h"
39 : #ifdef DEBUGGING
40 : #include "debug.h"
41 : #endif
42 : #include "cnst.h"
43 : #include "rom_com.h"
44 : #include "prot.h"
45 : #include "wmc_auto.h"
46 : #include <math.h> /* for cosf, sinf */
47 : #include <assert.h>
48 :
49 112431604 : static ivas_error get_edct_table(
50 : const float **edct_table,
51 : const int16_t length )
52 : {
53 : ivas_error error;
54 :
55 112431604 : error = IVAS_ERR_OK;
56 112431604 : *edct_table = NULL;
57 :
58 112431604 : switch ( length )
59 : {
60 336503 : case 1200:
61 336503 : *edct_table = edct_table_600;
62 336503 : break;
63 40644324 : case 960:
64 40644324 : *edct_table = edct_table_480;
65 40644324 : break;
66 25288572 : case 640:
67 25288572 : *edct_table = edct_table_320;
68 25288572 : break;
69 17247685 : case 320:
70 17247685 : *edct_table = edct_table_160;
71 17247685 : break;
72 10624378 : case 256:
73 10624378 : *edct_table = edct_table_128;
74 10624378 : break;
75 1562134 : case 240:
76 1562134 : *edct_table = edct_table_120;
77 1562134 : break;
78 79 : case 200:
79 79 : *edct_table = edct_table_100;
80 79 : break;
81 2034079 : case 160:
82 2034079 : *edct_table = edct_table_80;
83 2034079 : break;
84 392 : case 40:
85 392 : *edct_table = edct_table_20;
86 392 : break;
87 121535 : case 800:
88 121535 : *edct_table = edct_table_400;
89 121535 : break;
90 12301810 : case 512:
91 12301810 : *edct_table = edct_table_256;
92 12301810 : break;
93 715868 : case 480:
94 715868 : *edct_table = edct_table_240;
95 715868 : break;
96 442390 : case 400:
97 442390 : *edct_table = edct_table_200;
98 442390 : break;
99 720172 : case 128:
100 720172 : *edct_table = edct_table_64;
101 720172 : break;
102 391683 : case 80:
103 391683 : *edct_table = edct_table_40;
104 391683 : break;
105 : #ifdef DEBUGGING
106 : default:
107 : return IVAS_ERROR( IVAS_ERR_INTERNAL_FATAL, "edct/edst(): length is not in table!" );
108 : #endif
109 : }
110 :
111 112431604 : return error;
112 : }
113 :
114 :
115 : /*-----------------------------------------------------------------*
116 : * edct()
117 : *
118 : * DCT-IV transform
119 : *-----------------------------------------------------------------*/
120 :
121 100302680 : void edct(
122 : const float *x, /* i : input signal */
123 : float *y, /* o : output transform */
124 : const int16_t length, /* i : length */
125 : const int16_t element_mode /* i : IVAS element mode */
126 : )
127 : {
128 : int16_t i;
129 : float re[L_FRAME_PLUS / 2];
130 : float im[L_FRAME_PLUS / 2];
131 100302680 : const float *edct_table = 0;
132 : float local;
133 :
134 100302680 : get_edct_table( &edct_table, length );
135 :
136 : /* Twiddling and Pre-rotate */
137 30944122564 : for ( i = 0; i < length / 2; i++ )
138 : {
139 30843819884 : re[i] = x[2 * i] * edct_table[i] + x[length - 1 - 2 * i] * edct_table[length / 2 - 1 - i];
140 30843819884 : im[i] = x[length - 1 - 2 * i] * edct_table[i] - x[2 * i] * edct_table[length / 2 - 1 - i];
141 : }
142 :
143 100302680 : if ( element_mode == EVS_MONO )
144 : {
145 933985 : DoFFT( re, im, length / 2 );
146 : }
147 : else
148 : {
149 99368695 : fft( re, im, length / 2, 1 );
150 : }
151 :
152 100302680 : local = ( 0.75f * EVS_PI ) / length;
153 :
154 30944122564 : for ( i = 0; i < length / 2; i++ )
155 : {
156 : float tmp;
157 30843819884 : tmp = re[i] - im[i] * local; /* 3*pi/(4*length) is a constant */
158 30843819884 : im[i] = im[i] + re[i] * local;
159 30843819884 : re[i] = tmp;
160 : }
161 :
162 : /* Post-rotate and obtain the output data */
163 30944122564 : for ( i = 0; i < length / 2; i++ )
164 : {
165 30843819884 : y[2 * i] = re[i] * edct_table[i] + im[i] * edct_table[length / 2 - 1 - i];
166 30843819884 : y[length - 1 - 2 * i] = re[i] * edct_table[length / 2 - 1 - i] - im[i] * edct_table[i];
167 : }
168 :
169 100302680 : return;
170 : }
171 :
172 :
173 : /*-------------------------------------------------------------------------*
174 : * edst()
175 : *
176 : * DST-IV transform
177 : *-------------------------------------------------------------------------*/
178 :
179 12128924 : void edst(
180 : const float *x, /* i : input signal */
181 : float *y, /* o : output transform */
182 : const int16_t length, /* i : length */
183 : const int16_t element_mode /* i : IVAS element mode */
184 : )
185 : {
186 : int16_t i;
187 : float re[L_FRAME_PLUS / 2];
188 : float im[L_FRAME_PLUS / 2];
189 12128924 : const float *edct_table = 0;
190 : float local;
191 :
192 12128924 : get_edct_table( &edct_table, length );
193 :
194 : /* Twiddling and Pre-rotate */
195 4961499532 : for ( i = 0; i < length / 2; i++ )
196 : {
197 4949370608 : re[i] = x[length - 1 - 2 * i] * edct_table[i] + x[2 * i] * edct_table[length / 2 - 1 - i];
198 4949370608 : im[i] = x[2 * i] * edct_table[i] - x[length - 1 - 2 * i] * edct_table[length / 2 - 1 - i];
199 : }
200 :
201 12128924 : if ( element_mode == EVS_MONO )
202 : {
203 202449 : DoFFT( re, im, length / 2 );
204 : }
205 : else
206 : {
207 11926475 : fft( re, im, length / 2, 1 );
208 : }
209 :
210 12128924 : local = ( 0.75f * EVS_PI ) / length;
211 :
212 4961499532 : for ( i = 0; i < length / 2; i++ )
213 : {
214 : float tmp;
215 4949370608 : tmp = re[i] - im[i] * local; /* 3*pi/(4*length) is a constant */
216 4949370608 : im[i] = im[i] + re[i] * local;
217 4949370608 : re[i] = tmp;
218 : }
219 :
220 : /* Post-rotate and obtain the output data */
221 4961499532 : for ( i = 0; i < length / 2; i++ )
222 : {
223 4949370608 : y[2 * i] = re[i] * edct_table[i] + im[i] * edct_table[length / 2 - 1 - i];
224 4949370608 : y[length - 1 - 2 * i] = im[i] * edct_table[i] - re[i] * edct_table[length / 2 - 1 - i];
225 : }
226 :
227 12128924 : return;
228 : }
229 :
230 :
231 : /*-------------------------------------------------------------------------*
232 : * edxt()
233 : *
234 : * DCT/DST-II or III transform (currently also calculates DCT-IV and DST-IV)
235 : *-------------------------------------------------------------------------*/
236 :
237 23614 : void edxt(
238 : const float *x, /* i : input signal */
239 : float *y, /* o : output transform */
240 : const int16_t length, /* i : length */
241 : const uint16_t kernelType, /* i : kernel type (0 - 3) */
242 : const uint16_t synthesis /* i : nonzero for inverse */
243 : )
244 : {
245 23614 : const float pi_len = EVS_PI / length;
246 : int16_t k;
247 :
248 23614 : if ( kernelType == MDST_II || kernelType == MDCT_II )
249 23614 : {
250 23614 : const int16_t Nm1 = length - 1;
251 23614 : const float xSign = 2.f * ( kernelType >> 1 ) - 1.f;
252 23614 : const float scale = 0.5f * pi_len;
253 : float re[L_FRAME_PLUS];
254 : float im[L_FRAME_PLUS];
255 :
256 23614 : if ( !synthesis )
257 : {
258 992422 : for ( k = Nm1 >> 1; k >= 0; k-- ) /* pre-modulation of audio input */
259 : {
260 989040 : re[k] = x[2 * k];
261 989040 : re[Nm1 - k] = x[2 * k + 1] * xSign;
262 989040 : im[k] = im[Nm1 - k] = 0.f;
263 : }
264 :
265 3382 : if ( length == 512 )
266 : {
267 0 : DoRTFTn( re, im, 512 );
268 : }
269 : else /* fft() doesn't support 512 */
270 : {
271 3382 : fft( re, im, length, 1 );
272 : }
273 :
274 3382 : if ( kernelType >> 1 )
275 : {
276 493160 : for ( k = Nm1 >> 1; k > 0; k-- )
277 : {
278 : #ifdef NONBE_FIX_NONBE_BETWEEN_OPTIMIZATION_LEVELS_2
279 491477 : volatile float angle_tmp = scale * k;
280 491477 : const float wRe = cosf( angle_tmp );
281 491477 : const float wIm = sinf( angle_tmp );
282 : #else
283 : const float wRe = cosf( scale * k );
284 : const float wIm = sinf( scale * k );
285 : #endif
286 :
287 491477 : y[k] /*pt 1*/ = wRe * re[k] + wIm * im[k];
288 491477 : y[length - k] = wIm * re[k] - wRe * im[k];
289 : }
290 1683 : y[length >> 1] = re[length >> 1] * sqrtf( 0.5f );
291 : }
292 : else /* forw. DST-II */
293 : {
294 495880 : for ( k = Nm1 >> 1; k > 0; k-- )
295 : {
296 : #ifdef NONBE_FIX_NONBE_BETWEEN_OPTIMIZATION_LEVELS_2
297 494181 : volatile float angle_tmp = scale * k;
298 494181 : const float wRe = cosf( angle_tmp );
299 494181 : const float wIm = sinf( angle_tmp );
300 : #else
301 : const float wRe = cosf( scale * k );
302 : const float wIm = sinf( scale * k );
303 : #endif
304 :
305 494181 : y[Nm1 - k] = wRe * re[k] + wIm * im[k];
306 494181 : y[k - 1] = wIm * re[k] - wRe * im[k];
307 : }
308 1699 : y[Nm1 >> 1] = re[length >> 1] * sqrtf( 0.5f );
309 : }
310 :
311 3382 : y[Nm1 - Nm1 * ( kernelType >> 1 )] = re[0] * 0.5f;
312 : }
313 : else /* inverse II = III */
314 : {
315 20232 : if ( kernelType >> 1 )
316 : {
317 2491824 : for ( k = Nm1 >> 1; k > 0; k-- )
318 : {
319 : #ifdef NONBE_FIX_NONBE_BETWEEN_OPTIMIZATION_LEVELS_2
320 2481744 : volatile float angle_tmp = scale * k;
321 2481744 : const float wRe = cosf( angle_tmp ) * 0.5f;
322 2481744 : const float wIm = sinf( angle_tmp ) * 0.5f;
323 : #else
324 : const float wRe = cosf( scale * k ) * 0.5f;
325 : const float wIm = sinf( scale * k ) * 0.5f;
326 : #endif
327 :
328 2481744 : re[k] = wRe * x[k] + wIm * x[length - k];
329 2481744 : im[k] = wRe * x[length - k] - wIm * x[k];
330 : }
331 10080 : re[length >> 1] = x[length >> 1] * sqrtf( 0.5f );
332 : }
333 : else /* DST type III */
334 : {
335 2493096 : for ( k = Nm1 >> 1; k > 0; k-- )
336 : {
337 : #ifdef NONBE_FIX_NONBE_BETWEEN_OPTIMIZATION_LEVELS_2
338 2482944 : volatile float angle_tmp = scale * k;
339 2482944 : const float wRe = cosf( angle_tmp ) * 0.5f;
340 2482944 : const float wIm = sinf( angle_tmp ) * 0.5f;
341 : #else
342 : const float wRe = cosf( scale * k ) * 0.5f;
343 : const float wIm = sinf( scale * k ) * 0.5f;
344 : #endif
345 :
346 2482944 : re[k] = wRe * x[Nm1 - k] + wIm * x[k - 1];
347 2482944 : im[k] = wRe * x[k - 1] - wIm * x[Nm1 - k];
348 : }
349 10152 : re[length >> 1] = x[Nm1 >> 1] * sqrtf( 0.5f );
350 : }
351 :
352 20232 : re[0] = x[Nm1 - Nm1 * ( kernelType >> 1 )];
353 20232 : im[0] = im[length >> 1] = 0.f;
354 4984920 : for ( k = Nm1 >> 1; k > 0; k-- )
355 : {
356 4964688 : re[length - k] = re[k];
357 4964688 : im[length - k] = -im[k];
358 : }
359 :
360 20232 : if ( length == 512 )
361 : {
362 4698 : DoRTFTn( re, im, 512 );
363 : }
364 : else /* fft() doesn't support 512 */
365 : {
366 15534 : fft( re, im, length, 1 );
367 : }
368 :
369 5005152 : for ( k = Nm1 >> 1; k >= 0; k-- ) /* post-modulation of FFT output */
370 : {
371 4984920 : y[2 * k] = re[k];
372 4984920 : y[2 * k + 1] = re[Nm1 - k] * xSign;
373 : }
374 : }
375 : }
376 : else
377 : {
378 0 : assert( !"Unsupported Kernel type in edxt()" );
379 : }
380 :
381 23614 : v_multc( y, ( kernelType == MDCT_II ? -1.f : 1.f ) * sqrtf( 2.f / length ), y, length );
382 :
383 23614 : return;
384 : }
385 :
386 : /*-----------------------------------------------------------------*
387 : * iedct_short()
388 : *
389 : * Inverse EDCT for short frames
390 : *-----------------------------------------------------------------*/
391 :
392 74904 : void iedct_short(
393 : const float *in, /* i : input vector */
394 : float *out, /* o : output vector */
395 : const int16_t segment_length, /* i : length */
396 : const int16_t element_mode /* i : IVAS element mode */
397 : )
398 : {
399 : float alias[MAX_SEGMENT_LENGTH];
400 : int16_t i;
401 :
402 74904 : edct( in, alias, segment_length / 2, element_mode );
403 :
404 5324584 : for ( i = 0; i < segment_length / 4; i++ )
405 : {
406 5249680 : out[i] = alias[segment_length / 4 + i];
407 5249680 : out[segment_length / 4 + i] = -alias[segment_length / 2 - 1 - i];
408 5249680 : out[segment_length / 2 + i] = -alias[segment_length / 4 - 1 - i];
409 5249680 : out[3 * segment_length / 4 + i] = -alias[i];
410 : }
411 :
412 74904 : return;
413 : }
|
