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 :
48 102408318 : static ivas_error get_edct_table(
49 : const float **edct_table,
50 : const int16_t length )
51 : {
52 : ivas_error error;
53 :
54 102408318 : error = IVAS_ERR_OK;
55 102408318 : *edct_table = NULL;
56 :
57 102408318 : switch ( length )
58 : {
59 288267 : case 1200:
60 288267 : *edct_table = edct_table_600;
61 288267 : break;
62 33718404 : case 960:
63 33718404 : *edct_table = edct_table_480;
64 33718404 : break;
65 24452275 : case 640:
66 24452275 : *edct_table = edct_table_320;
67 24452275 : break;
68 21483255 : case 320:
69 21483255 : *edct_table = edct_table_160;
70 21483255 : break;
71 8448442 : case 256:
72 8448442 : *edct_table = edct_table_128;
73 8448442 : break;
74 1359124 : case 240:
75 1359124 : *edct_table = edct_table_120;
76 1359124 : break;
77 91 : case 200:
78 91 : *edct_table = edct_table_100;
79 91 : break;
80 1973064 : case 160:
81 1973064 : *edct_table = edct_table_80;
82 1973064 : break;
83 488 : case 40:
84 488 : *edct_table = edct_table_20;
85 488 : break;
86 117967 : case 800:
87 117967 : *edct_table = edct_table_400;
88 117967 : break;
89 8207475 : case 512:
90 8207475 : *edct_table = edct_table_256;
91 8207475 : break;
92 628913 : case 480:
93 628913 : *edct_table = edct_table_240;
94 628913 : break;
95 481317 : case 400:
96 481317 : *edct_table = edct_table_200;
97 481317 : break;
98 513665 : case 128:
99 513665 : *edct_table = edct_table_64;
100 513665 : break;
101 735571 : case 80:
102 735571 : *edct_table = edct_table_40;
103 735571 : break;
104 : #ifdef DEBUGGING
105 : default:
106 : return IVAS_ERROR( IVAS_ERR_INTERNAL_FATAL, "edct/edst(): length is not in table!" );
107 : #endif
108 : }
109 :
110 102408318 : return error;
111 : }
112 :
113 :
114 : /*-----------------------------------------------------------------*
115 : * edct()
116 : *
117 : * DCT-IV transform
118 : *-----------------------------------------------------------------*/
119 :
120 88937370 : void edct(
121 : const float *x, /* i : input signal */
122 : float *y, /* o : output transform */
123 : const int16_t length, /* i : length */
124 : const int16_t element_mode /* i : IVAS element mode */
125 : )
126 : {
127 : int16_t i;
128 : float re[L_FRAME_PLUS / 2];
129 : float im[L_FRAME_PLUS / 2];
130 88937370 : const float *edct_table = 0;
131 : float local;
132 :
133 88937370 : get_edct_table( &edct_table, length );
134 :
135 : /* Twiddling and Pre-rotate */
136 26099998474 : for ( i = 0; i < length / 2; i++ )
137 : {
138 26011061104 : re[i] = x[2 * i] * edct_table[i] + x[length - 1 - 2 * i] * edct_table[length / 2 - 1 - i];
139 26011061104 : im[i] = x[length - 1 - 2 * i] * edct_table[i] - x[2 * i] * edct_table[length / 2 - 1 - i];
140 : }
141 :
142 88937370 : if ( element_mode == EVS_MONO )
143 : {
144 873200 : DoFFT( re, im, length / 2 );
145 : }
146 : else
147 : {
148 88064170 : fft( re, im, length / 2, 1 );
149 : }
150 :
151 88937370 : local = ( 0.75f * EVS_PI ) / length;
152 :
153 26099998474 : for ( i = 0; i < length / 2; i++ )
154 : {
155 : float tmp;
156 26011061104 : tmp = re[i] - im[i] * local; /* 3*pi/(4*length) is a constant */
157 26011061104 : im[i] = im[i] + re[i] * local;
158 26011061104 : re[i] = tmp;
159 : }
160 :
161 : /* Post-rotate and obtain the output data */
162 26099998474 : for ( i = 0; i < length / 2; i++ )
163 : {
164 26011061104 : y[2 * i] = re[i] * edct_table[i] + im[i] * edct_table[length / 2 - 1 - i];
165 26011061104 : y[length - 1 - 2 * i] = re[i] * edct_table[length / 2 - 1 - i] - im[i] * edct_table[i];
166 : }
167 :
168 88937370 : return;
169 : }
170 :
171 :
172 : /*-------------------------------------------------------------------------*
173 : * edst()
174 : *
175 : * DST-IV transform
176 : *-------------------------------------------------------------------------*/
177 :
178 13470948 : void edst(
179 : const float *x, /* i : input signal */
180 : float *y, /* o : output transform */
181 : const int16_t length, /* i : length */
182 : const int16_t element_mode /* i : IVAS element mode */
183 : )
184 : {
185 : int16_t i;
186 : float re[L_FRAME_PLUS / 2];
187 : float im[L_FRAME_PLUS / 2];
188 13470948 : const float *edct_table = 0;
189 : float local;
190 :
191 13470948 : get_edct_table( &edct_table, length );
192 :
193 : /* Twiddling and Pre-rotate */
194 5482412520 : for ( i = 0; i < length / 2; i++ )
195 : {
196 5468941572 : re[i] = x[length - 1 - 2 * i] * edct_table[i] + x[2 * i] * edct_table[length / 2 - 1 - i];
197 5468941572 : im[i] = x[2 * i] * edct_table[i] - x[length - 1 - 2 * i] * edct_table[length / 2 - 1 - i];
198 : }
199 :
200 13470948 : if ( element_mode == EVS_MONO )
201 : {
202 196852 : DoFFT( re, im, length / 2 );
203 : }
204 : else
205 : {
206 13274096 : fft( re, im, length / 2, 1 );
207 : }
208 :
209 13470948 : local = ( 0.75f * EVS_PI ) / length;
210 :
211 5482412520 : for ( i = 0; i < length / 2; i++ )
212 : {
213 : float tmp;
214 5468941572 : tmp = re[i] - im[i] * local; /* 3*pi/(4*length) is a constant */
215 5468941572 : im[i] = im[i] + re[i] * local;
216 5468941572 : re[i] = tmp;
217 : }
218 :
219 : /* Post-rotate and obtain the output data */
220 5482412520 : for ( i = 0; i < length / 2; i++ )
221 : {
222 5468941572 : y[2 * i] = re[i] * edct_table[i] + im[i] * edct_table[length / 2 - 1 - i];
223 5468941572 : y[length - 1 - 2 * i] = im[i] * edct_table[i] - re[i] * edct_table[length / 2 - 1 - i];
224 : }
225 :
226 13470948 : return;
227 : }
228 :
229 : #define FAST_EDXT /* optimized FFT-based DCT/DST algorithm */
230 :
231 : /*-------------------------------------------------------------------------*
232 : * edxt()
233 : *
234 : * DCT/DST-II or III transform (currently also calculates DCT-IV and DST-IV)
235 : *-------------------------------------------------------------------------*/
236 :
237 16712 : 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 16712 : const float pi_len = EVS_PI / length;
246 : int16_t k, m;
247 :
248 : #ifdef FAST_EDXT
249 16712 : if ( kernelType == MDST_II || kernelType == MDCT_II )
250 16712 : {
251 16712 : const int16_t Nm1 = length - 1;
252 16712 : const float xSign = 2.f * ( kernelType >> 1 ) - 1.f;
253 16712 : const float scale = 0.5f * pi_len;
254 : float re[L_FRAME_PLUS];
255 : float im[L_FRAME_PLUS];
256 :
257 16712 : if ( !synthesis )
258 : {
259 1194200 : for ( k = Nm1 >> 1; k >= 0; k-- ) /* pre-modulation of audio input */
260 : {
261 1190160 : re[k] = x[2 * k];
262 1190160 : re[Nm1 - k] = x[2 * k + 1] * xSign;
263 1190160 : im[k] = im[Nm1 - k] = 0.f;
264 : }
265 :
266 4040 : if ( length == 512 )
267 : {
268 0 : DoRTFTn( re, im, 512 );
269 : }
270 : else /* fft() doesn't support 512 */
271 : {
272 4040 : fft( re, im, length, 1 );
273 : }
274 :
275 4040 : if ( kernelType >> 1 )
276 : {
277 593400 : for ( k = Nm1 >> 1; k > 0; k-- )
278 : {
279 591389 : const float wRe = cosf( scale * k );
280 591389 : const float wIm = sinf( scale * k );
281 :
282 591389 : y[k] /*pt 1*/ = wRe * re[k] + wIm * im[k];
283 591389 : y[length - k] = wIm * re[k] - wRe * im[k];
284 : }
285 2011 : y[length >> 1] = re[length >> 1] * sqrtf( 0.5f );
286 : }
287 : else /* forw. DST-II */
288 : {
289 596760 : for ( k = Nm1 >> 1; k > 0; k-- )
290 : {
291 594731 : const float wRe = cosf( scale * k );
292 594731 : const float wIm = sinf( scale * k );
293 :
294 594731 : y[Nm1 - k] = wRe * re[k] + wIm * im[k];
295 594731 : y[k - 1] = wIm * re[k] - wRe * im[k];
296 : }
297 2029 : y[Nm1 >> 1] = re[length >> 1] * sqrtf( 0.5f );
298 : }
299 :
300 4040 : y[Nm1 - Nm1 * ( kernelType >> 1 )] = re[0] * 0.5f;
301 : }
302 : else /* inverse II = III */
303 : {
304 12672 : if ( kernelType >> 1 )
305 : {
306 1681392 : for ( k = Nm1 >> 1; k > 0; k-- )
307 : {
308 1674926 : const float wRe = cosf( scale * k ) * 0.5f;
309 1674926 : const float wIm = sinf( scale * k ) * 0.5f;
310 :
311 1674926 : re[k] = wRe * x[k] + wIm * x[length - k];
312 1674926 : im[k] = wRe * x[length - k] - wIm * x[k];
313 : }
314 6466 : re[length >> 1] = x[length >> 1] * sqrtf( 0.5f );
315 : }
316 : else /* DST type III */
317 : {
318 1557928 : for ( k = Nm1 >> 1; k > 0; k-- )
319 : {
320 1551722 : const float wRe = cosf( scale * k ) * 0.5f;
321 1551722 : const float wIm = sinf( scale * k ) * 0.5f;
322 :
323 1551722 : re[k] = wRe * x[Nm1 - k] + wIm * x[k - 1];
324 1551722 : im[k] = wRe * x[k - 1] - wIm * x[Nm1 - k];
325 : }
326 6206 : re[length >> 1] = x[Nm1 >> 1] * sqrtf( 0.5f );
327 : }
328 :
329 12672 : re[0] = x[Nm1 - Nm1 * ( kernelType >> 1 )];
330 12672 : im[0] = im[length >> 1] = 0.f;
331 3239320 : for ( k = Nm1 >> 1; k > 0; k-- )
332 : {
333 3226648 : re[length - k] = re[k];
334 3226648 : im[length - k] = -im[k];
335 : }
336 :
337 12672 : if ( length == 512 )
338 : {
339 3436 : DoRTFTn( re, im, 512 );
340 : }
341 : else /* fft() doesn't support 512 */
342 : {
343 9236 : fft( re, im, length, 1 );
344 : }
345 :
346 3251992 : for ( k = Nm1 >> 1; k >= 0; k-- ) /* post-modulation of FFT output */
347 : {
348 3239320 : y[2 * k] = re[k];
349 3239320 : y[2 * k + 1] = re[Nm1 - k] * xSign;
350 : }
351 : }
352 : }
353 : else
354 : #endif
355 0 : if ( kernelType & 1 ) /* DST */
356 : {
357 0 : const float offK = ( kernelType == MDST_II && synthesis ? 0.5f : 1.0f - 0.5f * ( kernelType >> 1 ) );
358 0 : const float offM = ( kernelType == MDST_II && synthesis ? 1.0f : 0.5f );
359 :
360 0 : for ( k = 0; k < length; k++ )
361 : {
362 0 : y[k] = 0.f;
363 0 : for ( m = 0; m < length; m++ )
364 : {
365 0 : y[k] += x[m] * sinf( pi_len * ( m + offM ) * ( k + offK ) );
366 : }
367 : }
368 0 : if ( offK == 1.f )
369 : {
370 0 : y[length - 1] *= 0.5f; /* scale Nyquist sample */
371 : }
372 : }
373 : else /* kernelType 0, 2: DCT */
374 : {
375 0 : const float offK = ( kernelType == MDCT_II && synthesis ? 0.5f : 0.5f - 0.5f * ( kernelType >> 1 ) );
376 0 : const float offM = ( kernelType == MDCT_II && synthesis ? 0.0f : 0.5f );
377 :
378 0 : for ( k = 0; k < length; k++ )
379 : {
380 0 : y[k] = 0.f;
381 0 : for ( m = 0; m < length; m++ )
382 : {
383 0 : y[k] += x[m] * cosf( pi_len * ( m + offM ) * ( k + offK ) );
384 : }
385 : }
386 0 : if ( offK == 0.f )
387 : {
388 0 : y[0] *= 0.5f; /* scale lowest (i.e. DC) sample */
389 : }
390 : }
391 :
392 16712 : v_multc( y, ( kernelType == MDCT_II ? -1.f : 1.f ) * sqrtf( 2.f / length ), y, length );
393 :
394 16712 : return;
395 : }
396 :
397 : /*-----------------------------------------------------------------*
398 : * iedct_short()
399 : *
400 : * Inverse EDCT for short frames
401 : *-----------------------------------------------------------------*/
402 :
403 59684 : void iedct_short(
404 : const float *in, /* i : input vector */
405 : float *out, /* o : output vector */
406 : const int16_t segment_length, /* i : length */
407 : const int16_t element_mode /* i : IVAS element mode */
408 : )
409 : {
410 : float alias[MAX_SEGMENT_LENGTH];
411 : int16_t i;
412 :
413 59684 : edct( in, alias, segment_length / 2, element_mode );
414 :
415 4020724 : for ( i = 0; i < segment_length / 4; i++ )
416 : {
417 3961040 : out[i] = alias[segment_length / 4 + i];
418 3961040 : out[segment_length / 4 + i] = -alias[segment_length / 2 - 1 - i];
419 3961040 : out[segment_length / 2 + i] = -alias[segment_length / 4 - 1 - i];
420 3961040 : out[3 * segment_length / 4 + i] = -alias[i];
421 : }
422 :
423 59684 : return;
424 : }
|
