diff options
Diffstat (limited to 'faad2/src/libfaad/cfft.c')
| -rw-r--r-- | faad2/src/libfaad/cfft.c | 1005 |
1 files changed, 1005 insertions, 0 deletions
diff --git a/faad2/src/libfaad/cfft.c b/faad2/src/libfaad/cfft.c new file mode 100644 index 0000000..4235c11 --- /dev/null +++ b/faad2/src/libfaad/cfft.c @@ -0,0 +1,1005 @@ +/*
+** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding
+** Copyright (C) 2003-2005 M. Bakker, Nero AG, http://www.nero.com
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License as published by
+** the Free Software Foundation; either version 2 of the License, or
+** (at your option) any later version.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+**
+** Any non-GPL usage of this software or parts of this software is strictly
+** forbidden.
+**
+** The "appropriate copyright message" mentioned in section 2c of the GPLv2
+** must read: "Code from FAAD2 is copyright (c) Nero AG, www.nero.com"
+**
+** Commercial non-GPL licensing of this software is possible.
+** For more info contact Nero AG through Mpeg4AAClicense@nero.com.
+**
+** $Id: cfft.c,v 1.35 2007/11/01 12:33:29 menno Exp $
+**/
+
+/*
+ * Algorithmically based on Fortran-77 FFTPACK
+ * by Paul N. Swarztrauber(Version 4, 1985).
+ *
+ * Does even sized fft only
+ */
+
+/* isign is +1 for backward and -1 for forward transforms */
+
+#include "common.h"
+#include "structs.h"
+
+#include <stdlib.h>
+
+#include "cfft.h"
+#include "cfft_tab.h"
+
+
+/* static function declarations */
+static void passf2pos(const uint16_t ido, const uint16_t l1, const complex_t *cc,
+ complex_t *ch, const complex_t *wa);
+static void passf2neg(const uint16_t ido, const uint16_t l1, const complex_t *cc,
+ complex_t *ch, const complex_t *wa);
+static void passf3(const uint16_t ido, const uint16_t l1, const complex_t *cc,
+ complex_t *ch, const complex_t *wa1, const complex_t *wa2, const int8_t isign);
+static void passf4pos(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch,
+ const complex_t *wa1, const complex_t *wa2, const complex_t *wa3);
+static void passf4neg(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch,
+ const complex_t *wa1, const complex_t *wa2, const complex_t *wa3);
+static void passf5(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch,
+ const complex_t *wa1, const complex_t *wa2, const complex_t *wa3,
+ const complex_t *wa4, const int8_t isign);
+INLINE void cfftf1(uint16_t n, complex_t *c, complex_t *ch,
+ const uint16_t *ifac, const complex_t *wa, const int8_t isign);
+static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac);
+
+
+/*----------------------------------------------------------------------
+ passf2, passf3, passf4, passf5. Complex FFT passes fwd and bwd.
+ ----------------------------------------------------------------------*/
+
+static void passf2pos(const uint16_t ido, const uint16_t l1, const complex_t *cc,
+ complex_t *ch, const complex_t *wa)
+{
+ uint16_t i, k, ah, ac;
+
+ if (ido == 1)
+ {
+ for (k = 0; k < l1; k++)
+ {
+ ah = 2*k;
+ ac = 4*k;
+
+ RE(ch[ah]) = RE(cc[ac]) + RE(cc[ac+1]);
+ RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]);
+ IM(ch[ah]) = IM(cc[ac]) + IM(cc[ac+1]);
+ IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]);
+ }
+ } else {
+ for (k = 0; k < l1; k++)
+ {
+ ah = k*ido;
+ ac = 2*k*ido;
+
+ for (i = 0; i < ido; i++)
+ {
+ complex_t t2;
+
+ RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]);
+ RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]);
+
+ IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]);
+ IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]);
+
+#if 1
+ ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]),
+ IM(t2), RE(t2), RE(wa[i]), IM(wa[i]));
+#else
+ ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]),
+ RE(t2), IM(t2), RE(wa[i]), IM(wa[i]));
+#endif
+ }
+ }
+ }
+}
+
+static void passf2neg(const uint16_t ido, const uint16_t l1, const complex_t *cc,
+ complex_t *ch, const complex_t *wa)
+{
+ uint16_t i, k, ah, ac;
+
+ if (ido == 1)
+ {
+ for (k = 0; k < l1; k++)
+ {
+ ah = 2*k;
+ ac = 4*k;
+
+ RE(ch[ah]) = RE(cc[ac]) + RE(cc[ac+1]);
+ RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]);
+ IM(ch[ah]) = IM(cc[ac]) + IM(cc[ac+1]);
+ IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]);
+ }
+ } else {
+ for (k = 0; k < l1; k++)
+ {
+ ah = k*ido;
+ ac = 2*k*ido;
+
+ for (i = 0; i < ido; i++)
+ {
+ complex_t t2;
+
+ RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]);
+ RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]);
+
+ IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]);
+ IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]);
+
+#if 1
+ ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]),
+ RE(t2), IM(t2), RE(wa[i]), IM(wa[i]));
+#else
+ ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]),
+ IM(t2), RE(t2), RE(wa[i]), IM(wa[i]));
+#endif
+ }
+ }
+ }
+}
+
+
+static void passf3(const uint16_t ido, const uint16_t l1, const complex_t *cc,
+ complex_t *ch, const complex_t *wa1, const complex_t *wa2,
+ const int8_t isign)
+{
+ static real_t taur = FRAC_CONST(-0.5);
+ static real_t taui = FRAC_CONST(0.866025403784439);
+ uint16_t i, k, ac, ah;
+ complex_t c2, c3, d2, d3, t2;
+
+ if (ido == 1)
+ {
+ if (isign == 1)
+ {
+ for (k = 0; k < l1; k++)
+ {
+ ac = 3*k+1;
+ ah = k;
+
+ RE(t2) = RE(cc[ac]) + RE(cc[ac+1]);
+ IM(t2) = IM(cc[ac]) + IM(cc[ac+1]);
+ RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur);
+ IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur);
+
+ RE(ch[ah]) = RE(cc[ac-1]) + RE(t2);
+ IM(ch[ah]) = IM(cc[ac-1]) + IM(t2);
+
+ RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui);
+ IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui);
+
+ RE(ch[ah+l1]) = RE(c2) - IM(c3);
+ IM(ch[ah+l1]) = IM(c2) + RE(c3);
+ RE(ch[ah+2*l1]) = RE(c2) + IM(c3);
+ IM(ch[ah+2*l1]) = IM(c2) - RE(c3);
+ }
+ } else {
+ for (k = 0; k < l1; k++)
+ {
+ ac = 3*k+1;
+ ah = k;
+
+ RE(t2) = RE(cc[ac]) + RE(cc[ac+1]);
+ IM(t2) = IM(cc[ac]) + IM(cc[ac+1]);
+ RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur);
+ IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur);
+
+ RE(ch[ah]) = RE(cc[ac-1]) + RE(t2);
+ IM(ch[ah]) = IM(cc[ac-1]) + IM(t2);
+
+ RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui);
+ IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui);
+
+ RE(ch[ah+l1]) = RE(c2) + IM(c3);
+ IM(ch[ah+l1]) = IM(c2) - RE(c3);
+ RE(ch[ah+2*l1]) = RE(c2) - IM(c3);
+ IM(ch[ah+2*l1]) = IM(c2) + RE(c3);
+ }
+ }
+ } else {
+ if (isign == 1)
+ {
+ for (k = 0; k < l1; k++)
+ {
+ for (i = 0; i < ido; i++)
+ {
+ ac = i + (3*k+1)*ido;
+ ah = i + k * ido;
+
+ RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]);
+ RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur);
+ IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]);
+ IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur);
+
+ RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
+ IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
+
+ RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui);
+ IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui);
+
+ RE(d2) = RE(c2) - IM(c3);
+ IM(d3) = IM(c2) - RE(c3);
+ RE(d3) = RE(c2) + IM(c3);
+ IM(d2) = IM(c2) + RE(c3);
+
+#if 1
+ ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
+ IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
+ ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
+ IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
+#else
+ ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
+ RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
+ ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
+ RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
+#endif
+ }
+ }
+ } else {
+ for (k = 0; k < l1; k++)
+ {
+ for (i = 0; i < ido; i++)
+ {
+ ac = i + (3*k+1)*ido;
+ ah = i + k * ido;
+
+ RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]);
+ RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur);
+ IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]);
+ IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur);
+
+ RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
+ IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
+
+ RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui);
+ IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui);
+
+ RE(d2) = RE(c2) + IM(c3);
+ IM(d3) = IM(c2) + RE(c3);
+ RE(d3) = RE(c2) - IM(c3);
+ IM(d2) = IM(c2) - RE(c3);
+
+#if 1
+ ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
+ RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
+ ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
+ RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
+#else
+ ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
+ IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
+ ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
+ IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
+#endif
+ }
+ }
+ }
+ }
+}
+
+
+static void passf4pos(const uint16_t ido, const uint16_t l1, const complex_t *cc,
+ complex_t *ch, const complex_t *wa1, const complex_t *wa2,
+ const complex_t *wa3)
+{
+ uint16_t i, k, ac, ah;
+
+ if (ido == 1)
+ {
+ for (k = 0; k < l1; k++)
+ {
+ complex_t t1, t2, t3, t4;
+
+ ac = 4*k;
+ ah = k;
+
+ RE(t2) = RE(cc[ac]) + RE(cc[ac+2]);
+ RE(t1) = RE(cc[ac]) - RE(cc[ac+2]);
+ IM(t2) = IM(cc[ac]) + IM(cc[ac+2]);
+ IM(t1) = IM(cc[ac]) - IM(cc[ac+2]);
+ RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]);
+ IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]);
+ IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]);
+ RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]);
+
+ RE(ch[ah]) = RE(t2) + RE(t3);
+ RE(ch[ah+2*l1]) = RE(t2) - RE(t3);
+
+ IM(ch[ah]) = IM(t2) + IM(t3);
+ IM(ch[ah+2*l1]) = IM(t2) - IM(t3);
+
+ RE(ch[ah+l1]) = RE(t1) + RE(t4);
+ RE(ch[ah+3*l1]) = RE(t1) - RE(t4);
+
+ IM(ch[ah+l1]) = IM(t1) + IM(t4);
+ IM(ch[ah+3*l1]) = IM(t1) - IM(t4);
+ }
+ } else {
+ for (k = 0; k < l1; k++)
+ {
+ ac = 4*k*ido;
+ ah = k*ido;
+
+ for (i = 0; i < ido; i++)
+ {
+ complex_t c2, c3, c4, t1, t2, t3, t4;
+
+ RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]);
+ RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]);
+ IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]);
+ IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]);
+ RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]);
+ IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]);
+ IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]);
+ RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]);
+
+ RE(c2) = RE(t1) + RE(t4);
+ RE(c4) = RE(t1) - RE(t4);
+
+ IM(c2) = IM(t1) + IM(t4);
+ IM(c4) = IM(t1) - IM(t4);
+
+ RE(ch[ah+i]) = RE(t2) + RE(t3);
+ RE(c3) = RE(t2) - RE(t3);
+
+ IM(ch[ah+i]) = IM(t2) + IM(t3);
+ IM(c3) = IM(t2) - IM(t3);
+
+#if 1
+ ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]),
+ IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i]));
+ ComplexMult(&IM(ch[ah+i+2*l1*ido]), &RE(ch[ah+i+2*l1*ido]),
+ IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i]));
+ ComplexMult(&IM(ch[ah+i+3*l1*ido]), &RE(ch[ah+i+3*l1*ido]),
+ IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i]));
+#else
+ ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]),
+ RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i]));
+ ComplexMult(&RE(ch[ah+i+2*l1*ido]), &IM(ch[ah+i+2*l1*ido]),
+ RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i]));
+ ComplexMult(&RE(ch[ah+i+3*l1*ido]), &IM(ch[ah+i+3*l1*ido]),
+ RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i]));
+#endif
+ }
+ }
+ }
+}
+
+static void passf4neg(const uint16_t ido, const uint16_t l1, const complex_t *cc,
+ complex_t *ch, const complex_t *wa1, const complex_t *wa2,
+ const complex_t *wa3)
+{
+ uint16_t i, k, ac, ah;
+
+ if (ido == 1)
+ {
+ for (k = 0; k < l1; k++)
+ {
+ complex_t t1, t2, t3, t4;
+
+ ac = 4*k;
+ ah = k;
+
+ RE(t2) = RE(cc[ac]) + RE(cc[ac+2]);
+ RE(t1) = RE(cc[ac]) - RE(cc[ac+2]);
+ IM(t2) = IM(cc[ac]) + IM(cc[ac+2]);
+ IM(t1) = IM(cc[ac]) - IM(cc[ac+2]);
+ RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]);
+ IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]);
+ IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]);
+ RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]);
+
+ RE(ch[ah]) = RE(t2) + RE(t3);
+ RE(ch[ah+2*l1]) = RE(t2) - RE(t3);
+
+ IM(ch[ah]) = IM(t2) + IM(t3);
+ IM(ch[ah+2*l1]) = IM(t2) - IM(t3);
+
+ RE(ch[ah+l1]) = RE(t1) - RE(t4);
+ RE(ch[ah+3*l1]) = RE(t1) + RE(t4);
+
+ IM(ch[ah+l1]) = IM(t1) - IM(t4);
+ IM(ch[ah+3*l1]) = IM(t1) + IM(t4);
+ }
+ } else {
+ for (k = 0; k < l1; k++)
+ {
+ ac = 4*k*ido;
+ ah = k*ido;
+
+ for (i = 0; i < ido; i++)
+ {
+ complex_t c2, c3, c4, t1, t2, t3, t4;
+
+ RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]);
+ RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]);
+ IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]);
+ IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]);
+ RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]);
+ IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]);
+ IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]);
+ RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]);
+
+ RE(c2) = RE(t1) - RE(t4);
+ RE(c4) = RE(t1) + RE(t4);
+
+ IM(c2) = IM(t1) - IM(t4);
+ IM(c4) = IM(t1) + IM(t4);
+
+ RE(ch[ah+i]) = RE(t2) + RE(t3);
+ RE(c3) = RE(t2) - RE(t3);
+
+ IM(ch[ah+i]) = IM(t2) + IM(t3);
+ IM(c3) = IM(t2) - IM(t3);
+
+#if 1
+ ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]),
+ RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i]));
+ ComplexMult(&RE(ch[ah+i+2*l1*ido]), &IM(ch[ah+i+2*l1*ido]),
+ RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i]));
+ ComplexMult(&RE(ch[ah+i+3*l1*ido]), &IM(ch[ah+i+3*l1*ido]),
+ RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i]));
+#else
+ ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]),
+ IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i]));
+ ComplexMult(&IM(ch[ah+i+2*l1*ido]), &RE(ch[ah+i+2*l1*ido]),
+ IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i]));
+ ComplexMult(&IM(ch[ah+i+3*l1*ido]), &RE(ch[ah+i+3*l1*ido]),
+ IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i]));
+#endif
+ }
+ }
+ }
+}
+
+static void passf5(const uint16_t ido, const uint16_t l1, const complex_t *cc,
+ complex_t *ch, const complex_t *wa1, const complex_t *wa2, const complex_t *wa3,
+ const complex_t *wa4, const int8_t isign)
+{
+ static real_t tr11 = FRAC_CONST(0.309016994374947);
+ static real_t ti11 = FRAC_CONST(0.951056516295154);
+ static real_t tr12 = FRAC_CONST(-0.809016994374947);
+ static real_t ti12 = FRAC_CONST(0.587785252292473);
+ uint16_t i, k, ac, ah;
+ complex_t c2, c3, c4, c5, d3, d4, d5, d2, t2, t3, t4, t5;
+
+ if (ido == 1)
+ {
+ if (isign == 1)
+ {
+ for (k = 0; k < l1; k++)
+ {
+ ac = 5*k + 1;
+ ah = k;
+
+ RE(t2) = RE(cc[ac]) + RE(cc[ac+3]);
+ IM(t2) = IM(cc[ac]) + IM(cc[ac+3]);
+ RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]);
+ IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]);
+ RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]);
+ IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]);
+ RE(t5) = RE(cc[ac]) - RE(cc[ac+3]);
+ IM(t5) = IM(cc[ac]) - IM(cc[ac+3]);
+
+ RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
+ IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
+
+ RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
+ IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
+ RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
+ IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
+
+ ComplexMult(&RE(c5), &RE(c4),
+ ti11, ti12, RE(t5), RE(t4));
+ ComplexMult(&IM(c5), &IM(c4),
+ ti11, ti12, IM(t5), IM(t4));
+
+ RE(ch[ah+l1]) = RE(c2) - IM(c5);
+ IM(ch[ah+l1]) = IM(c2) + RE(c5);
+ RE(ch[ah+2*l1]) = RE(c3) - IM(c4);
+ IM(ch[ah+2*l1]) = IM(c3) + RE(c4);
+ RE(ch[ah+3*l1]) = RE(c3) + IM(c4);
+ IM(ch[ah+3*l1]) = IM(c3) - RE(c4);
+ RE(ch[ah+4*l1]) = RE(c2) + IM(c5);
+ IM(ch[ah+4*l1]) = IM(c2) - RE(c5);
+ }
+ } else {
+ for (k = 0; k < l1; k++)
+ {
+ ac = 5*k + 1;
+ ah = k;
+
+ RE(t2) = RE(cc[ac]) + RE(cc[ac+3]);
+ IM(t2) = IM(cc[ac]) + IM(cc[ac+3]);
+ RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]);
+ IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]);
+ RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]);
+ IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]);
+ RE(t5) = RE(cc[ac]) - RE(cc[ac+3]);
+ IM(t5) = IM(cc[ac]) - IM(cc[ac+3]);
+
+ RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
+ IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
+
+ RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
+ IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
+ RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
+ IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
+
+ ComplexMult(&RE(c4), &RE(c5),
+ ti12, ti11, RE(t5), RE(t4));
+ ComplexMult(&IM(c4), &IM(c5),
+ ti12, ti11, IM(t5), IM(t4));
+
+ RE(ch[ah+l1]) = RE(c2) + IM(c5);
+ IM(ch[ah+l1]) = IM(c2) - RE(c5);
+ RE(ch[ah+2*l1]) = RE(c3) + IM(c4);
+ IM(ch[ah+2*l1]) = IM(c3) - RE(c4);
+ RE(ch[ah+3*l1]) = RE(c3) - IM(c4);
+ IM(ch[ah+3*l1]) = IM(c3) + RE(c4);
+ RE(ch[ah+4*l1]) = RE(c2) - IM(c5);
+ IM(ch[ah+4*l1]) = IM(c2) + RE(c5);
+ }
+ }
+ } else {
+ if (isign == 1)
+ {
+ for (k = 0; k < l1; k++)
+ {
+ for (i = 0; i < ido; i++)
+ {
+ ac = i + (k*5 + 1) * ido;
+ ah = i + k * ido;
+
+ RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]);
+ IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]);
+ RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]);
+ IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]);
+ RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]);
+ IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]);
+ RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]);
+ IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]);
+
+ RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3);
+ IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3);
+
+ RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
+ IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
+ RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
+ IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
+
+ ComplexMult(&RE(c5), &RE(c4),
+ ti11, ti12, RE(t5), RE(t4));
+ ComplexMult(&IM(c5), &IM(c4),
+ ti11, ti12, IM(t5), IM(t4));
+
+ IM(d2) = IM(c2) + RE(c5);
+ IM(d3) = IM(c3) + RE(c4);
+ RE(d4) = RE(c3) + IM(c4);
+ RE(d5) = RE(c2) + IM(c5);
+ RE(d2) = RE(c2) - IM(c5);
+ IM(d5) = IM(c2) - RE(c5);
+ RE(d3) = RE(c3) - IM(c4);
+ IM(d4) = IM(c3) - RE(c4);
+
+#if 1
+ ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
+ IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
+ ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
+ IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
+ ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]),
+ IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i]));
+ ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]),
+ IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i]));
+#else
+ ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
+ RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
+ ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
+ RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
+ ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]),
+ RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i]));
+ ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]),
+ RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i]));
+#endif
+ }
+ }
+ } else {
+ for (k = 0; k < l1; k++)
+ {
+ for (i = 0; i < ido; i++)
+ {
+ ac = i + (k*5 + 1) * ido;
+ ah = i + k * ido;
+
+ RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]);
+ IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]);
+ RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]);
+ IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]);
+ RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]);
+ IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]);
+ RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]);
+ IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]);
+
+ RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3);
+ IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3);
+
+ RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
+ IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
+ RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
+ IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
+
+ ComplexMult(&RE(c4), &RE(c5),
+ ti12, ti11, RE(t5), RE(t4));
+ ComplexMult(&IM(c4), &IM(c5),
+ ti12, ti11, IM(t5), IM(t4));
+
+ IM(d2) = IM(c2) - RE(c5);
+ IM(d3) = IM(c3) - RE(c4);
+ RE(d4) = RE(c3) - IM(c4);
+ RE(d5) = RE(c2) - IM(c5);
+ RE(d2) = RE(c2) + IM(c5);
+ IM(d5) = IM(c2) + RE(c5);
+ RE(d3) = RE(c3) + IM(c4);
+ IM(d4) = IM(c3) + RE(c4);
+
+#if 1
+ ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
+ RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
+ ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
+ RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
+ ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]),
+ RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i]));
+ ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]),
+ RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i]));
+#else
+ ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
+ IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
+ ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
+ IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
+ ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]),
+ IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i]));
+ ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]),
+ IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i]));
+#endif
+ }
+ }
+ }
+ }
+}
+
+
+/*----------------------------------------------------------------------
+ cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
+ ----------------------------------------------------------------------*/
+
+static INLINE void cfftf1pos(uint16_t n, complex_t *c, complex_t *ch,
+ const uint16_t *ifac, const complex_t *wa,
+ const int8_t isign)
+{
+ uint16_t i;
+ uint16_t k1, l1, l2;
+ uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
+
+ nf = ifac[1];
+ na = 0;
+ l1 = 1;
+ iw = 0;
+
+ for (k1 = 2; k1 <= nf+1; k1++)
+ {
+ ip = ifac[k1];
+ l2 = ip*l1;
+ ido = n / l2;
+ idl1 = ido*l1;
+
+ switch (ip)
+ {
+ case 4:
+ ix2 = iw + ido;
+ ix3 = ix2 + ido;
+
+ if (na == 0)
+ passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]);
+ else
+ passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]);
+
+ na = 1 - na;
+ break;
+ case 2:
+ if (na == 0)
+ passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]);
+ else
+ passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]);
+
+ na = 1 - na;
+ break;
+ case 3:
+ ix2 = iw + ido;
+
+ if (na == 0)
+ passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign);
+ else
+ passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign);
+
+ na = 1 - na;
+ break;
+ case 5:
+ ix2 = iw + ido;
+ ix3 = ix2 + ido;
+ ix4 = ix3 + ido;
+
+ if (na == 0)
+ passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
+ else
+ passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
+
+ na = 1 - na;
+ break;
+ }
+
+ l1 = l2;
+ iw += (ip-1) * ido;
+ }
+
+ if (na == 0)
+ return;
+
+ for (i = 0; i < n; i++)
+ {
+ RE(c[i]) = RE(ch[i]);
+ IM(c[i]) = IM(ch[i]);
+ }
+}
+
+static INLINE void cfftf1neg(uint16_t n, complex_t *c, complex_t *ch,
+ const uint16_t *ifac, const complex_t *wa,
+ const int8_t isign)
+{
+ uint16_t i;
+ uint16_t k1, l1, l2;
+ uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
+
+ nf = ifac[1];
+ na = 0;
+ l1 = 1;
+ iw = 0;
+
+ for (k1 = 2; k1 <= nf+1; k1++)
+ {
+ ip = ifac[k1];
+ l2 = ip*l1;
+ ido = n / l2;
+ idl1 = ido*l1;
+
+ switch (ip)
+ {
+ case 4:
+ ix2 = iw + ido;
+ ix3 = ix2 + ido;
+
+ if (na == 0)
+ passf4neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]);
+ else
+ passf4neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]);
+
+ na = 1 - na;
+ break;
+ case 2:
+ if (na == 0)
+ passf2neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]);
+ else
+ passf2neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]);
+
+ na = 1 - na;
+ break;
+ case 3:
+ ix2 = iw + ido;
+
+ if (na == 0)
+ passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign);
+ else
+ passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign);
+
+ na = 1 - na;
+ break;
+ case 5:
+ ix2 = iw + ido;
+ ix3 = ix2 + ido;
+ ix4 = ix3 + ido;
+
+ if (na == 0)
+ passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
+ else
+ passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
+
+ na = 1 - na;
+ break;
+ }
+
+ l1 = l2;
+ iw += (ip-1) * ido;
+ }
+
+ if (na == 0)
+ return;
+
+ for (i = 0; i < n; i++)
+ {
+ RE(c[i]) = RE(ch[i]);
+ IM(c[i]) = IM(ch[i]);
+ }
+}
+
+void cfftf(cfft_info *cfft, complex_t *c)
+{
+ cfftf1neg(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, -1);
+}
+
+void cfftb(cfft_info *cfft, complex_t *c)
+{
+ cfftf1pos(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, +1);
+}
+
+static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac)
+{
+ static uint16_t ntryh[4] = {3, 4, 2, 5};
+#ifndef FIXED_POINT
+ real_t arg, argh, argld, fi;
+ uint16_t ido, ipm;
+ uint16_t i1, k1, l1, l2;
+ uint16_t ld, ii, ip;
+#endif
+ uint16_t ntry = 0, i, j;
+ uint16_t ib;
+ uint16_t nf, nl, nq, nr;
+
+ nl = n;
+ nf = 0;
+ j = 0;
+
+startloop:
+ j++;
+
+ if (j <= 4)
+ ntry = ntryh[j-1];
+ else
+ ntry += 2;
+
+ do
+ {
+ nq = nl / ntry;
+ nr = nl - ntry*nq;
+
+ if (nr != 0)
+ goto startloop;
+
+ nf++;
+ ifac[nf+1] = ntry;
+ nl = nq;
+
+ if (ntry == 2 && nf != 1)
+ {
+ for (i = 2; i <= nf; i++)
+ {
+ ib = nf - i + 2;
+ ifac[ib+1] = ifac[ib];
+ }
+ ifac[2] = 2;
+ }
+ } while (nl != 1);
+
+ ifac[0] = n;
+ ifac[1] = nf;
+
+#ifndef FIXED_POINT
+ argh = (real_t)2.0*(real_t)M_PI / (real_t)n;
+ i = 0;
+ l1 = 1;
+
+ for (k1 = 1; k1 <= nf; k1++)
+ {
+ ip = ifac[k1+1];
+ ld = 0;
+ l2 = l1*ip;
+ ido = n / l2;
+ ipm = ip - 1;
+
+ for (j = 0; j < ipm; j++)
+ {
+ i1 = i;
+ RE(wa[i]) = 1.0;
+ IM(wa[i]) = 0.0;
+ ld += l1;
+ fi = 0;
+ argld = ld*argh;
+
+ for (ii = 0; ii < ido; ii++)
+ {
+ i++;
+ fi++;
+ arg = fi * argld;
+ RE(wa[i]) = (real_t)cos(arg);
+#if 1
+ IM(wa[i]) = (real_t)sin(arg);
+#else
+ IM(wa[i]) = (real_t)-sin(arg);
+#endif
+ }
+
+ if (ip > 5)
+ {
+ RE(wa[i1]) = RE(wa[i]);
+ IM(wa[i1]) = IM(wa[i]);
+ }
+ }
+ l1 = l2;
+ }
+#endif
+}
+
+cfft_info *cffti(uint16_t n)
+{
+ cfft_info *cfft = (cfft_info*)faad_malloc(sizeof(cfft_info));
+
+ cfft->n = n;
+ cfft->work = (complex_t*)faad_malloc(n*sizeof(complex_t));
+
+#ifndef FIXED_POINT
+ cfft->tab = (complex_t*)faad_malloc(n*sizeof(complex_t));
+
+ cffti1(n, cfft->tab, cfft->ifac);
+#else
+ cffti1(n, NULL, cfft->ifac);
+
+ switch (n)
+ {
+ case 64: cfft->tab = (complex_t*)cfft_tab_64; break;
+ case 512: cfft->tab = (complex_t*)cfft_tab_512; break;
+#ifdef LD_DEC
+ case 256: cfft->tab = (complex_t*)cfft_tab_256; break;
+#endif
+
+#ifdef ALLOW_SMALL_FRAMELENGTH
+ case 60: cfft->tab = (complex_t*)cfft_tab_60; break;
+ case 480: cfft->tab = (complex_t*)cfft_tab_480; break;
+#ifdef LD_DEC
+ case 240: cfft->tab = (complex_t*)cfft_tab_240; break;
+#endif
+#endif
+ case 128: cfft->tab = (complex_t*)cfft_tab_128; break;
+ }
+#endif
+
+ return cfft;
+}
+
+void cfftu(cfft_info *cfft)
+{
+ if (cfft->work) faad_free(cfft->work);
+#ifndef FIXED_POINT
+ if (cfft->tab) faad_free(cfft->tab);
+#endif
+
+ if (cfft) faad_free(cfft);
+}
+
|