sync with ffmpeg
[mplayer/glamo.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.
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.
19 ** Any non-GPL usage of this software or parts of this software is strictly
20 ** forbidden.
22 ** Commercial non-GPL licensing of this software is possible.
23 ** For more info contact Ahead Software through Mpeg4AAClicense@nero.com.
25 ** $Id: cfft.c,v 1.30 2004/09/08 09:43:11 gcp Exp $
26 **/
29 * Algorithmically based on Fortran-77 FFTPACK
30 * by Paul N. Swarztrauber(Version 4, 1985).
32 * Does even sized fft only
35 /* isign is +1 for backward and -1 for forward transforms */
37 #include "common.h"
38 #include "structs.h"
40 #include <stdlib.h>
42 #include "cfft.h"
43 #include "cfft_tab.h"
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);
65 /*----------------------------------------------------------------------
66 passf2, passf3, passf4, passf5. Complex FFT passes fwd and bwd.
67 ----------------------------------------------------------------------*/
69 static void passf2pos(const uint16_t ido, const uint16_t l1, const complex_t *cc,
70 complex_t *ch, const complex_t *wa)
72 uint16_t i, k, ah, ac;
74 if (ido == 1)
76 for (k = 0; k < l1; k++)
78 ah = 2*k;
79 ac = 4*k;
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]);
86 } else {
87 for (k = 0; k < l1; k++)
89 ah = k*ido;
90 ac = 2*k*ido;
92 for (i = 0; i < ido; i++)
94 complex_t t2;
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]);
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]);
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
114 static void passf2neg(const uint16_t ido, const uint16_t l1, const complex_t *cc,
115 complex_t *ch, const complex_t *wa)
117 uint16_t i, k, ah, ac;
119 if (ido == 1)
121 for (k = 0; k < l1; k++)
123 ah = 2*k;
124 ac = 4*k;
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]);
131 } else {
132 for (k = 0; k < l1; k++)
134 ah = k*ido;
135 ac = 2*k*ido;
137 for (i = 0; i < ido; i++)
139 complex_t t2;
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]);
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]);
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
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)
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;
169 if (ido == 1)
171 if (isign == 1)
173 for (k = 0; k < l1; k++)
175 ac = 3*k+1;
176 ah = k;
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);
183 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2);
184 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2);
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);
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);
194 } else {
195 for (k = 0; k < l1; k++)
197 ac = 3*k+1;
198 ah = k;
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);
205 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2);
206 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2);
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);
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);
217 } else {
218 if (isign == 1)
220 for (k = 0; k < l1; k++)
222 for (i = 0; i < ido; i++)
224 ac = i + (3*k+1)*ido;
225 ah = i + k * ido;
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);
232 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
233 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
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);
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);
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
256 } else {
257 for (k = 0; k < l1; k++)
259 for (i = 0; i < ido; i++)
261 ac = i + (3*k+1)*ido;
262 ah = i + k * ido;
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);
269 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
270 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
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);
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);
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
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)
302 uint16_t i, k, ac, ah;
304 if (ido == 1)
306 for (k = 0; k < l1; k++)
308 complex_t t1, t2, t3, t4;
310 ac = 4*k;
311 ah = k;
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]);
322 RE(ch[ah]) = RE(t2) + RE(t3);
323 RE(ch[ah+2*l1]) = RE(t2) - RE(t3);
325 IM(ch[ah]) = IM(t2) + IM(t3);
326 IM(ch[ah+2*l1]) = IM(t2) - IM(t3);
328 RE(ch[ah+l1]) = RE(t1) + RE(t4);
329 RE(ch[ah+3*l1]) = RE(t1) - RE(t4);
331 IM(ch[ah+l1]) = IM(t1) + IM(t4);
332 IM(ch[ah+3*l1]) = IM(t1) - IM(t4);
334 } else {
335 for (k = 0; k < l1; k++)
337 ac = 4*k*ido;
338 ah = k*ido;
340 for (i = 0; i < ido; i++)
342 complex_t c2, c3, c4, t1, t2, t3, t4;
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]);
353 RE(c2) = RE(t1) + RE(t4);
354 RE(c4) = RE(t1) - RE(t4);
356 IM(c2) = IM(t1) + IM(t4);
357 IM(c4) = IM(t1) - IM(t4);
359 RE(ch[ah+i]) = RE(t2) + RE(t3);
360 RE(c3) = RE(t2) - RE(t3);
362 IM(ch[ah+i]) = IM(t2) + IM(t3);
363 IM(c3) = IM(t2) - IM(t3);
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
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)
389 uint16_t i, k, ac, ah;
391 if (ido == 1)
393 for (k = 0; k < l1; k++)
395 complex_t t1, t2, t3, t4;
397 ac = 4*k;
398 ah = k;
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]);
409 RE(ch[ah]) = RE(t2) + RE(t3);
410 RE(ch[ah+2*l1]) = RE(t2) - RE(t3);
412 IM(ch[ah]) = IM(t2) + IM(t3);
413 IM(ch[ah+2*l1]) = IM(t2) - IM(t3);
415 RE(ch[ah+l1]) = RE(t1) - RE(t4);
416 RE(ch[ah+3*l1]) = RE(t1) + RE(t4);
418 IM(ch[ah+l1]) = IM(t1) - IM(t4);
419 IM(ch[ah+3*l1]) = IM(t1) + IM(t4);
421 } else {
422 for (k = 0; k < l1; k++)
424 ac = 4*k*ido;
425 ah = k*ido;
427 for (i = 0; i < ido; i++)
429 complex_t c2, c3, c4, t1, t2, t3, t4;
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]);
440 RE(c2) = RE(t1) - RE(t4);
441 RE(c4) = RE(t1) + RE(t4);
443 IM(c2) = IM(t1) - IM(t4);
444 IM(c4) = IM(t1) + IM(t4);
446 RE(ch[ah+i]) = RE(t2) + RE(t3);
447 RE(c3) = RE(t2) - RE(t3);
449 IM(ch[ah+i]) = IM(t2) + IM(t3);
450 IM(c3) = IM(t2) - IM(t3);
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
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)
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;
483 if (ido == 1)
485 if (isign == 1)
487 for (k = 0; k < l1; k++)
489 ac = 5*k + 1;
490 ah = k;
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]);
501 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
502 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
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);
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));
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);
523 } else {
524 for (k = 0; k < l1; k++)
526 ac = 5*k + 1;
527 ah = k;
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]);
538 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
539 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
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);
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));
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);
561 } else {
562 if (isign == 1)
564 for (k = 0; k < l1; k++)
566 for (i = 0; i < ido; i++)
568 ac = i + (k*5 + 1) * ido;
569 ah = i + k * ido;
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]);
580 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3);
581 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3);
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);
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));
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);
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
623 } else {
624 for (k = 0; k < l1; k++)
626 for (i = 0; i < ido; i++)
628 ac = i + (k*5 + 1) * ido;
629 ah = i + k * ido;
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]);
640 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3);
641 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3);
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);
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));
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);
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
688 /*----------------------------------------------------------------------
689 cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
690 ----------------------------------------------------------------------*/
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)
696 uint16_t i;
697 uint16_t k1, l1, l2;
698 uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
700 nf = ifac[1];
701 na = 0;
702 l1 = 1;
703 iw = 0;
705 for (k1 = 2; k1 <= nf+1; k1++)
707 ip = ifac[k1];
708 l2 = ip*l1;
709 ido = n / l2;
710 idl1 = ido*l1;
712 switch (ip)
714 case 4:
715 ix2 = iw + ido;
716 ix3 = ix2 + ido;
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]);
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]);
731 na = 1 - na;
732 break;
733 case 3:
734 ix2 = iw + ido;
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);
741 na = 1 - na;
742 break;
743 case 5:
744 ix2 = iw + ido;
745 ix3 = ix2 + ido;
746 ix4 = ix3 + ido;
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);
753 na = 1 - na;
754 break;
757 l1 = l2;
758 iw += (ip-1) * ido;
761 if (na == 0)
762 return;
764 for (i = 0; i < n; i++)
766 RE(c[i]) = RE(ch[i]);
767 IM(c[i]) = IM(ch[i]);
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)
775 uint16_t i;
776 uint16_t k1, l1, l2;
777 uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
779 nf = ifac[1];
780 na = 0;
781 l1 = 1;
782 iw = 0;
784 for (k1 = 2; k1 <= nf+1; k1++)
786 ip = ifac[k1];
787 l2 = ip*l1;
788 ido = n / l2;
789 idl1 = ido*l1;
791 switch (ip)
793 case 4:
794 ix2 = iw + ido;
795 ix3 = ix2 + ido;
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]);
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]);
810 na = 1 - na;
811 break;
812 case 3:
813 ix2 = iw + ido;
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);
820 na = 1 - na;
821 break;
822 case 5:
823 ix2 = iw + ido;
824 ix3 = ix2 + ido;
825 ix4 = ix3 + ido;
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);
832 na = 1 - na;
833 break;
836 l1 = l2;
837 iw += (ip-1) * ido;
840 if (na == 0)
841 return;
843 for (i = 0; i < n; i++)
845 RE(c[i]) = RE(ch[i]);
846 IM(c[i]) = IM(ch[i]);
850 void cfftf(cfft_info *cfft, complex_t *c)
852 cfftf1neg(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, -1);
855 void cfftb(cfft_info *cfft, complex_t *c)
857 cfftf1pos(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, +1);
860 static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac)
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;
873 nl = n;
874 nf = 0;
875 j = 0;
877 startloop:
878 j++;
880 if (j <= 4)
881 ntry = ntryh[j-1];
882 else
883 ntry += 2;
887 nq = nl / ntry;
888 nr = nl - ntry*nq;
890 if (nr != 0)
891 goto startloop;
893 nf++;
894 ifac[nf+1] = ntry;
895 nl = nq;
897 if (ntry == 2 && nf != 1)
899 for (i = 2; i <= nf; i++)
901 ib = nf - i + 2;
902 ifac[ib+1] = ifac[ib];
904 ifac[2] = 2;
906 } while (nl != 1);
908 ifac[0] = n;
909 ifac[1] = nf;
911 #ifndef FIXED_POINT
912 argh = (real_t)2.0*(real_t)M_PI / (real_t)n;
913 i = 0;
914 l1 = 1;
916 for (k1 = 1; k1 <= nf; k1++)
918 ip = ifac[k1+1];
919 ld = 0;
920 l2 = l1*ip;
921 ido = n / l2;
922 ipm = ip - 1;
924 for (j = 0; j < ipm; j++)
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;
933 for (ii = 0; ii < ido; ii++)
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
946 if (ip > 5)
948 RE(wa[i1]) = RE(wa[i]);
949 IM(wa[i1]) = IM(wa[i]);
952 l1 = l2;
954 #endif
957 cfft_info *cffti(uint16_t n)
959 cfft_info *cfft = (cfft_info*)faad_malloc(sizeof(cfft_info));
961 cfft->n = n;
962 cfft->work = (complex_t*)faad_malloc(n*sizeof(complex_t));
964 #ifndef FIXED_POINT
965 cfft->tab = (complex_t*)faad_malloc(n*sizeof(complex_t));
967 cffti1(n, cfft->tab, cfft->ifac);
968 #else
969 cffti1(n, NULL, cfft->ifac);
971 switch (n)
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
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;
988 #endif
990 return cfft;
993 void cfftu(cfft_info *cfft)
995 if (cfft->work) faad_free(cfft->work);
996 #ifndef FIXED_POINT
997 if (cfft->tab) faad_free(cfft->tab);
998 #endif
1000 if (cfft) faad_free(cfft);