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 "prot.h"
43 : #include "rom_com.h"
44 : #include "wmc_auto.h"
45 :
46 : /*---------------------------------------------------------------------*
47 : * Local constants
48 : *---------------------------------------------------------------------*/
49 :
50 : #define N_MAX_FFT 1024
51 : #define N_MAX_DIV2 ( N_MAX_FFT >> 1 )
52 : #define N_MAX_DIV4 ( N_MAX_DIV2 >> 1 )
53 :
54 : /*---------------------------------------------------------------------*
55 : * fft_rel()
56 : *
57 : * Computes the split-radix FFT in place for the real-valued
58 : * signal x of length n. The algorithm has been ported from
59 : * the Fortran code of [1].
60 : *
61 : * The function needs sine and cosine tables t_sin and t_cos,
62 : * and the constant N_MAX_FFT. The table entries are defined as
63 : * sin(2*pi*i) and cos(2*pi*i) for i = 0, 1, ..., N_MAX_FFT-1. The
64 : * implementation assumes that any entry will not be needed
65 : * outside the tables. Therefore, N_MAX_FFT and n must be properly
66 : * set. The function has been tested with the values n = 16,
67 : * 32, 64, 128, 256, and N_MAX_FFT = 1280.
68 : *
69 : * References
70 : * [1] H.V. Sorensen, D.L. Jones, M.T. Heideman, C.S. Burrus,
71 : * "Real-valued fast Fourier transform algorithm," IEEE
72 : * Trans. on Signal Processing, Vol.35, No.6, pp 849-863,
73 : * 1987.
74 : *
75 : * OUTPUT
76 : * x[0:n-1] Transform coeffients in the order re[0], re[1],
77 : * ..., re[n/2], im[n/2-1], ..., im[1].
78 : *---------------------------------------------------------------------*/
79 :
80 40860621 : void fft_rel(
81 : float x[], /* i/o: input/output vector */
82 : const int16_t n, /* i : vector length */
83 : const int16_t m /* i : log2 of vector length */
84 : )
85 : {
86 : int16_t i, j, k, n1, n2, n4;
87 : int16_t step;
88 : float xt, t1, t2;
89 : float *x0, *x1, *x2;
90 : float *xi2, *xi3, *xi4, *xi1;
91 : const float *s, *c;
92 : const int16_t *idx;
93 :
94 : /* !!!! NMAX = 256 is hardcoded here (similar optimizations should be done for NMAX > 256) !!! */
95 :
96 : float *x2even, *x2odd;
97 : float temp[512];
98 :
99 40860621 : if ( n == 128 || n == 256 || n == 512 )
100 : {
101 39316097 : idx = fft256_read_indexes;
102 :
103 : /* Combined Digit reverse counter & Length two butterflies */
104 39316097 : if ( n == 128 )
105 : {
106 13822 : x2 = temp;
107 898430 : for ( i = 0; i < 64; i++ )
108 : {
109 884608 : j = *idx++;
110 884608 : k = *idx++;
111 :
112 884608 : *x2++ = x[j >> 1] + x[k >> 1];
113 884608 : *x2++ = x[j >> 1] - x[k >> 1];
114 : }
115 : }
116 39302275 : else if ( n == 256 )
117 : {
118 37709277 : x2 = temp;
119 4864496733 : for ( i = 0; i < 128; i++ )
120 : {
121 4826787456 : j = *idx++;
122 4826787456 : k = *idx++;
123 :
124 4826787456 : *x2++ = x[j] + x[k];
125 4826787456 : *x2++ = x[j] - x[k];
126 : }
127 : }
128 1592998 : else if ( n == 512 )
129 : {
130 1592998 : x2even = temp;
131 1592998 : x2odd = temp + 256;
132 :
133 205496742 : for ( i = 0; i < 128; i++ )
134 : {
135 203903744 : j = 2 * *idx++;
136 203903744 : k = 2 * *idx++;
137 :
138 203903744 : *x2even++ = x[j] + x[k];
139 203903744 : *x2even++ = x[j] - x[k];
140 203903744 : *x2odd++ = x[++j] + x[++k];
141 203903744 : *x2odd++ = x[j] - x[k];
142 : }
143 : }
144 :
145 : /*-----------------------------------------------------------------*
146 : * 1st Stage Loop has been Unrolled because n4 is '1' and that
147 : * allows the elimination of the 'for_ (j = 1; j < n4; j++)' loop
148 : * and the associated pointers initialization.
149 : * Also, it allows to Put the Data from 'temp' back into 'x' due
150 : * to the previous Combined Digit Reverse and Length two butterflies
151 : *-----------------------------------------------------------------*/
152 :
153 : /*for_ (k = 2; k < 3; k++)*/
154 : {
155 39316097 : x0 = temp;
156 39316097 : x1 = x0 + 2;
157 39316097 : x2 = x;
158 :
159 2657055873 : for ( i = 0; i < n; i += 4 )
160 : {
161 2617739776 : *x2++ = *x0++ + *x1; /* x[i] = xt + x[i+n2]; */
162 2617739776 : *x2++ = *x0;
163 2617739776 : *x2++ = *--x0 - *x1++; /* x[i+n2] = xt - x[i+n2]; */
164 2617739776 : *x2++ = -*x1; /* x[i+n2+n4] = -x[i+n2+n4]; */
165 :
166 2617739776 : x0 += 4;
167 2617739776 : x1 += 3; /* x1 has already advanced */
168 : }
169 : }
170 : }
171 : else
172 : {
173 : /*-----------------------------------------------------------------*
174 : * Digit reverse counter
175 : *-----------------------------------------------------------------*/
176 :
177 1544524 : j = 0;
178 1544524 : x0 = &x[0];
179 59559796 : for ( i = 0; i < n - 1; i++ )
180 : {
181 58015272 : if ( i < j )
182 : {
183 25627885 : xt = x[j];
184 25627885 : x[j] = *x0;
185 25627885 : *x0 = xt;
186 : }
187 58015272 : x0++;
188 58015272 : k = n / 2;
189 109485376 : while ( k <= j )
190 : {
191 51470104 : j -= k;
192 51470104 : k = k >> 1;
193 : }
194 58015272 : j += k;
195 : }
196 :
197 : /*-----------------------------------------------------------------*
198 : * Length two butterflies
199 : *-----------------------------------------------------------------*/
200 :
201 1544524 : x0 = &x[0];
202 1544524 : x1 = &x[1];
203 31324422 : for ( i = 0; i < n / 2; i++ )
204 : {
205 29779898 : *x1 = *x0 - *x1;
206 29779898 : *x0 = *x0 * 2 - *x1;
207 :
208 29779898 : x0++;
209 29779898 : x0++;
210 29779898 : x1++;
211 29779898 : x1++;
212 : }
213 :
214 : /*-----------------------------------------------------------------*
215 : * 1st Stage Loop has been Unrolled because n4 is '1' and that
216 : * allows the elimination of the 'for_ (j = 1; j < n4; j++)' loop
217 : * and the associated pointers initialization.
218 : *-----------------------------------------------------------------*/
219 :
220 : /* for_ (k = 2; k < 3; k++) */
221 : {
222 1544524 : x0 = x;
223 1544524 : x1 = x0 + 2;
224 :
225 16434473 : for ( i = 0; i < n; i += 4 )
226 : {
227 14889949 : *x1 = *x0 - *x1; /* x[i+n2] = xt - x[i+n2]; */
228 14889949 : *x0 = *x0 * 2 - *x1++; /* x[i] = xt + x[i+n2]; */
229 14889949 : *x1 = -*x1; /* x[i+n2+n4] = -x[i+n2+n4]; */
230 :
231 14889949 : x0 += 4;
232 14889949 : x1 += 3; /* x1 has already advanced */
233 : }
234 : }
235 : }
236 :
237 : /*-----------------------------------------------------------------*
238 : * Other butterflies
239 : *
240 : * The implementation described in [1] has been changed by using
241 : * table lookup for evaluating sine and cosine functions. The
242 : * variable ind and its increment step are needed to access table
243 : * entries. Note that this implementation assumes n4 to be so
244 : * small that ind will never exceed the table. Thus the input
245 : * argument n and the constant N_MAX_FFT must be set properly.
246 : *-----------------------------------------------------------------*/
247 :
248 40860621 : n4 = 1;
249 40860621 : n2 = 2;
250 40860621 : n1 = 4;
251 :
252 40860621 : step = N_MAX_DIV4;
253 :
254 281792499 : for ( k = 3; k <= m; k++ )
255 : {
256 240931878 : step >>= 1;
257 240931878 : n4 <<= 1;
258 240931878 : n2 <<= 1;
259 240931878 : n1 <<= 1;
260 :
261 240931878 : x0 = x;
262 240931878 : x1 = x0 + n2;
263 240931878 : x2 = x1 + n4;
264 :
265 2832700982 : for ( i = 0; i < n; i += n1 )
266 : {
267 2591769104 : *x1 = *x0 - *x1; /* x[i+n2] = xt - x[i+n2]; */
268 2591769104 : *x0 = *x0 * 2 - *x1; /* x[i] = xt + x[i+n2]; */
269 2591769104 : *x2 = -*x2; /* x[i+n2+n4] = -x[i+n2+n4]; */
270 :
271 2591769104 : s = sincos_t_ext;
272 2591769104 : c = s + N_MAX_FFT / 4; /* 1024/4 = 256, 256/4=64 */
273 2591769104 : xi1 = x0;
274 2591769104 : xi3 = xi1 + n2;
275 2591769104 : xi2 = xi3;
276 2591769104 : x0 += n1;
277 2591769104 : xi4 = x0;
278 :
279 15968483136 : for ( j = 1; j < n4; j++ )
280 : {
281 13376714032 : xi3++;
282 13376714032 : xi1++;
283 13376714032 : xi4--;
284 13376714032 : xi2--;
285 13376714032 : c += step;
286 13376714032 : s += step; /* autoincrement by ar0 */
287 :
288 13376714032 : t1 = *xi3 * *c + *xi4 * *s; /* t1 = *xi3**(pt_c+ind) + *xi4**(pt_s+ind); */
289 13376714032 : t2 = *xi3 * *s - *xi4 * *c; /* t2 = *xi3**(pt_s+ind) - *xi4**(pt_c+ind); */
290 :
291 13376714032 : *xi4 = *xi2 - t2;
292 13376714032 : *xi2 = *xi1 - t1;
293 13376714032 : *xi1 = *xi1 * 2 - *xi2;
294 13376714032 : *xi3 = -2 * t2 - *xi4;
295 : }
296 :
297 2591769104 : x1 += n1;
298 2591769104 : x2 += n1;
299 : }
300 : }
301 :
302 40860621 : return;
303 : }
|
