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 3194115 : 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 3194115 : if ( n == 128 || n == 256 || n == 512 )
100 : {
101 3027355 : idx = fft256_read_indexes;
102 :
103 : /* Combined Digit reverse counter & Length two butterflies */
104 3027355 : if ( n == 128 )
105 : {
106 1224 : x2 = temp;
107 79560 : for ( i = 0; i < 64; i++ )
108 : {
109 78336 : j = *idx++;
110 78336 : k = *idx++;
111 :
112 78336 : *x2++ = x[j >> 1] + x[k >> 1];
113 78336 : *x2++ = x[j >> 1] - x[k >> 1];
114 : }
115 : }
116 3026131 : else if ( n == 256 )
117 : {
118 2880301 : x2 = temp;
119 371558829 : for ( i = 0; i < 128; i++ )
120 : {
121 368678528 : j = *idx++;
122 368678528 : k = *idx++;
123 :
124 368678528 : *x2++ = x[j] + x[k];
125 368678528 : *x2++ = x[j] - x[k];
126 : }
127 : }
128 145830 : else if ( n == 512 )
129 : {
130 145830 : x2even = temp;
131 145830 : x2odd = temp + 256;
132 :
133 18812070 : for ( i = 0; i < 128; i++ )
134 : {
135 18666240 : j = 2 * *idx++;
136 18666240 : k = 2 * *idx++;
137 :
138 18666240 : *x2even++ = x[j] + x[k];
139 18666240 : *x2even++ = x[j] - x[k];
140 18666240 : *x2odd++ = x[++j] + x[++k];
141 18666240 : *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 3027355 : x0 = temp;
156 3027355 : x1 = x0 + 2;
157 3027355 : x2 = x;
158 :
159 206072027 : for ( i = 0; i < n; i += 4 )
160 : {
161 203044672 : *x2++ = *x0++ + *x1; /* x[i] = xt + x[i+n2]; */
162 203044672 : *x2++ = *x0;
163 203044672 : *x2++ = *--x0 - *x1++; /* x[i+n2] = xt - x[i+n2]; */
164 203044672 : *x2++ = -*x1; /* x[i+n2+n4] = -x[i+n2+n4]; */
165 :
166 203044672 : x0 += 4;
167 203044672 : x1 += 3; /* x1 has already advanced */
168 : }
169 : }
170 : }
171 : else
172 : {
173 : /*-----------------------------------------------------------------*
174 : * Digit reverse counter
175 : *-----------------------------------------------------------------*/
176 :
177 166760 : j = 0;
178 166760 : x0 = &x[0];
179 7324880 : for ( i = 0; i < n - 1; i++ )
180 : {
181 7158120 : if ( i < j )
182 : {
183 3168404 : xt = x[j];
184 3168404 : x[j] = *x0;
185 3168404 : *x0 = xt;
186 : }
187 7158120 : x0++;
188 7158120 : k = n / 2;
189 13548224 : while ( k <= j )
190 : {
191 6390104 : j -= k;
192 6390104 : k = k >> 1;
193 : }
194 7158120 : j += k;
195 : }
196 :
197 : /*-----------------------------------------------------------------*
198 : * Length two butterflies
199 : *-----------------------------------------------------------------*/
200 :
201 166760 : x0 = &x[0];
202 166760 : x1 = &x[1];
203 3829200 : for ( i = 0; i < n / 2; i++ )
204 : {
205 3662440 : *x1 = *x0 - *x1;
206 3662440 : *x0 = *x0 * 2 - *x1;
207 :
208 3662440 : x0++;
209 3662440 : x0++;
210 3662440 : x1++;
211 3662440 : 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 166760 : x0 = x;
223 166760 : x1 = x0 + 2;
224 :
225 1997980 : for ( i = 0; i < n; i += 4 )
226 : {
227 1831220 : *x1 = *x0 - *x1; /* x[i+n2] = xt - x[i+n2]; */
228 1831220 : *x0 = *x0 * 2 - *x1++; /* x[i] = xt + x[i+n2]; */
229 1831220 : *x1 = -*x1; /* x[i+n2+n4] = -x[i+n2+n4]; */
230 :
231 1831220 : x0 += 4;
232 1831220 : 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 3194115 : n4 = 1;
249 3194115 : n2 = 2;
250 3194115 : n1 = 4;
251 :
252 3194115 : step = N_MAX_DIV4;
253 :
254 21937347 : for ( k = 3; k <= m; k++ )
255 : {
256 18743232 : step >>= 1;
257 18743232 : n4 <<= 1;
258 18743232 : n2 <<= 1;
259 18743232 : n1 <<= 1;
260 :
261 18743232 : x0 = x;
262 18743232 : x1 = x0 + n2;
263 18743232 : x2 = x1 + n4;
264 :
265 220425009 : for ( i = 0; i < n; i += n1 )
266 : {
267 201681777 : *x1 = *x0 - *x1; /* x[i+n2] = xt - x[i+n2]; */
268 201681777 : *x0 = *x0 * 2 - *x1; /* x[i] = xt + x[i+n2]; */
269 201681777 : *x2 = -*x2; /* x[i+n2+n4] = -x[i+n2+n4]; */
270 :
271 201681777 : s = sincos_t_ext;
272 201681777 : c = s + N_MAX_FFT / 4; /* 1024/4 = 256, 256/4=64 */
273 201681777 : xi1 = x0;
274 201681777 : xi3 = xi1 + n2;
275 201681777 : xi2 = xi3;
276 201681777 : x0 += n1;
277 201681777 : xi4 = x0;
278 :
279 1244146560 : for ( j = 1; j < n4; j++ )
280 : {
281 1042464783 : xi3++;
282 1042464783 : xi1++;
283 1042464783 : xi4--;
284 1042464783 : xi2--;
285 1042464783 : c += step;
286 1042464783 : s += step; /* autoincrement by ar0 */
287 :
288 1042464783 : t1 = *xi3 * *c + *xi4 * *s; /* t1 = *xi3**(pt_c+ind) + *xi4**(pt_s+ind); */
289 1042464783 : t2 = *xi3 * *s - *xi4 * *c; /* t2 = *xi3**(pt_s+ind) - *xi4**(pt_c+ind); */
290 :
291 1042464783 : *xi4 = *xi2 - t2;
292 1042464783 : *xi2 = *xi1 - t1;
293 1042464783 : *xi1 = *xi1 * 2 - *xi2;
294 1042464783 : *xi3 = -2 * t2 - *xi4;
295 : }
296 :
297 201681777 : x1 += n1;
298 201681777 : x2 += n1;
299 : }
300 : }
301 :
302 3194115 : return;
303 : }
|
