/proc/self/cwd/libfaad/cfft.c
Line | Count | Source (jump to first uncovered line) |
1 | | /* |
2 | | ** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding |
3 | | ** Copyright (C) 2003-2005 M. Bakker, Nero AG, http://www.nero.com |
4 | | ** |
5 | | ** This program is free software; you can redistribute it and/or modify |
6 | | ** it under the terms of the GNU General Public License as published by |
7 | | ** the Free Software Foundation; either version 2 of the License, or |
8 | | ** (at your option) any later version. |
9 | | ** |
10 | | ** This program is distributed in the hope that it will be useful, |
11 | | ** but WITHOUT ANY WARRANTY; without even the implied warranty of |
12 | | ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
13 | | ** GNU General Public License for more details. |
14 | | ** |
15 | | ** You should have received a copy of the GNU General Public License |
16 | | ** along with this program; if not, write to the Free Software |
17 | | ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
18 | | ** |
19 | | ** Any non-GPL usage of this software or parts of this software is strictly |
20 | | ** forbidden. |
21 | | ** |
22 | | ** The "appropriate copyright message" mentioned in section 2c of the GPLv2 |
23 | | ** must read: "Code from FAAD2 is copyright (c) Nero AG, www.nero.com" |
24 | | ** |
25 | | ** Commercial non-GPL licensing of this software is possible. |
26 | | ** For more info contact Nero AG through Mpeg4AAClicense@nero.com. |
27 | | ** |
28 | | ** $Id: cfft.c,v 1.35 2007/11/01 12:33:29 menno Exp $ |
29 | | **/ |
30 | | |
31 | | /* |
32 | | * Algorithmically based on Fortran-77 FFTPACK |
33 | | * by Paul N. Swarztrauber(Version 4, 1985). |
34 | | * |
35 | | * Does even sized fft only |
36 | | */ |
37 | | |
38 | | /* isign is +1 for backward and -1 for forward transforms */ |
39 | | |
40 | | #include "common.h" |
41 | | #include "structs.h" |
42 | | |
43 | | #include <stdlib.h> |
44 | | |
45 | | #include "cfft.h" |
46 | | #include "cfft_tab.h" |
47 | | |
48 | | |
49 | | /* static function declarations */ |
50 | | static void passf2pos(const uint16_t ido, const uint16_t l1, const complex_t *cc, |
51 | | complex_t *ch, const complex_t *wa); |
52 | | static void passf2neg(const uint16_t ido, const uint16_t l1, const complex_t *cc, |
53 | | complex_t *ch, const complex_t *wa); |
54 | | static void passf3(const uint16_t ido, const uint16_t l1, const complex_t *cc, |
55 | | complex_t *ch, const complex_t *wa1, const complex_t *wa2, const int8_t isign); |
56 | | static void passf4pos(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch, |
57 | | const complex_t *wa1, const complex_t *wa2, const complex_t *wa3); |
58 | | static void passf4neg(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch, |
59 | | const complex_t *wa1, const complex_t *wa2, const complex_t *wa3); |
60 | | static void passf5(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch, |
61 | | const complex_t *wa1, const complex_t *wa2, const complex_t *wa3, |
62 | | const complex_t *wa4, const int8_t isign); |
63 | | static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac); |
64 | | |
65 | | |
66 | | /*---------------------------------------------------------------------- |
67 | | passf2, passf3, passf4, passf5. Complex FFT passes fwd and bwd. |
68 | | ----------------------------------------------------------------------*/ |
69 | | |
70 | | static void passf2pos(const uint16_t ido, const uint16_t l1, const complex_t *cc, |
71 | | complex_t *ch, const complex_t *wa) |
72 | 180k | { |
73 | 180k | uint16_t i, k, ah, ac; |
74 | | |
75 | 180k | if (ido == 1) |
76 | 0 | { |
77 | 0 | #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION |
78 | | // TODO: remove this code once fuzzer proves it is totally unreahable |
79 | | // For supported frame lengths that have odd number of factor 2 it is |
80 | | // never the last factor; consequently `ido` should never be 1. |
81 | 0 | __builtin_trap(); |
82 | | /* |
83 | | #endif |
84 | | for (k = 0; k < l1; k++) |
85 | | { |
86 | | ah = 2*k; |
87 | | ac = 4*k; |
88 | | |
89 | | RE(ch[ah]) = RE(cc[ac]) + RE(cc[ac+1]); |
90 | | RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]); |
91 | | IM(ch[ah]) = IM(cc[ac]) + IM(cc[ac+1]); |
92 | | IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]); |
93 | | } |
94 | | #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION |
95 | | */ |
96 | 0 | #endif |
97 | 180k | } else { |
98 | 361k | for (k = 0; k < l1; k++) |
99 | 180k | { |
100 | 180k | ah = k*ido; |
101 | 180k | ac = 2*k*ido; |
102 | | |
103 | 46.0M | for (i = 0; i < ido; i++) |
104 | 45.8M | { |
105 | 45.8M | complex_t t2; |
106 | | |
107 | 45.8M | RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]); |
108 | 45.8M | RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]); |
109 | | |
110 | 45.8M | IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]); |
111 | 45.8M | IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]); |
112 | | |
113 | 45.8M | #if 1 |
114 | 45.8M | ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]), |
115 | 45.8M | IM(t2), RE(t2), RE(wa[i]), IM(wa[i])); |
116 | | #else |
117 | | ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]), |
118 | | RE(t2), IM(t2), RE(wa[i]), IM(wa[i])); |
119 | | #endif |
120 | 45.8M | } |
121 | 180k | } |
122 | 180k | } |
123 | 180k | } |
124 | | |
125 | | static void passf2neg(const uint16_t ido, const uint16_t l1, const complex_t *cc, |
126 | | complex_t *ch, const complex_t *wa) |
127 | 5.67k | { |
128 | 5.67k | uint16_t i, k, ah, ac; |
129 | | |
130 | 5.67k | if (ido == 1) |
131 | 0 | { |
132 | 0 | #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION |
133 | | // TODO: remove this code once fuzzer proves it is totally unreahable |
134 | | // For supported frame lengths that have odd number of factor 2 it is |
135 | | // never the last factor; consequently `ido` should never be 1. |
136 | 0 | __builtin_trap(); |
137 | | /* |
138 | | #endif |
139 | | for (k = 0; k < l1; k++) |
140 | | { |
141 | | ah = 2*k; |
142 | | ac = 4*k; |
143 | | |
144 | | RE(ch[ah]) = RE(cc[ac]) + RE(cc[ac+1]); |
145 | | RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]); |
146 | | IM(ch[ah]) = IM(cc[ac]) + IM(cc[ac+1]); |
147 | | IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]); |
148 | | } |
149 | | #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION |
150 | | */ |
151 | 0 | #endif |
152 | 5.67k | } else { |
153 | 11.3k | for (k = 0; k < l1; k++) |
154 | 5.67k | { |
155 | 5.67k | ah = k*ido; |
156 | 5.67k | ac = 2*k*ido; |
157 | | |
158 | 1.40M | for (i = 0; i < ido; i++) |
159 | 1.40M | { |
160 | 1.40M | complex_t t2; |
161 | | |
162 | 1.40M | RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]); |
163 | 1.40M | RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]); |
164 | | |
165 | 1.40M | IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]); |
166 | 1.40M | IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]); |
167 | | |
168 | 1.40M | #if 1 |
169 | 1.40M | ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]), |
170 | 1.40M | RE(t2), IM(t2), RE(wa[i]), IM(wa[i])); |
171 | | #else |
172 | | ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]), |
173 | | IM(t2), RE(t2), RE(wa[i]), IM(wa[i])); |
174 | | #endif |
175 | 1.40M | } |
176 | 5.67k | } |
177 | 5.67k | } |
178 | 5.67k | } |
179 | | |
180 | | |
181 | | static void passf3(const uint16_t ido, const uint16_t l1, const complex_t *cc, |
182 | | complex_t *ch, const complex_t *wa1, const complex_t *wa2, |
183 | | const int8_t isign) |
184 | 45.2k | { |
185 | 45.2k | static real_t taur = FRAC_CONST(-0.5); |
186 | 45.2k | static real_t taui = FRAC_CONST(0.866025403784439); |
187 | 45.2k | uint16_t i, k, ac, ah; |
188 | 45.2k | complex_t c2, c3, d2, d3, t2; |
189 | | |
190 | 45.2k | if (ido == 1) |
191 | 0 | { |
192 | 0 | #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION |
193 | | // TODO: remove this code once fuzzer proves it is totally unreahable |
194 | | // 3 is never the the biggest factor for supported frame lengths; |
195 | | // consequently `ido` should never be 1. |
196 | 0 | __builtin_trap(); |
197 | | /* |
198 | | #endif |
199 | | if (isign == 1) |
200 | | { |
201 | | for (k = 0; k < l1; k++) |
202 | | { |
203 | | ac = 3*k+1; |
204 | | ah = k; |
205 | | |
206 | | RE(t2) = RE(cc[ac]) + RE(cc[ac+1]); |
207 | | IM(t2) = IM(cc[ac]) + IM(cc[ac+1]); |
208 | | RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur); |
209 | | IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur); |
210 | | |
211 | | RE(ch[ah]) = RE(cc[ac-1]) + RE(t2); |
212 | | IM(ch[ah]) = IM(cc[ac-1]) + IM(t2); |
213 | | |
214 | | RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui); |
215 | | IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui); |
216 | | |
217 | | RE(ch[ah+l1]) = RE(c2) - IM(c3); |
218 | | IM(ch[ah+l1]) = IM(c2) + RE(c3); |
219 | | RE(ch[ah+2*l1]) = RE(c2) + IM(c3); |
220 | | IM(ch[ah+2*l1]) = IM(c2) - RE(c3); |
221 | | } |
222 | | } else { |
223 | | for (k = 0; k < l1; k++) |
224 | | { |
225 | | ac = 3*k+1; |
226 | | ah = k; |
227 | | |
228 | | RE(t2) = RE(cc[ac]) + RE(cc[ac+1]); |
229 | | IM(t2) = IM(cc[ac]) + IM(cc[ac+1]); |
230 | | RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur); |
231 | | IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur); |
232 | | |
233 | | RE(ch[ah]) = RE(cc[ac-1]) + RE(t2); |
234 | | IM(ch[ah]) = IM(cc[ac-1]) + IM(t2); |
235 | | |
236 | | RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui); |
237 | | IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui); |
238 | | |
239 | | RE(ch[ah+l1]) = RE(c2) + IM(c3); |
240 | | IM(ch[ah+l1]) = IM(c2) - RE(c3); |
241 | | RE(ch[ah+2*l1]) = RE(c2) - IM(c3); |
242 | | IM(ch[ah+2*l1]) = IM(c2) + RE(c3); |
243 | | } |
244 | | } |
245 | | #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION |
246 | | */ |
247 | 0 | #endif |
248 | 45.2k | } else { |
249 | 45.2k | if (isign == 1) |
250 | 42.0k | { |
251 | 106k | for (k = 0; k < l1; k++) |
252 | 64.7k | { |
253 | 4.13M | for (i = 0; i < ido; i++) |
254 | 4.06M | { |
255 | 4.06M | ac = i + (3*k+1)*ido; |
256 | 4.06M | ah = i + k * ido; |
257 | | |
258 | 4.06M | RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]); |
259 | 4.06M | RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur); |
260 | 4.06M | IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]); |
261 | 4.06M | IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur); |
262 | | |
263 | 4.06M | RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2); |
264 | 4.06M | IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2); |
265 | | |
266 | 4.06M | RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui); |
267 | 4.06M | IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui); |
268 | | |
269 | 4.06M | RE(d2) = RE(c2) - IM(c3); |
270 | 4.06M | IM(d3) = IM(c2) - RE(c3); |
271 | 4.06M | RE(d3) = RE(c2) + IM(c3); |
272 | 4.06M | IM(d2) = IM(c2) + RE(c3); |
273 | | |
274 | 4.06M | #if 1 |
275 | 4.06M | ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]), |
276 | 4.06M | IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i])); |
277 | 4.06M | ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]), |
278 | 4.06M | IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i])); |
279 | | #else |
280 | | ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]), |
281 | | RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i])); |
282 | | ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]), |
283 | | RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i])); |
284 | | #endif |
285 | 4.06M | } |
286 | 64.7k | } |
287 | 42.0k | } else { |
288 | 9.51k | for (k = 0; k < l1; k++) |
289 | 6.33k | { |
290 | 512k | for (i = 0; i < ido; i++) |
291 | 506k | { |
292 | 506k | ac = i + (3*k+1)*ido; |
293 | 506k | ah = i + k * ido; |
294 | | |
295 | 506k | RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]); |
296 | 506k | RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur); |
297 | 506k | IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]); |
298 | 506k | IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur); |
299 | | |
300 | 506k | RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2); |
301 | 506k | IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2); |
302 | | |
303 | 506k | RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui); |
304 | 506k | IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui); |
305 | | |
306 | 506k | RE(d2) = RE(c2) + IM(c3); |
307 | 506k | IM(d3) = IM(c2) + RE(c3); |
308 | 506k | RE(d3) = RE(c2) - IM(c3); |
309 | 506k | IM(d2) = IM(c2) - RE(c3); |
310 | | |
311 | 506k | #if 1 |
312 | 506k | ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]), |
313 | 506k | RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i])); |
314 | 506k | ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]), |
315 | 506k | RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i])); |
316 | | #else |
317 | | ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]), |
318 | | IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i])); |
319 | | ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]), |
320 | | IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i])); |
321 | | #endif |
322 | 506k | } |
323 | 6.33k | } |
324 | 3.18k | } |
325 | 45.2k | } |
326 | 45.2k | } |
327 | | |
328 | | |
329 | | static void passf4pos(const uint16_t ido, const uint16_t l1, const complex_t *cc, |
330 | | complex_t *ch, const complex_t *wa1, const complex_t *wa2, |
331 | | const complex_t *wa3) |
332 | 911k | { |
333 | 911k | uint16_t i, k, ac, ah; |
334 | | |
335 | 911k | if (ido == 1) |
336 | 229k | { |
337 | 21.6M | for (k = 0; k < l1; k++) |
338 | 21.3M | { |
339 | 21.3M | complex_t t1, t2, t3, t4; |
340 | | |
341 | 21.3M | ac = 4*k; |
342 | 21.3M | ah = k; |
343 | | |
344 | 21.3M | RE(t2) = RE(cc[ac]) + RE(cc[ac+2]); |
345 | 21.3M | RE(t1) = RE(cc[ac]) - RE(cc[ac+2]); |
346 | 21.3M | IM(t2) = IM(cc[ac]) + IM(cc[ac+2]); |
347 | 21.3M | IM(t1) = IM(cc[ac]) - IM(cc[ac+2]); |
348 | 21.3M | RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]); |
349 | 21.3M | IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]); |
350 | 21.3M | IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]); |
351 | 21.3M | RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]); |
352 | | |
353 | 21.3M | RE(ch[ah]) = RE(t2) + RE(t3); |
354 | 21.3M | RE(ch[ah+2*l1]) = RE(t2) - RE(t3); |
355 | | |
356 | 21.3M | IM(ch[ah]) = IM(t2) + IM(t3); |
357 | 21.3M | IM(ch[ah+2*l1]) = IM(t2) - IM(t3); |
358 | | |
359 | 21.3M | RE(ch[ah+l1]) = RE(t1) + RE(t4); |
360 | 21.3M | RE(ch[ah+3*l1]) = RE(t1) - RE(t4); |
361 | | |
362 | 21.3M | IM(ch[ah+l1]) = IM(t1) + IM(t4); |
363 | 21.3M | IM(ch[ah+3*l1]) = IM(t1) - IM(t4); |
364 | 21.3M | } |
365 | 682k | } else { |
366 | 8.43M | for (k = 0; k < l1; k++) |
367 | 7.74M | { |
368 | 7.74M | ac = 4*k*ido; |
369 | 7.74M | ah = k*ido; |
370 | | |
371 | 76.6M | for (i = 0; i < ido; i++) |
372 | 68.8M | { |
373 | 68.8M | complex_t c2, c3, c4, t1, t2, t3, t4; |
374 | | |
375 | 68.8M | RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]); |
376 | 68.8M | RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]); |
377 | 68.8M | IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]); |
378 | 68.8M | IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]); |
379 | 68.8M | RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]); |
380 | 68.8M | IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]); |
381 | 68.8M | IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]); |
382 | 68.8M | RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]); |
383 | | |
384 | 68.8M | RE(c2) = RE(t1) + RE(t4); |
385 | 68.8M | RE(c4) = RE(t1) - RE(t4); |
386 | | |
387 | 68.8M | IM(c2) = IM(t1) + IM(t4); |
388 | 68.8M | IM(c4) = IM(t1) - IM(t4); |
389 | | |
390 | 68.8M | RE(ch[ah+i]) = RE(t2) + RE(t3); |
391 | 68.8M | RE(c3) = RE(t2) - RE(t3); |
392 | | |
393 | 68.8M | IM(ch[ah+i]) = IM(t2) + IM(t3); |
394 | 68.8M | IM(c3) = IM(t2) - IM(t3); |
395 | | |
396 | 68.8M | #if 1 |
397 | 68.8M | ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]), |
398 | 68.8M | IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i])); |
399 | 68.8M | ComplexMult(&IM(ch[ah+i+2*l1*ido]), &RE(ch[ah+i+2*l1*ido]), |
400 | 68.8M | IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i])); |
401 | 68.8M | ComplexMult(&IM(ch[ah+i+3*l1*ido]), &RE(ch[ah+i+3*l1*ido]), |
402 | 68.8M | IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i])); |
403 | | #else |
404 | | ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]), |
405 | | RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i])); |
406 | | ComplexMult(&RE(ch[ah+i+2*l1*ido]), &IM(ch[ah+i+2*l1*ido]), |
407 | | RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i])); |
408 | | ComplexMult(&RE(ch[ah+i+3*l1*ido]), &IM(ch[ah+i+3*l1*ido]), |
409 | | RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i])); |
410 | | #endif |
411 | 68.8M | } |
412 | 7.74M | } |
413 | 682k | } |
414 | 911k | } |
415 | | |
416 | | static void passf4neg(const uint16_t ido, const uint16_t l1, const complex_t *cc, |
417 | | complex_t *ch, const complex_t *wa1, const complex_t *wa2, |
418 | | const complex_t *wa3) |
419 | 16.8k | { |
420 | 16.8k | uint16_t i, k, ac, ah; |
421 | | |
422 | 16.8k | if (ido == 1) |
423 | 2.61k | { |
424 | 331k | for (k = 0; k < l1; k++) |
425 | 328k | { |
426 | 328k | complex_t t1, t2, t3, t4; |
427 | | |
428 | 328k | ac = 4*k; |
429 | 328k | ah = k; |
430 | | |
431 | 328k | RE(t2) = RE(cc[ac]) + RE(cc[ac+2]); |
432 | 328k | RE(t1) = RE(cc[ac]) - RE(cc[ac+2]); |
433 | 328k | IM(t2) = IM(cc[ac]) + IM(cc[ac+2]); |
434 | 328k | IM(t1) = IM(cc[ac]) - IM(cc[ac+2]); |
435 | 328k | RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]); |
436 | 328k | IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]); |
437 | 328k | IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]); |
438 | 328k | RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]); |
439 | | |
440 | 328k | RE(ch[ah]) = RE(t2) + RE(t3); |
441 | 328k | RE(ch[ah+2*l1]) = RE(t2) - RE(t3); |
442 | | |
443 | 328k | IM(ch[ah]) = IM(t2) + IM(t3); |
444 | 328k | IM(ch[ah+2*l1]) = IM(t2) - IM(t3); |
445 | | |
446 | 328k | RE(ch[ah+l1]) = RE(t1) - RE(t4); |
447 | 328k | RE(ch[ah+3*l1]) = RE(t1) + RE(t4); |
448 | | |
449 | 328k | IM(ch[ah+l1]) = IM(t1) - IM(t4); |
450 | 328k | IM(ch[ah+3*l1]) = IM(t1) + IM(t4); |
451 | 328k | } |
452 | 14.1k | } else { |
453 | 217k | for (k = 0; k < l1; k++) |
454 | 202k | { |
455 | 202k | ac = 4*k*ido; |
456 | 202k | ah = k*ido; |
457 | | |
458 | 1.94M | for (i = 0; i < ido; i++) |
459 | 1.74M | { |
460 | 1.74M | complex_t c2, c3, c4, t1, t2, t3, t4; |
461 | | |
462 | 1.74M | RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]); |
463 | 1.74M | RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]); |
464 | 1.74M | IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]); |
465 | 1.74M | IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]); |
466 | 1.74M | RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]); |
467 | 1.74M | IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]); |
468 | 1.74M | IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]); |
469 | 1.74M | RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]); |
470 | | |
471 | 1.74M | RE(c2) = RE(t1) - RE(t4); |
472 | 1.74M | RE(c4) = RE(t1) + RE(t4); |
473 | | |
474 | 1.74M | IM(c2) = IM(t1) - IM(t4); |
475 | 1.74M | IM(c4) = IM(t1) + IM(t4); |
476 | | |
477 | 1.74M | RE(ch[ah+i]) = RE(t2) + RE(t3); |
478 | 1.74M | RE(c3) = RE(t2) - RE(t3); |
479 | | |
480 | 1.74M | IM(ch[ah+i]) = IM(t2) + IM(t3); |
481 | 1.74M | IM(c3) = IM(t2) - IM(t3); |
482 | | |
483 | 1.74M | #if 1 |
484 | 1.74M | ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]), |
485 | 1.74M | RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i])); |
486 | 1.74M | ComplexMult(&RE(ch[ah+i+2*l1*ido]), &IM(ch[ah+i+2*l1*ido]), |
487 | 1.74M | RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i])); |
488 | 1.74M | ComplexMult(&RE(ch[ah+i+3*l1*ido]), &IM(ch[ah+i+3*l1*ido]), |
489 | 1.74M | RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i])); |
490 | | #else |
491 | | ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]), |
492 | | IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i])); |
493 | | ComplexMult(&IM(ch[ah+i+2*l1*ido]), &RE(ch[ah+i+2*l1*ido]), |
494 | | IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i])); |
495 | | ComplexMult(&IM(ch[ah+i+3*l1*ido]), &RE(ch[ah+i+3*l1*ido]), |
496 | | IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i])); |
497 | | #endif |
498 | 1.74M | } |
499 | 202k | } |
500 | 14.1k | } |
501 | 16.8k | } |
502 | | |
503 | | static void passf5(const uint16_t ido, const uint16_t l1, const complex_t *cc, |
504 | | complex_t *ch, const complex_t *wa1, const complex_t *wa2, const complex_t *wa3, |
505 | | const complex_t *wa4, const int8_t isign) |
506 | 45.2k | { |
507 | 45.2k | static real_t tr11 = FRAC_CONST(0.309016994374947); |
508 | 45.2k | static real_t ti11 = FRAC_CONST(0.951056516295154); |
509 | 45.2k | static real_t tr12 = FRAC_CONST(-0.809016994374947); |
510 | 45.2k | static real_t ti12 = FRAC_CONST(0.587785252292473); |
511 | 45.2k | uint16_t k, ac, ah; |
512 | 45.2k | complex_t c2, c3, c4, c5, t2, t3, t4, t5; |
513 | 45.2k | (void)wa1; |
514 | 45.2k | (void)wa2; |
515 | 45.2k | (void)wa3; |
516 | 45.2k | (void)wa4; |
517 | | |
518 | 45.2k | if (ido == 1) |
519 | 45.2k | { |
520 | 45.2k | if (isign == 1) |
521 | 42.0k | { |
522 | 2.48M | for (k = 0; k < l1; k++) |
523 | 2.44M | { |
524 | 2.44M | ac = 5*k + 1; |
525 | 2.44M | ah = k; |
526 | | |
527 | 2.44M | RE(t2) = RE(cc[ac]) + RE(cc[ac+3]); |
528 | 2.44M | IM(t2) = IM(cc[ac]) + IM(cc[ac+3]); |
529 | 2.44M | RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]); |
530 | 2.44M | IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]); |
531 | 2.44M | RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]); |
532 | 2.44M | IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]); |
533 | 2.44M | RE(t5) = RE(cc[ac]) - RE(cc[ac+3]); |
534 | 2.44M | IM(t5) = IM(cc[ac]) - IM(cc[ac+3]); |
535 | | |
536 | 2.44M | RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3); |
537 | 2.44M | IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3); |
538 | | |
539 | 2.44M | RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12); |
540 | 2.44M | IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12); |
541 | 2.44M | RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11); |
542 | 2.44M | IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11); |
543 | | |
544 | 2.44M | ComplexMult(&RE(c5), &RE(c4), |
545 | 2.44M | ti11, ti12, RE(t5), RE(t4)); |
546 | 2.44M | ComplexMult(&IM(c5), &IM(c4), |
547 | 2.44M | ti11, ti12, IM(t5), IM(t4)); |
548 | | |
549 | 2.44M | RE(ch[ah+l1]) = RE(c2) - IM(c5); |
550 | 2.44M | IM(ch[ah+l1]) = IM(c2) + RE(c5); |
551 | 2.44M | RE(ch[ah+2*l1]) = RE(c3) - IM(c4); |
552 | 2.44M | IM(ch[ah+2*l1]) = IM(c3) + RE(c4); |
553 | 2.44M | RE(ch[ah+3*l1]) = RE(c3) + IM(c4); |
554 | 2.44M | IM(ch[ah+3*l1]) = IM(c3) - RE(c4); |
555 | 2.44M | RE(ch[ah+4*l1]) = RE(c2) + IM(c5); |
556 | 2.44M | IM(ch[ah+4*l1]) = IM(c2) - RE(c5); |
557 | 2.44M | } |
558 | 42.0k | } else { |
559 | 307k | for (k = 0; k < l1; k++) |
560 | 303k | { |
561 | 303k | ac = 5*k + 1; |
562 | 303k | ah = k; |
563 | | |
564 | 303k | RE(t2) = RE(cc[ac]) + RE(cc[ac+3]); |
565 | 303k | IM(t2) = IM(cc[ac]) + IM(cc[ac+3]); |
566 | 303k | RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]); |
567 | 303k | IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]); |
568 | 303k | RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]); |
569 | 303k | IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]); |
570 | 303k | RE(t5) = RE(cc[ac]) - RE(cc[ac+3]); |
571 | 303k | IM(t5) = IM(cc[ac]) - IM(cc[ac+3]); |
572 | | |
573 | 303k | RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3); |
574 | 303k | IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3); |
575 | | |
576 | 303k | RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12); |
577 | 303k | IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12); |
578 | 303k | RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11); |
579 | 303k | IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11); |
580 | | |
581 | 303k | ComplexMult(&RE(c4), &RE(c5), |
582 | 303k | ti12, ti11, RE(t5), RE(t4)); |
583 | 303k | ComplexMult(&IM(c4), &IM(c5), |
584 | 303k | ti12, ti11, IM(t5), IM(t4)); |
585 | | |
586 | 303k | RE(ch[ah+l1]) = RE(c2) + IM(c5); |
587 | 303k | IM(ch[ah+l1]) = IM(c2) - RE(c5); |
588 | 303k | RE(ch[ah+2*l1]) = RE(c3) + IM(c4); |
589 | 303k | IM(ch[ah+2*l1]) = IM(c3) - RE(c4); |
590 | 303k | RE(ch[ah+3*l1]) = RE(c3) - IM(c4); |
591 | 303k | IM(ch[ah+3*l1]) = IM(c3) + RE(c4); |
592 | 303k | RE(ch[ah+4*l1]) = RE(c2) - IM(c5); |
593 | 303k | IM(ch[ah+4*l1]) = IM(c2) + RE(c5); |
594 | 303k | } |
595 | 3.18k | } |
596 | 45.2k | } else { |
597 | 0 | #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION |
598 | | // TODO: remove this code once fuzzer proves it is totally unreahable |
599 | | // 5 if the biggest factor and it never repeated for supported frame |
600 | | // lengths; consequently `ido` should always be 1. |
601 | 0 | __builtin_trap(); |
602 | | /* |
603 | | #else |
604 | | uint16_t i; |
605 | | complex_t d3, d4, d5, d2; |
606 | | #endif |
607 | | if (isign == 1) |
608 | | { |
609 | | for (k = 0; k < l1; k++) |
610 | | { |
611 | | for (i = 0; i < ido; i++) |
612 | | { |
613 | | ac = i + (k*5 + 1) * ido; |
614 | | ah = i + k * ido; |
615 | | |
616 | | RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]); |
617 | | IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]); |
618 | | RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]); |
619 | | IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]); |
620 | | RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]); |
621 | | IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]); |
622 | | RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]); |
623 | | IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]); |
624 | | |
625 | | RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3); |
626 | | IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3); |
627 | | |
628 | | RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12); |
629 | | IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12); |
630 | | RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11); |
631 | | IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11); |
632 | | |
633 | | ComplexMult(&RE(c5), &RE(c4), |
634 | | ti11, ti12, RE(t5), RE(t4)); |
635 | | ComplexMult(&IM(c5), &IM(c4), |
636 | | ti11, ti12, IM(t5), IM(t4)); |
637 | | |
638 | | IM(d2) = IM(c2) + RE(c5); |
639 | | IM(d3) = IM(c3) + RE(c4); |
640 | | RE(d4) = RE(c3) + IM(c4); |
641 | | RE(d5) = RE(c2) + IM(c5); |
642 | | RE(d2) = RE(c2) - IM(c5); |
643 | | IM(d5) = IM(c2) - RE(c5); |
644 | | RE(d3) = RE(c3) - IM(c4); |
645 | | IM(d4) = IM(c3) - RE(c4); |
646 | | |
647 | | #if 1 |
648 | | ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]), |
649 | | IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i])); |
650 | | ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]), |
651 | | IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i])); |
652 | | ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]), |
653 | | IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i])); |
654 | | ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]), |
655 | | IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i])); |
656 | | #else |
657 | | ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]), |
658 | | RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i])); |
659 | | ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]), |
660 | | RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i])); |
661 | | ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]), |
662 | | RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i])); |
663 | | ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]), |
664 | | RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i])); |
665 | | #endif |
666 | | } |
667 | | } |
668 | | } else { |
669 | | for (k = 0; k < l1; k++) |
670 | | { |
671 | | for (i = 0; i < ido; i++) |
672 | | { |
673 | | ac = i + (k*5 + 1) * ido; |
674 | | ah = i + k * ido; |
675 | | |
676 | | RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]); |
677 | | IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]); |
678 | | RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]); |
679 | | IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]); |
680 | | RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]); |
681 | | IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]); |
682 | | RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]); |
683 | | IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]); |
684 | | |
685 | | RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3); |
686 | | IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3); |
687 | | |
688 | | RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12); |
689 | | IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12); |
690 | | RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11); |
691 | | IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11); |
692 | | |
693 | | ComplexMult(&RE(c4), &RE(c5), |
694 | | ti12, ti11, RE(t5), RE(t4)); |
695 | | ComplexMult(&IM(c4), &IM(c5), |
696 | | ti12, ti11, IM(t5), IM(t4)); |
697 | | |
698 | | IM(d2) = IM(c2) - RE(c5); |
699 | | IM(d3) = IM(c3) - RE(c4); |
700 | | RE(d4) = RE(c3) - IM(c4); |
701 | | RE(d5) = RE(c2) - IM(c5); |
702 | | RE(d2) = RE(c2) + IM(c5); |
703 | | IM(d5) = IM(c2) + RE(c5); |
704 | | RE(d3) = RE(c3) + IM(c4); |
705 | | IM(d4) = IM(c3) + RE(c4); |
706 | | |
707 | | #if 1 |
708 | | ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]), |
709 | | RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i])); |
710 | | ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]), |
711 | | RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i])); |
712 | | ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]), |
713 | | RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i])); |
714 | | ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]), |
715 | | RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i])); |
716 | | #else |
717 | | ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]), |
718 | | IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i])); |
719 | | ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]), |
720 | | IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i])); |
721 | | ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]), |
722 | | IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i])); |
723 | | ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]), |
724 | | IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i])); |
725 | | #endif |
726 | | } |
727 | | } |
728 | | } |
729 | | #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION |
730 | | */ |
731 | 0 | #endif |
732 | 0 | } |
733 | 45.2k | } |
734 | | |
735 | | |
736 | | /*---------------------------------------------------------------------- |
737 | | cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs. |
738 | | ----------------------------------------------------------------------*/ |
739 | | |
740 | | static INLINE void cfftf1pos(uint16_t n, complex_t *c, complex_t *ch, |
741 | | const uint16_t *ifac, const complex_t *wa, |
742 | | const int8_t isign) |
743 | 271k | { |
744 | 271k | uint16_t i; |
745 | 271k | uint16_t k1, l1, l2; |
746 | 271k | uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido; |
747 | | |
748 | 271k | nf = ifac[1]; |
749 | 271k | na = 0; |
750 | 271k | l1 = 1; |
751 | 271k | iw = 0; |
752 | | |
753 | 1.44M | for (k1 = 2; k1 <= nf+1; k1++) |
754 | 1.17M | { |
755 | 1.17M | ip = ifac[k1]; |
756 | 1.17M | l2 = ip*l1; |
757 | 1.17M | ido = n / l2; |
758 | | |
759 | 1.17M | switch (ip) |
760 | 1.17M | { |
761 | 911k | case 4: |
762 | 911k | ix2 = iw + ido; |
763 | 911k | ix3 = ix2 + ido; |
764 | | |
765 | 911k | if (na == 0) |
766 | 481k | passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]); |
767 | 429k | else |
768 | 429k | passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]); |
769 | | |
770 | 911k | na = 1 - na; |
771 | 911k | break; |
772 | 180k | case 2: |
773 | 180k | if (na == 0) |
774 | 180k | passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]); |
775 | 0 | else |
776 | 0 | passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]); |
777 | | |
778 | 180k | na = 1 - na; |
779 | 180k | break; |
780 | 42.0k | case 3: |
781 | 42.0k | ix2 = iw + ido; |
782 | | |
783 | 42.0k | if (na == 0) |
784 | 19.4k | passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign); |
785 | 22.6k | else |
786 | 22.6k | passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign); |
787 | | |
788 | 42.0k | na = 1 - na; |
789 | 42.0k | break; |
790 | 42.0k | case 5: |
791 | 42.0k | ix2 = iw + ido; |
792 | 42.0k | ix3 = ix2 + ido; |
793 | 42.0k | ix4 = ix3 + ido; |
794 | | |
795 | 42.0k | if (na == 0) |
796 | 41.0k | passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); |
797 | 1.00k | else |
798 | 1.00k | passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); |
799 | | |
800 | 42.0k | na = 1 - na; |
801 | 42.0k | break; |
802 | 1.17M | } |
803 | | |
804 | 1.17M | l1 = l2; |
805 | 1.17M | iw += (ip-1) * ido; |
806 | 1.17M | } |
807 | | |
808 | 271k | if (na == 0) |
809 | 1.70k | return; |
810 | | |
811 | 97.5M | for (i = 0; i < n; i++) |
812 | 97.3M | { |
813 | 97.3M | RE(c[i]) = RE(ch[i]); |
814 | 97.3M | IM(c[i]) = IM(ch[i]); |
815 | 97.3M | } |
816 | 269k | } |
817 | | |
818 | | static INLINE void cfftf1neg(uint16_t n, complex_t *c, complex_t *ch, |
819 | | const uint16_t *ifac, const complex_t *wa, |
820 | | const int8_t isign) |
821 | 5.79k | { |
822 | 5.79k | uint16_t i; |
823 | 5.79k | uint16_t k1, l1, l2; |
824 | 5.79k | uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido; |
825 | | |
826 | 5.79k | nf = ifac[1]; |
827 | 5.79k | na = 0; |
828 | 5.79k | l1 = 1; |
829 | 5.79k | iw = 0; |
830 | | |
831 | 34.6k | for (k1 = 2; k1 <= nf+1; k1++) |
832 | 28.8k | { |
833 | 28.8k | ip = ifac[k1]; |
834 | 28.8k | l2 = ip*l1; |
835 | 28.8k | ido = n / l2; |
836 | | |
837 | 28.8k | switch (ip) |
838 | 28.8k | { |
839 | 16.8k | case 4: |
840 | 16.8k | ix2 = iw + ido; |
841 | 16.8k | ix3 = ix2 + ido; |
842 | | |
843 | 16.8k | if (na == 0) |
844 | 8.40k | passf4neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]); |
845 | 8.40k | else |
846 | 8.40k | passf4neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]); |
847 | | |
848 | 16.8k | na = 1 - na; |
849 | 16.8k | break; |
850 | 5.67k | case 2: |
851 | 5.67k | if (na == 0) |
852 | 5.67k | passf2neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]); |
853 | 0 | else |
854 | 0 | passf2neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]); |
855 | | |
856 | 5.67k | na = 1 - na; |
857 | 5.67k | break; |
858 | 3.18k | case 3: |
859 | 3.18k | ix2 = iw + ido; |
860 | | |
861 | 3.18k | if (na == 0) |
862 | 28 | passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign); |
863 | 3.15k | else |
864 | 3.15k | passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign); |
865 | | |
866 | 3.18k | na = 1 - na; |
867 | 3.18k | break; |
868 | 3.18k | case 5: |
869 | 3.18k | ix2 = iw + ido; |
870 | 3.18k | ix3 = ix2 + ido; |
871 | 3.18k | ix4 = ix3 + ido; |
872 | | |
873 | 3.18k | if (na == 0) |
874 | 3.15k | passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); |
875 | 28 | else |
876 | 28 | passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); |
877 | | |
878 | 3.18k | na = 1 - na; |
879 | 3.18k | break; |
880 | 28.8k | } |
881 | | |
882 | 28.8k | l1 = l2; |
883 | 28.8k | iw += (ip-1) * ido; |
884 | 28.8k | } |
885 | | |
886 | 5.79k | if (na == 0) |
887 | 117 | return; |
888 | | |
889 | 2.81M | for (i = 0; i < n; i++) |
890 | 2.80M | { |
891 | 2.80M | RE(c[i]) = RE(ch[i]); |
892 | 2.80M | IM(c[i]) = IM(ch[i]); |
893 | 2.80M | } |
894 | 5.67k | } |
895 | | |
896 | | void cfftf(cfft_info *cfft, complex_t *c) |
897 | 5.79k | { |
898 | 5.79k | cfftf1neg(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, -1); |
899 | 5.79k | } |
900 | | |
901 | | void cfftb(cfft_info *cfft, complex_t *c) |
902 | 271k | { |
903 | 271k | cfftf1pos(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, +1); |
904 | 271k | } |
905 | | |
906 | | static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac) |
907 | 20.4k | { |
908 | 20.4k | static uint16_t ntryh[4] = {3, 4, 2, 5}; |
909 | | #ifndef FIXED_POINT |
910 | | real_t arg, argh, argld, fi; |
911 | | uint16_t ido, ipm; |
912 | | uint16_t i1, k1, l1, l2; |
913 | | uint16_t ld, ii, ip; |
914 | | #endif |
915 | 20.4k | uint16_t ntry = 0, i, j; |
916 | 20.4k | uint16_t ib; |
917 | 20.4k | uint16_t nf, nl, nq, nr; |
918 | | |
919 | 20.4k | nl = n; |
920 | 20.4k | nf = 0; |
921 | 20.4k | j = 0; |
922 | | |
923 | 54.0k | startloop: |
924 | 54.0k | j++; |
925 | | |
926 | 54.0k | if (j <= 4) |
927 | 54.0k | ntry = ntryh[j-1]; |
928 | 0 | else |
929 | 0 | ntry += 2; |
930 | | |
931 | 54.0k | do |
932 | 115k | { |
933 | 115k | nq = nl / ntry; |
934 | 115k | nr = nl - ntry*nq; |
935 | | |
936 | 115k | if (nr != 0) |
937 | 33.6k | goto startloop; |
938 | | |
939 | 81.6k | nf++; |
940 | 81.6k | ifac[nf+1] = ntry; |
941 | 81.6k | nl = nq; |
942 | | |
943 | 81.6k | if (ntry == 2 && nf != 1) |
944 | 6.80k | { |
945 | 32.7k | for (i = 2; i <= nf; i++) |
946 | 25.9k | { |
947 | 25.9k | ib = nf - i + 2; |
948 | 25.9k | ifac[ib+1] = ifac[ib]; |
949 | 25.9k | } |
950 | 6.80k | ifac[2] = 2; |
951 | 6.80k | } |
952 | 81.6k | } while (nl != 1); |
953 | | |
954 | 20.4k | ifac[0] = n; |
955 | 20.4k | ifac[1] = nf; |
956 | | |
957 | | #ifndef FIXED_POINT |
958 | | argh = (real_t)2.0*(real_t)M_PI / (real_t)n; |
959 | | i = 0; |
960 | | l1 = 1; |
961 | | |
962 | | for (k1 = 1; k1 <= nf; k1++) |
963 | | { |
964 | | ip = ifac[k1+1]; |
965 | | ld = 0; |
966 | | l2 = l1*ip; |
967 | | ido = n / l2; |
968 | | ipm = ip - 1; |
969 | | |
970 | | for (j = 0; j < ipm; j++) |
971 | | { |
972 | | i1 = i; |
973 | | RE(wa[i]) = 1.0; |
974 | | IM(wa[i]) = 0.0; |
975 | | ld += l1; |
976 | | fi = 0; |
977 | | argld = ld*argh; |
978 | | |
979 | | for (ii = 0; ii < ido; ii++) |
980 | | { |
981 | | i++; |
982 | | fi++; |
983 | | arg = fi * argld; |
984 | | RE(wa[i]) = (real_t)cos(arg); |
985 | | #if 1 |
986 | | IM(wa[i]) = (real_t)sin(arg); |
987 | | #else |
988 | | IM(wa[i]) = (real_t)-sin(arg); |
989 | | #endif |
990 | | } |
991 | | |
992 | | if (ip > 5) |
993 | | { |
994 | | RE(wa[i1]) = RE(wa[i]); |
995 | | IM(wa[i1]) = IM(wa[i]); |
996 | | } |
997 | | } |
998 | | l1 = l2; |
999 | | } |
1000 | | #else // FIXED_POINT |
1001 | 20.4k | (void)wa; /* whoa! does it work for FIXED_POINT? */ |
1002 | 20.4k | #endif |
1003 | 20.4k | } |
1004 | | |
1005 | | cfft_info *cffti(uint16_t n) |
1006 | 20.4k | { |
1007 | 20.4k | cfft_info *cfft = (cfft_info*)faad_malloc(sizeof(cfft_info)); |
1008 | | |
1009 | 20.4k | cfft->n = n; |
1010 | 20.4k | cfft->work = (complex_t*)faad_malloc(n*sizeof(complex_t)); |
1011 | | |
1012 | | #ifndef FIXED_POINT |
1013 | | cfft->tab = (complex_t*)faad_malloc(n*sizeof(complex_t)); |
1014 | | |
1015 | | cffti1(n, cfft->tab, cfft->ifac); |
1016 | | #else |
1017 | 20.4k | cffti1(n, NULL, cfft->ifac); |
1018 | | |
1019 | 20.4k | switch (n) |
1020 | 20.4k | { |
1021 | 5.51k | case 64: cfft->tab = (complex_t*)cfft_tab_64; break; |
1022 | 5.51k | case 512: cfft->tab = (complex_t*)cfft_tab_512; break; |
1023 | 0 | #ifdef LD_DEC |
1024 | 5.51k | case 256: cfft->tab = (complex_t*)cfft_tab_256; break; |
1025 | 0 | #endif |
1026 | | |
1027 | 0 | #ifdef ALLOW_SMALL_FRAMELENGTH |
1028 | 1.29k | case 60: cfft->tab = (complex_t*)cfft_tab_60; break; |
1029 | 1.29k | case 480: cfft->tab = (complex_t*)cfft_tab_480; break; |
1030 | 0 | #ifdef LD_DEC |
1031 | 1.29k | case 240: cfft->tab = (complex_t*)cfft_tab_240; break; |
1032 | 0 | #endif |
1033 | 0 | #endif |
1034 | 0 | case 128: cfft->tab = (complex_t*)cfft_tab_128; break; |
1035 | 20.4k | } |
1036 | 20.4k | #endif |
1037 | | |
1038 | 20.4k | return cfft; |
1039 | 20.4k | } |
1040 | | |
1041 | | void cfftu(cfft_info *cfft) |
1042 | 20.4k | { |
1043 | 20.4k | if (cfft->work) faad_free(cfft->work); |
1044 | | #ifndef FIXED_POINT |
1045 | | if (cfft->tab) faad_free(cfft->tab); |
1046 | | #endif |
1047 | | |
1048 | 20.4k | if (cfft) faad_free(cfft); |
1049 | 20.4k | } |
1050 | | |