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