2 ** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding
3 ** Copyright (C) 2003-2005 M. Bakker, Nero AG, http://www.nero.com
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.
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
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"
25 ** Commercial non-GPL licensing of this software is possible.
26 ** For more info contact Nero AG through Mpeg4AAClicense@nero.com.
28 ** $Id: cfft.c,v 1.35 2007/11/01 12:33:29 menno Exp $
32 * Algorithmically based on Fortran-77 FFTPACK
33 * by Paul N. Swarztrauber(Version 4, 1985).
35 * Does even sized fft only
38 /* isign is +1 for backward and -1 for forward transforms */
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 INLINE void cfftf1(uint16_t n, complex_t *c, complex_t *ch,
64 const uint16_t *ifac, const complex_t *wa, const int8_t isign);
65 static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac);
68 /*----------------------------------------------------------------------
69 passf2, passf3, passf4, passf5. Complex FFT passes fwd and bwd.
70 ----------------------------------------------------------------------*/
72 static void passf2pos(const uint16_t ido, const uint16_t l1, const complex_t *cc,
73 complex_t *ch, const complex_t *wa)
75 uint16_t i, k, ah, ac;
79 for (k = 0; k < l1; k++)
84 RE(ch[ah]) = RE(cc[ac]) + RE(cc[ac+1]);
85 RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]);
86 IM(ch[ah]) = IM(cc[ac]) + IM(cc[ac+1]);
87 IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]);
90 for (k = 0; k < l1; k++)
95 for (i = 0; i < ido; i++)
99 RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]);
100 RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]);
102 IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]);
103 IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]);
106 ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]),
107 IM(t2), RE(t2), RE(wa[i]), IM(wa[i]));
109 ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]),
110 RE(t2), IM(t2), RE(wa[i]), IM(wa[i]));
117 static void passf2neg(const uint16_t ido, const uint16_t l1, const complex_t *cc,
118 complex_t *ch, const complex_t *wa)
120 uint16_t i, k, ah, ac;
124 for (k = 0; k < l1; k++)
129 RE(ch[ah]) = RE(cc[ac]) + RE(cc[ac+1]);
130 RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]);
131 IM(ch[ah]) = IM(cc[ac]) + IM(cc[ac+1]);
132 IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]);
135 for (k = 0; k < l1; k++)
140 for (i = 0; i < ido; i++)
144 RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]);
145 RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]);
147 IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]);
148 IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]);
151 ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]),
152 RE(t2), IM(t2), RE(wa[i]), IM(wa[i]));
154 ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]),
155 IM(t2), RE(t2), RE(wa[i]), IM(wa[i]));
163 static void passf3(const uint16_t ido, const uint16_t l1, const complex_t *cc,
164 complex_t *ch, const complex_t *wa1, const complex_t *wa2,
167 static real_t taur = FRAC_CONST(-0.5);
168 static real_t taui = FRAC_CONST(0.866025403784439);
169 uint16_t i, k, ac, ah;
170 complex_t c2, c3, d2, d3, t2;
176 for (k = 0; k < l1; k++)
181 RE(t2) = RE(cc[ac]) + RE(cc[ac+1]);
182 IM(t2) = IM(cc[ac]) + IM(cc[ac+1]);
183 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur);
184 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur);
186 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2);
187 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2);
189 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui);
190 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui);
192 RE(ch[ah+l1]) = RE(c2) - IM(c3);
193 IM(ch[ah+l1]) = IM(c2) + RE(c3);
194 RE(ch[ah+2*l1]) = RE(c2) + IM(c3);
195 IM(ch[ah+2*l1]) = IM(c2) - RE(c3);
198 for (k = 0; k < l1; k++)
203 RE(t2) = RE(cc[ac]) + RE(cc[ac+1]);
204 IM(t2) = IM(cc[ac]) + IM(cc[ac+1]);
205 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur);
206 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur);
208 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2);
209 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2);
211 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui);
212 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui);
214 RE(ch[ah+l1]) = RE(c2) + IM(c3);
215 IM(ch[ah+l1]) = IM(c2) - RE(c3);
216 RE(ch[ah+2*l1]) = RE(c2) - IM(c3);
217 IM(ch[ah+2*l1]) = IM(c2) + RE(c3);
223 for (k = 0; k < l1; k++)
225 for (i = 0; i < ido; i++)
227 ac = i + (3*k+1)*ido;
230 RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]);
231 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur);
232 IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]);
233 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur);
235 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
236 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
238 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui);
239 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui);
241 RE(d2) = RE(c2) - IM(c3);
242 IM(d3) = IM(c2) - RE(c3);
243 RE(d3) = RE(c2) + IM(c3);
244 IM(d2) = IM(c2) + RE(c3);
247 ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
248 IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
249 ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
250 IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
252 ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
253 RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
254 ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
255 RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
260 for (k = 0; k < l1; k++)
262 for (i = 0; i < ido; i++)
264 ac = i + (3*k+1)*ido;
267 RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]);
268 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur);
269 IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]);
270 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur);
272 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
273 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
275 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui);
276 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui);
278 RE(d2) = RE(c2) + IM(c3);
279 IM(d3) = IM(c2) + RE(c3);
280 RE(d3) = RE(c2) - IM(c3);
281 IM(d2) = IM(c2) - RE(c3);
284 ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
285 RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
286 ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
287 RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
289 ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
290 IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
291 ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
292 IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
301 static void passf4pos(const uint16_t ido, const uint16_t l1, const complex_t *cc,
302 complex_t *ch, const complex_t *wa1, const complex_t *wa2,
303 const complex_t *wa3)
305 uint16_t i, k, ac, ah;
309 for (k = 0; k < l1; k++)
311 complex_t t1, t2, t3, t4;
316 RE(t2) = RE(cc[ac]) + RE(cc[ac+2]);
317 RE(t1) = RE(cc[ac]) - RE(cc[ac+2]);
318 IM(t2) = IM(cc[ac]) + IM(cc[ac+2]);
319 IM(t1) = IM(cc[ac]) - IM(cc[ac+2]);
320 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]);
321 IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]);
322 IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]);
323 RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]);
325 RE(ch[ah]) = RE(t2) + RE(t3);
326 RE(ch[ah+2*l1]) = RE(t2) - RE(t3);
328 IM(ch[ah]) = IM(t2) + IM(t3);
329 IM(ch[ah+2*l1]) = IM(t2) - IM(t3);
331 RE(ch[ah+l1]) = RE(t1) + RE(t4);
332 RE(ch[ah+3*l1]) = RE(t1) - RE(t4);
334 IM(ch[ah+l1]) = IM(t1) + IM(t4);
335 IM(ch[ah+3*l1]) = IM(t1) - IM(t4);
338 for (k = 0; k < l1; k++)
343 for (i = 0; i < ido; i++)
345 complex_t c2, c3, c4, t1, t2, t3, t4;
347 RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]);
348 RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]);
349 IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]);
350 IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]);
351 RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]);
352 IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]);
353 IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]);
354 RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]);
356 RE(c2) = RE(t1) + RE(t4);
357 RE(c4) = RE(t1) - RE(t4);
359 IM(c2) = IM(t1) + IM(t4);
360 IM(c4) = IM(t1) - IM(t4);
362 RE(ch[ah+i]) = RE(t2) + RE(t3);
363 RE(c3) = RE(t2) - RE(t3);
365 IM(ch[ah+i]) = IM(t2) + IM(t3);
366 IM(c3) = IM(t2) - IM(t3);
369 ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]),
370 IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i]));
371 ComplexMult(&IM(ch[ah+i+2*l1*ido]), &RE(ch[ah+i+2*l1*ido]),
372 IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i]));
373 ComplexMult(&IM(ch[ah+i+3*l1*ido]), &RE(ch[ah+i+3*l1*ido]),
374 IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i]));
376 ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]),
377 RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i]));
378 ComplexMult(&RE(ch[ah+i+2*l1*ido]), &IM(ch[ah+i+2*l1*ido]),
379 RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i]));
380 ComplexMult(&RE(ch[ah+i+3*l1*ido]), &IM(ch[ah+i+3*l1*ido]),
381 RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i]));
388 static void passf4neg(const uint16_t ido, const uint16_t l1, const complex_t *cc,
389 complex_t *ch, const complex_t *wa1, const complex_t *wa2,
390 const complex_t *wa3)
392 uint16_t i, k, ac, ah;
396 for (k = 0; k < l1; k++)
398 complex_t t1, t2, t3, t4;
403 RE(t2) = RE(cc[ac]) + RE(cc[ac+2]);
404 RE(t1) = RE(cc[ac]) - RE(cc[ac+2]);
405 IM(t2) = IM(cc[ac]) + IM(cc[ac+2]);
406 IM(t1) = IM(cc[ac]) - IM(cc[ac+2]);
407 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]);
408 IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]);
409 IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]);
410 RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]);
412 RE(ch[ah]) = RE(t2) + RE(t3);
413 RE(ch[ah+2*l1]) = RE(t2) - RE(t3);
415 IM(ch[ah]) = IM(t2) + IM(t3);
416 IM(ch[ah+2*l1]) = IM(t2) - IM(t3);
418 RE(ch[ah+l1]) = RE(t1) - RE(t4);
419 RE(ch[ah+3*l1]) = RE(t1) + RE(t4);
421 IM(ch[ah+l1]) = IM(t1) - IM(t4);
422 IM(ch[ah+3*l1]) = IM(t1) + IM(t4);
425 for (k = 0; k < l1; k++)
430 for (i = 0; i < ido; i++)
432 complex_t c2, c3, c4, t1, t2, t3, t4;
434 RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]);
435 RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]);
436 IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]);
437 IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]);
438 RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]);
439 IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]);
440 IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]);
441 RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]);
443 RE(c2) = RE(t1) - RE(t4);
444 RE(c4) = RE(t1) + RE(t4);
446 IM(c2) = IM(t1) - IM(t4);
447 IM(c4) = IM(t1) + IM(t4);
449 RE(ch[ah+i]) = RE(t2) + RE(t3);
450 RE(c3) = RE(t2) - RE(t3);
452 IM(ch[ah+i]) = IM(t2) + IM(t3);
453 IM(c3) = IM(t2) - IM(t3);
456 ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]),
457 RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i]));
458 ComplexMult(&RE(ch[ah+i+2*l1*ido]), &IM(ch[ah+i+2*l1*ido]),
459 RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i]));
460 ComplexMult(&RE(ch[ah+i+3*l1*ido]), &IM(ch[ah+i+3*l1*ido]),
461 RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i]));
463 ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]),
464 IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i]));
465 ComplexMult(&IM(ch[ah+i+2*l1*ido]), &RE(ch[ah+i+2*l1*ido]),
466 IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i]));
467 ComplexMult(&IM(ch[ah+i+3*l1*ido]), &RE(ch[ah+i+3*l1*ido]),
468 IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i]));
475 static void passf5(const uint16_t ido, const uint16_t l1, const complex_t *cc,
476 complex_t *ch, const complex_t *wa1, const complex_t *wa2, const complex_t *wa3,
477 const complex_t *wa4, const int8_t isign)
479 static real_t tr11 = FRAC_CONST(0.309016994374947);
480 static real_t ti11 = FRAC_CONST(0.951056516295154);
481 static real_t tr12 = FRAC_CONST(-0.809016994374947);
482 static real_t ti12 = FRAC_CONST(0.587785252292473);
483 uint16_t i, k, ac, ah;
484 complex_t c2, c3, c4, c5, d3, d4, d5, d2, t2, t3, t4, t5;
490 for (k = 0; k < l1; k++)
495 RE(t2) = RE(cc[ac]) + RE(cc[ac+3]);
496 IM(t2) = IM(cc[ac]) + IM(cc[ac+3]);
497 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]);
498 IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]);
499 RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]);
500 IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]);
501 RE(t5) = RE(cc[ac]) - RE(cc[ac+3]);
502 IM(t5) = IM(cc[ac]) - IM(cc[ac+3]);
504 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
505 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
507 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
508 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
509 RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
510 IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
512 ComplexMult(&RE(c5), &RE(c4),
513 ti11, ti12, RE(t5), RE(t4));
514 ComplexMult(&IM(c5), &IM(c4),
515 ti11, ti12, IM(t5), IM(t4));
517 RE(ch[ah+l1]) = RE(c2) - IM(c5);
518 IM(ch[ah+l1]) = IM(c2) + RE(c5);
519 RE(ch[ah+2*l1]) = RE(c3) - IM(c4);
520 IM(ch[ah+2*l1]) = IM(c3) + RE(c4);
521 RE(ch[ah+3*l1]) = RE(c3) + IM(c4);
522 IM(ch[ah+3*l1]) = IM(c3) - RE(c4);
523 RE(ch[ah+4*l1]) = RE(c2) + IM(c5);
524 IM(ch[ah+4*l1]) = IM(c2) - RE(c5);
527 for (k = 0; k < l1; k++)
532 RE(t2) = RE(cc[ac]) + RE(cc[ac+3]);
533 IM(t2) = IM(cc[ac]) + IM(cc[ac+3]);
534 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]);
535 IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]);
536 RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]);
537 IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]);
538 RE(t5) = RE(cc[ac]) - RE(cc[ac+3]);
539 IM(t5) = IM(cc[ac]) - IM(cc[ac+3]);
541 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
542 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
544 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
545 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
546 RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
547 IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
549 ComplexMult(&RE(c4), &RE(c5),
550 ti12, ti11, RE(t5), RE(t4));
551 ComplexMult(&IM(c4), &IM(c5),
552 ti12, ti11, IM(t5), IM(t4));
554 RE(ch[ah+l1]) = RE(c2) + IM(c5);
555 IM(ch[ah+l1]) = IM(c2) - RE(c5);
556 RE(ch[ah+2*l1]) = RE(c3) + IM(c4);
557 IM(ch[ah+2*l1]) = IM(c3) - RE(c4);
558 RE(ch[ah+3*l1]) = RE(c3) - IM(c4);
559 IM(ch[ah+3*l1]) = IM(c3) + RE(c4);
560 RE(ch[ah+4*l1]) = RE(c2) - IM(c5);
561 IM(ch[ah+4*l1]) = IM(c2) + RE(c5);
567 for (k = 0; k < l1; k++)
569 for (i = 0; i < ido; i++)
571 ac = i + (k*5 + 1) * ido;
574 RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]);
575 IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]);
576 RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]);
577 IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]);
578 RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]);
579 IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]);
580 RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]);
581 IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]);
583 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3);
584 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3);
586 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
587 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
588 RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
589 IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
591 ComplexMult(&RE(c5), &RE(c4),
592 ti11, ti12, RE(t5), RE(t4));
593 ComplexMult(&IM(c5), &IM(c4),
594 ti11, ti12, IM(t5), IM(t4));
596 IM(d2) = IM(c2) + RE(c5);
597 IM(d3) = IM(c3) + RE(c4);
598 RE(d4) = RE(c3) + IM(c4);
599 RE(d5) = RE(c2) + IM(c5);
600 RE(d2) = RE(c2) - IM(c5);
601 IM(d5) = IM(c2) - RE(c5);
602 RE(d3) = RE(c3) - IM(c4);
603 IM(d4) = IM(c3) - RE(c4);
606 ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
607 IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
608 ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
609 IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
610 ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]),
611 IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i]));
612 ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]),
613 IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i]));
615 ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
616 RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
617 ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
618 RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
619 ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]),
620 RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i]));
621 ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]),
622 RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i]));
627 for (k = 0; k < l1; k++)
629 for (i = 0; i < ido; i++)
631 ac = i + (k*5 + 1) * ido;
634 RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]);
635 IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]);
636 RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]);
637 IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]);
638 RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]);
639 IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]);
640 RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]);
641 IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]);
643 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3);
644 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3);
646 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
647 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
648 RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
649 IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
651 ComplexMult(&RE(c4), &RE(c5),
652 ti12, ti11, RE(t5), RE(t4));
653 ComplexMult(&IM(c4), &IM(c5),
654 ti12, ti11, IM(t5), IM(t4));
656 IM(d2) = IM(c2) - RE(c5);
657 IM(d3) = IM(c3) - RE(c4);
658 RE(d4) = RE(c3) - IM(c4);
659 RE(d5) = RE(c2) - IM(c5);
660 RE(d2) = RE(c2) + IM(c5);
661 IM(d5) = IM(c2) + RE(c5);
662 RE(d3) = RE(c3) + IM(c4);
663 IM(d4) = IM(c3) + RE(c4);
666 ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
667 RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
668 ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
669 RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
670 ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]),
671 RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i]));
672 ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]),
673 RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i]));
675 ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
676 IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
677 ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
678 IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
679 ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]),
680 IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i]));
681 ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]),
682 IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i]));
691 /*----------------------------------------------------------------------
692 cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
693 ----------------------------------------------------------------------*/
695 static INLINE void cfftf1pos(uint16_t n, complex_t *c, complex_t *ch,
696 const uint16_t *ifac, const complex_t *wa,
701 uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
708 for (k1 = 2; k1 <= nf+1; k1++)
722 passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]);
724 passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]);
730 passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]);
732 passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]);
740 passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign);
742 passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign);
752 passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
754 passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
767 for (i = 0; i < n; i++)
769 RE(c[i]) = RE(ch[i]);
770 IM(c[i]) = IM(ch[i]);
774 static INLINE void cfftf1neg(uint16_t n, complex_t *c, complex_t *ch,
775 const uint16_t *ifac, const complex_t *wa,
780 uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
787 for (k1 = 2; k1 <= nf+1; k1++)
801 passf4neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]);
803 passf4neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]);
809 passf2neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]);
811 passf2neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]);
819 passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign);
821 passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign);
831 passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
833 passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
846 for (i = 0; i < n; i++)
848 RE(c[i]) = RE(ch[i]);
849 IM(c[i]) = IM(ch[i]);
853 void cfftf(cfft_info *cfft, complex_t *c)
855 cfftf1neg(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, -1);
858 void cfftb(cfft_info *cfft, complex_t *c)
860 cfftf1pos(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, +1);
863 static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac)
865 static uint16_t ntryh[4] = {3, 4, 2, 5};
867 real_t arg, argh, argld, fi;
869 uint16_t i1, k1, l1, l2;
872 uint16_t ntry = 0, i, j;
874 uint16_t nf, nl, nq, nr;
900 if (ntry == 2 && nf != 1)
902 for (i = 2; i <= nf; i++)
905 ifac[ib+1] = ifac[ib];
915 argh = (real_t)2.0*(real_t)M_PI / (real_t)n;
919 for (k1 = 1; k1 <= nf; k1++)
927 for (j = 0; j < ipm; j++)
936 for (ii = 0; ii < ido; ii++)
941 RE(wa[i]) = (real_t)cos(arg);
943 IM(wa[i]) = (real_t)sin(arg);
945 IM(wa[i]) = (real_t)-sin(arg);
951 RE(wa[i1]) = RE(wa[i]);
952 IM(wa[i1]) = IM(wa[i]);
960 cfft_info *cffti(uint16_t n)
962 cfft_info *cfft = (cfft_info*)faad_malloc(sizeof(cfft_info));
965 cfft->work = (complex_t*)faad_malloc(n*sizeof(complex_t));
968 cfft->tab = (complex_t*)faad_malloc(n*sizeof(complex_t));
970 cffti1(n, cfft->tab, cfft->ifac);
972 cffti1(n, NULL, cfft->ifac);
976 case 64: cfft->tab = (complex_t*)cfft_tab_64; break;
977 case 512: cfft->tab = (complex_t*)cfft_tab_512; break;
979 case 256: cfft->tab = (complex_t*)cfft_tab_256; break;
982 #ifdef ALLOW_SMALL_FRAMELENGTH
983 case 60: cfft->tab = (complex_t*)cfft_tab_60; break;
984 case 480: cfft->tab = (complex_t*)cfft_tab_480; break;
986 case 240: cfft->tab = (complex_t*)cfft_tab_240; break;
989 case 128: cfft->tab = (complex_t*)cfft_tab_128; break;
996 void cfftu(cfft_info *cfft)
998 if (cfft->work) faad_free(cfft->work);
1000 if (cfft->tab) faad_free(cfft->tab);
1003 if (cfft) faad_free(cfft);