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