ref: 1a8846b50eb297b9faa06e1273b1589a7ee833ca
parent: 1f194f58d0de59460c8fdeda7eac01a9573f2836
author: menno <menno>
date: Sun Jul 14 15:11:12 EDT 2002
Changed fft; use of the FFTW library is not needed anymore, but stays as an option
--- a/libfaad/Makefile.am
+++ b/libfaad/Makefile.am
@@ -3,13 +3,13 @@
VPATH = @srcdir@:codebook:
-libfaad_la_SOURCES = bits.c data.c decoder.c drc.c error.c filtbank.c \
+libfaad_la_SOURCES = bits.c cfft.c data.c decoder.c drc.c error.c filtbank.c \
ic_predict.c is.c lt_predict.c mdct.c mp4.c ms.c output.c pns.c \
pulse.c specrec.c syntax.c tns.c reordered_spectral_data.c \
hcb_1.c hcb_2.c hcb_3.c hcb_4.c hcb_5.c hcb_6.c hcb_7.c hcb_8.c \
hcb_9.c hcb_10.c hcb_11.c hcb_sf.c
-libfaad_la_INCLUDES = analysis.h bits.h data.h decoder.h drc.h error.h \
+libfaad_la_INCLUDES = analysis.h bits.h cfft.h data.h decoder.h drc.h error.h \
filtbank.h huffman.h ic_predict.h is.h kbd_win.h lt_predict.h mdct.h mp4.h \
ms.h output.h pns.h pulse.h specrec.h syntax.h tns.h \
codebook/hcb.h
--- /dev/null
+++ b/libfaad/cfft.c
@@ -1,0 +1,671 @@
+/*
+** FAAD - Freeware Advanced Audio Decoder
+** Copyright (C) 2002 M. Bakker
+**
+** 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.
+**
+** $Id: cfft.c,v 1.1 2002/07/14 19:11:11 menno Exp $
+**/
+
+/*
+ * Algorithmically based on Fortran-77 FFTPACK
+ * by Paul N. Swarztrauber(Version 4, 1985).
+ */
+
+/* isign is +1 for backward and -1 for forward transforms */
+
+
+#include "common.h"
+#include <stdlib.h>
+
+#include "cfft.h"
+
+
+/*----------------------------------------------------------------------
+ passf2, passf3, passf4, passf5, passf. Complex FFT passes fwd and bwd.
+ ----------------------------------------------------------------------*/
+
+static void passf2(uint16_t ido, uint16_t l1, real_t *cc, real_t *ch,
+ real_t *wa1, int8_t isign)
+{+ uint16_t i, k, ah, ac;
+ real_t ti2, tr2;
+
+ if (ido <= 2)
+ {+ for (k = 0; k < l1; k++)
+ {+ ah = k*ido;
+ ac = 2*k*ido;
+ ch[ah] = cc[ac] + cc[ac+ido];
+ ch[ah+ido*l1] = cc[ac] - cc[ac+ido];
+ ch[ah+1] = cc[ac+1] + cc[ac+ido+1];
+ ch[ah+ido*l1+1] = cc[ac+1] - cc[ac+ido+1];
+ }
+ } else {+ for (k = 0; k < l1; k++)
+ {+ for(i = 0; i < ido-1; i += 2)
+ {+ ah = i + k*ido;
+ ac = i + 2*k*ido;
+ ch[ah] = cc[ac] + cc[ac+ido];
+ tr2 = cc[ac] - cc[ac+ido];
+ ch[ah+1] = cc[ac+1] + cc[ac+1+ido];
+ ti2 = cc[ac+1] - cc[ac+1+ido];
+ ch[ah+l1*ido+1] = wa1[i]*ti2 + isign*wa1[i+1]*tr2;
+ ch[ah+l1*ido] = wa1[i]*tr2 - isign*wa1[i+1]*ti2;
+ }
+ }
+ }
+}
+
+
+static void passf3(uint16_t ido, uint16_t l1, real_t *cc, real_t *ch,
+ real_t *wa1, real_t *wa2, int8_t isign)
+{+ static real_t taur = -0.5;
+ static real_t taui = 0.866025403784439;
+ uint16_t i, k, ac, ah;
+ real_t ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
+
+ if (ido == 2)
+ {+ for (k = 1; k <= l1; k++)
+ {+ ac = (3*k-2) * ido;
+ tr2 = cc[ac] + cc[ac+ido];
+ cr2 = cc[ac-ido] + taur*tr2;
+ ah = (k-1) * ido;
+ ch[ah] = cc[ac-ido] + tr2;
+
+ ti2 = cc[ac+1] + cc[ac+ido+1];
+ ci2 = cc[ac-ido+1] + taur*ti2;
+ ch[ah+1] = cc[ac-ido+1] + ti2;
+
+ cr3 = isign * taui * (cc[ac] - cc[ac+ido]);
+ ci3 = isign * taui * (cc[ac+1] - cc[ac+ido+1]);
+ ch[ah+l1*ido] = cr2 - ci3;
+ ch[ah+2*l1*ido] = cr2 + ci3;
+ ch[ah+l1*ido+1] = ci2 + cr3;
+ ch[ah+2*l1*ido+1] = ci2 - cr3;
+ }
+ } else {+ for (k = 1; k <= l1; k++)
+ {+ for (i = 0; i < ido-1; i += 2)
+ {+ ac = i + (3*k-2) * ido;
+ tr2 = cc[ac] + cc[ac+ido];
+ cr2 = cc[ac-ido] + taur*tr2;
+ ah = i + (k-1) * ido;
+ ch[ah] = cc[ac-ido] + tr2;
+ ti2 = cc[ac+1] + cc[ac+ido+1];
+ ci2 = cc[ac-ido+1] + taur*ti2;
+ ch[ah+1] = cc[ac-ido+1] + ti2;
+ cr3 = isign * taui * (cc[ac] - cc[ac+ido]);
+ ci3 = isign * taui * (cc[ac+1] - cc[ac+ido+1]);
+ dr2 = cr2 - ci3;
+ dr3 = cr2 + ci3;
+ di2 = ci2 + cr3;
+ di3 = ci2 - cr3;
+ ch[ah+l1*ido+1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
+ ch[ah+l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
+ ch[ah+2*l1*ido+1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
+ ch[ah+2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
+ }
+ }
+ }
+}
+
+
+static void passf4(uint16_t ido, uint16_t l1, real_t *cc, real_t *ch,
+ real_t *wa1, real_t *wa2, real_t *wa3, int8_t isign)
+{+ uint16_t i, k, ac, ah;
+ real_t ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2,
+ tr3, tr4;
+
+ if (ido == 2)
+ {+ for (k = 0; k < l1; k++)
+ {+ ac = 4*k*ido + 1;
+ ti1 = cc[ac] - cc[ac+2*ido];
+ ti2 = cc[ac] + cc[ac+2*ido];
+ tr4 = cc[ac+3*ido] - cc[ac+ido];
+ ti3 = cc[ac+ido] + cc[ac+3*ido];
+ tr1 = cc[ac-1] - cc[ac+2*ido-1];
+ tr2 = cc[ac-1] + cc[ac+2*ido-1];
+ ti4 = cc[ac+ido-1] - cc[ac+3*ido-1];
+ tr3 = cc[ac+ido-1] + cc[ac+3*ido-1];
+ ah = k*ido;
+ ch[ah] = tr2 + tr3;
+ ch[ah+2*l1*ido] = tr2 - tr3;
+ ch[ah+1] = ti2 + ti3;
+ ch[ah+2*l1*ido+1] = ti2 - ti3;
+ ch[ah+l1*ido] = tr1 + isign*tr4;
+ ch[ah+3*l1*ido] = tr1 - isign*tr4;
+ ch[ah+l1*ido+1] = ti1 + isign*ti4;
+ ch[ah+3*l1*ido+1] = ti1 - isign*ti4;
+ }
+ } else {+ for (k = 0; k < l1; k++)
+ {+ for (i = 0; i < ido-1; i += 2)
+ {+ ac = i + 1 + 4*k*ido;
+ ti1 = cc[ac] - cc[ac+2*ido];
+ ti2 = cc[ac] + cc[ac+2*ido];
+ ti3 = cc[ac+ido] + cc[ac+3*ido];
+ tr4 = cc[ac+3*ido] - cc[ac+ido];
+ tr1 = cc[ac-1] - cc[ac+2*ido-1];
+ tr2 = cc[ac-1] + cc[ac+2*ido-1];
+ ti4 = cc[ac+ido-1] - cc[ac+3*ido-1];
+ tr3 = cc[ac+ido-1] + cc[ac+3*ido-1];
+ ah = i + k*ido;
+ ch[ah] = tr2 + tr3;
+ cr3 = tr2 - tr3;
+ ch[ah+1] = ti2 + ti3;
+ ci3 = ti2 - ti3;
+ cr2 = tr1 + isign*tr4;
+ cr4 = tr1 - isign*tr4;
+ ci2 = ti1 + isign*ti4;
+ ci4 = ti1 - isign*ti4;
+ ch[ah+l1*ido] = wa1[i]*cr2 - isign*wa1[i+1]*ci2;
+ ch[ah+l1*ido+1] = wa1[i]*ci2 + isign*wa1[i+1]*cr2;
+ ch[ah+2*l1*ido] = wa2[i]*cr3 - isign*wa2[i+1]*ci3;
+ ch[ah+2*l1*ido+1] = wa2[i]*ci3 + isign*wa2[i+1]*cr3;
+ ch[ah+3*l1*ido] = wa3[i]*cr4 - isign*wa3[i+1]*ci4;
+ ch[ah+3*l1*ido+1] = wa3[i]*ci4 + isign*wa3[i+1]*cr4;
+ }
+ }
+ }
+}
+
+
+static void passf5(uint16_t ido, uint16_t l1, real_t *cc, real_t *ch,
+ real_t *wa1, real_t *wa2, real_t *wa3, real_t *wa4,
+ int8_t isign)
+{+ static real_t tr11 = 0.309016994374947;
+ static real_t ti11 = 0.951056516295154;
+ static real_t tr12 = -0.809016994374947;
+ static real_t ti12 = 0.587785252292473;
+ uint16_t i, k, ac, ah;
+ real_t ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4,
+ ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
+
+ if (ido == 2)
+ {+ for (k = 1; k <= l1; ++k)
+ {+ ac = (5*k-4) * ido + 1;
+ ti5 = cc[ac] - cc[ac+3*ido];
+ ti2 = cc[ac] + cc[ac+3*ido];
+ ti4 = cc[ac+ido] - cc[ac+2*ido];
+ ti3 = cc[ac+ido] + cc[ac+2*ido];
+ tr5 = cc[ac-1] - cc[ac+3*ido-1];
+ tr2 = cc[ac-1] + cc[ac+3*ido-1];
+ tr4 = cc[ac+ido-1] - cc[ac+2*ido-1];
+ tr3 = cc[ac+ido-1] + cc[ac+2*ido-1];
+ ah = (k-1) * ido;
+ ch[ah] = cc[ac-ido-1] + tr2 + tr3;
+ ch[ah+1] = cc[ac-ido] + ti2 + ti3;
+ cr2 = cc[ac-ido-1] + tr11*tr2 + tr12*tr3;
+ ci2 = cc[ac-ido] + tr11*ti2 + tr12*ti3;
+ cr3 = cc[ac-ido-1] + tr12*tr2 + tr11*tr3;
+ ci3 = cc[ac-ido] + tr12*ti2 + tr11*ti3;
+ cr5 = isign * (ti11*tr5 + ti12*tr4);
+ ci5 = isign * (ti11*ti5 + ti12*ti4);
+ cr4 = isign * (ti12*tr5 - ti11*tr4);
+ ci4 = isign * (ti12*ti5 - ti11*ti4);
+ ch[ah+l1*ido] = cr2 - ci5;
+ ch[ah+4*l1*ido] = cr2 + ci5;
+ ch[ah+l1*ido+1] = ci2 + cr5;
+ ch[ah+2*l1*ido+1]=ci3 + cr4;
+ ch[ah+2*l1*ido] = cr3 - ci4;
+ ch[ah+3*l1*ido] = cr3 + ci4;
+ ch[ah+3*l1*ido+1] = ci3 - cr4;
+ ch[ah+4*l1*ido+1] = ci2 - cr5;
+ }
+ } else {+ for (k = 1; k <= l1; k++)
+ {+ for (i = 0; i < ido-1; i += 2)
+ {+ ac = i + 1 + (k*5-4) * ido;
+ ti5 = cc[ac] - cc[ac+3*ido];
+ ti2 = cc[ac] + cc[ac+3*ido];
+ ti4 = cc[ac+ido] - cc[ac+2*ido];
+ ti3 = cc[ac+ido] + cc[ac+2*ido];
+ tr5 = cc[ac-1] - cc[ac+3*ido-1];
+ tr2 = cc[ac-1] + cc[ac+3*ido-1];
+ tr4 = cc[ac+ido-1] - cc[ac+2*ido-1];
+ tr3 = cc[ac+ido-1] + cc[ac+2*ido-1];
+ ah = i + (k-1) * ido;
+ ch[ah] = cc[ac-ido-1] + tr2 + tr3;
+ ch[ah+1] = cc[ac-ido] + ti2 + ti3;
+ cr2 = cc[ac-ido-1] + tr11*tr2 + tr12*tr3;
+ ci2 = cc[ac-ido] + tr11*ti2 + tr12*ti3;
+ cr3 = cc[ac-ido-1] + tr12*tr2 + tr11*tr3;
+ ci3 = cc[ac-ido] + tr12*ti2 + tr11*ti3;
+ cr5 = isign * (ti11*tr5 + ti12*tr4);
+ ci5 = isign * (ti11*ti5 + ti12*ti4);
+ cr4 = isign * (ti12*tr5 - ti11*tr4);
+ ci4 = isign * (ti12*ti5 - ti11*ti4);
+ dr3 = cr3 - ci4;
+ dr4 = cr3 + ci4;
+ di3 = ci3 + cr4;
+ di4 = ci3 - cr4;
+ dr5 = cr2 + ci5;
+ dr2 = cr2 - ci5;
+ di5 = ci2 - cr5;
+ di2 = ci2 + cr5;
+ ch[ah+l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
+ ch[ah+l1*ido+1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
+ ch[ah+2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
+ ch[ah+2*l1*ido+1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
+ ch[ah+3*l1*ido] = wa3[i]*dr4 - isign*wa3[i+1]*di4;
+ ch[ah+3*l1*ido+1] = wa3[i]*di4 + isign*wa3[i+1]*dr4;
+ ch[ah+4*l1*ido] = wa4[i]*dr5 - isign*wa4[i+1]*di5;
+ ch[ah+4*l1*ido+1] = wa4[i]*di5 + isign*wa4[i+1]*dr5;
+ }
+ }
+ }
+}
+
+
+static void passf(uint16_t *nac, uint16_t ido, uint16_t ip, uint16_t l1,
+ uint16_t idl1, real_t *cc, real_t *ch, real_t *wa,
+ int8_t isign)
+{+ uint16_t idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, nt, idj, idl;
+ uint16_t inc, idp;
+ real_t wai, war;
+
+ idot = ido / 2;
+ nt = ip*idl1;
+ ipph = (ip+1) / 2;
+ idp = ip*ido;
+
+ if (ido >= l1)
+ {+ for (j = 1; j < ipph; j++)
+ {+ jc = ip - j;
+
+ for (k = 0; k < l1; k++)
+ {+ for (i = 0; i < ido; i++)
+ {+ ch[i+(k+j*l1)*ido] = cc[i+(j+k*ip)*ido] + cc[i+(jc+k*ip)*ido];
+ ch[i+(k+jc*l1)*ido] = cc[i+(j+k*ip)*ido] - cc[i+(jc+k*ip)*ido];
+ }
+ }
+ }
+
+ for (k = 0; k < l1; k++)
+ {+ for (i = 0; i < ido; i++)
+ ch[i+k*ido] = cc[i+k*ip*ido];
+ }
+ } else {+ for (j = 1; j < ipph; j++)
+ {+ jc = ip - j;
+
+ for (i = 0; i < ido; i++)
+ {+ for (k = 0; k < l1; k++)
+ {+ ch[i+(k+j*l1)*ido] = cc[i+(j+k*ip)*ido] + cc[i+(jc+k*ip)*ido];
+ ch[i+(k+jc*l1)*ido] = cc[i+(j+k*ip)*ido] - cc[i+(jc+k*ip)*ido];
+ }
+ }
+ }
+
+ for (i = 0; i < ido; i++)
+ {+ for (k = 0; k < l1; k++)
+ ch[i+k*ido] = cc[i+k*ip*ido];
+ }
+ }
+
+ idl = 2 - ido;
+ inc = 0;
+
+ for (l = 1; l < ipph; l++)
+ {+ lc = ip - l;
+ idl += ido;
+
+ for (ik = 0; ik < idl1; ik++)
+ {+ cc[ik+l*idl1] = ch[ik] + wa[idl-2]*ch[ik+idl1];
+ cc[ik+lc*idl1] = isign*wa[idl-1]*ch[ik+(ip-1)*idl1];
+ }
+
+ idlj = idl;
+ inc += ido;
+
+ for (j = 2; j < ipph; j++)
+ {+ jc = ip - j;
+ idlj += inc;
+
+ if (idlj > idp)
+ idlj -= idp;
+
+ war = wa[idlj-2];
+ wai = wa[idlj-1];
+
+ for (ik = 0; ik < idl1; ik++)
+ {+ cc[ik+l*idl1] += war*ch[ik+j*idl1];
+ cc[ik+lc*idl1] += isign*wai*ch[ik+jc*idl1];
+ }
+ }
+ }
+
+ for (j = 1; j < ipph; j++)
+ {+ for (ik = 0; ik < idl1; ik++)
+ ch[ik] += ch[ik+j*idl1];
+ }
+
+ for (j = 1; j < ipph; j++)
+ {+ jc = ip - j;
+
+ for (ik = 1; ik < idl1; ik += 2)
+ {+ ch[ik-1+j*idl1] = cc[ik-1+j*idl1] - cc[ik+jc*idl1];
+ ch[ik-1+jc*idl1] = cc[ik-1+j*idl1] + cc[ik+jc*idl1];
+ ch[ik+j*idl1] = cc[ik+j*idl1] + cc[ik-1+jc*idl1];
+ ch[ik+jc*idl1] = cc[ik+j*idl1] - cc[ik-1+jc*idl1];
+ }
+ }
+
+ *nac = 1;
+
+ if (ido == 2)
+ return;
+
+ *nac = 0;
+
+ for (ik = 0; ik < idl1; ik++)
+ cc[ik] = ch[ik];
+
+ for (j = 1; j < ip; j++)
+ {+ for (k = 0; k < l1; k++)
+ {+ cc[(k+j*l1)*ido+0] = ch[(k+j*l1)*ido+0];
+ cc[(k+j*l1)*ido+1] = ch[(k+j*l1)*ido+1];
+ }
+ }
+
+ if (idot <= l1)
+ {+ idij = 0;
+
+ for (j = 1; j < ip; j++)
+ {+ idij += 2;
+
+ for (i = 3; i < ido; i += 2)
+ {+ idij += 2;
+
+ for (k = 0; k < l1; k++)
+ {+ cc[i-1+(k+j*l1)*ido] = wa[idij-2] * ch[i-1+(k+j*l1)*ido] -
+ isign * wa[idij-1] * ch[i+(k+j*l1)*ido];
+ cc[i+(k+j*l1)*ido] = wa[idij-2] * ch[i+(k+j*l1)*ido] +
+ isign * wa[idij-1] * ch[i-1+(k+j*l1)*ido];
+ }
+ }
+ }
+ } else {+ idj = 2 - ido;
+
+ for (j = 1; j < ip; j++)
+ {+ idj += ido;
+
+ for (k = 0; k < l1; k++)
+ {+ idij = idj;
+
+ for (i = 3; i < ido; i += 2)
+ {+ idij += 2;
+ cc[i-1+(k+j*l1)*ido] = wa[idij-2] * ch[i-1+(k+j*l1)*ido] -
+ isign * wa[idij-1] * ch[i+(k+j*l1)*ido];
+ cc[i+(k+j*l1)*ido] = wa[idij-2] * ch[i+(k+j*l1)*ido] +
+ isign * wa[idij-1] * ch[i-1+(k+j*l1)*ido];
+ }
+ }
+ }
+ }
+}
+
+
+
+/*----------------------------------------------------------------------
+ cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
+ ----------------------------------------------------------------------*/
+
+INLINE void cfftf1(uint16_t n, real_t *c, real_t *ch, real_t *wa,
+ uint16_t *ifac, int8_t isign)
+{+ uint16_t idot, i;
+ uint16_t k1, l1, l2;
+ uint16_t na, nf, ip, iw, ix2, ix3, ix4, nac, 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;
+ idot = ido+ido;
+ idl1 = idot*l1;
+
+ switch (ip)
+ {+ case 2:
+ if (na == 0)
+ passf2(idot, l1, c, ch, &wa[iw], isign);
+ else
+ passf2(idot, l1, ch, c, &wa[iw], isign);
+
+ na = 1 - na;
+ break;
+ case 3:
+ ix2 = iw + idot;
+
+ if (na == 0)
+ passf3(idot, l1, c, ch, &wa[iw], &wa[ix2], isign);
+ else
+ passf3(idot, l1, ch, c, &wa[iw], &wa[ix2], isign);
+
+ na = 1 - na;
+ break;
+ case 4:
+ ix2 = iw + idot;
+ ix3 = ix2 + idot;
+
+ if (na == 0)
+ passf4(idot, l1, c, ch, &wa[iw], &wa[ix2], &wa[ix3], isign);
+ else
+ passf4(idot, l1, ch, c, &wa[iw], &wa[ix2], &wa[ix3], isign);
+
+ na = 1 - na;
+ break;
+ case 5:
+ ix2 = iw + idot;
+ ix3 = ix2 + idot;
+ ix4 = ix3 + idot;
+
+ if (na == 0)
+ passf5(idot, l1, c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
+ else
+ passf5(idot, l1, ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
+
+ na = 1 - na;
+ break;
+ default:
+ if (na == 0)
+ passf(&nac, idot, ip, l1, idl1, c, ch, &wa[iw], isign);
+ else
+ passf(&nac, idot, ip, l1, idl1, ch, c, &wa[iw], isign);
+
+ if (nac != 0)
+ na = 1 - na;
+ break;
+ }
+
+ l1 = l2;
+ iw += (ip-1) * idot;
+ }
+
+ if (na == 0)
+ return;
+
+ for (i = 0; i < 2*n; i++)
+ c[i] = ch[i];
+}
+
+void cfftf(cfft_info cfft, real_t *c)
+{+ cfftf1(cfft.n, c, cfft.work, cfft.tab, cfft.ifac, -1);
+}
+
+void cfftb(cfft_info cfft, real_t *c)
+{+ cfftf1(cfft.n, c, cfft.work, cfft.tab, cfft.ifac, +1);
+}
+
+static void cffti1(uint16_t n, real_t *wa, uint16_t *ifac)
+{+ static uint16_t ntryh[4] = {3, 4, 2, 5};+ real_t arg, argh, argld, fi;
+ uint16_t idot, ntry, i, j;
+ uint16_t i1, k1, l1, l2, ib;
+ uint16_t ld, ii, nf, ip, nl, nq, nr;
+ uint16_t ido, ipm;
+
+ 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;
+ argh = 2*M_PI / (real_t)n;
+ i = 1;
+ l1 = 1;
+
+ for (k1 = 1; k1 <= nf; k1++)
+ {+ ip = ifac[k1+1];
+ ld = 0;
+ l2 = l1*ip;
+ ido = n / l2;
+ idot = ido + ido + 2;
+ ipm = ip - 1;
+
+ for (j = 1; j <= ipm; j++)
+ {+ i1 = i;
+ wa[i-1] = 1;
+ wa[i] = 0;
+ ld += l1;
+ fi = 0;
+ argld = ld*argh;
+
+ for (ii = 4; ii <= idot; ii += 2)
+ {+ i += 2;
+ fi += 1;
+ arg = fi*argld;
+ wa[i-1] = cos(arg);
+ wa[i] = sin(arg);
+ }
+
+ if (ip > 5)
+ {+ wa[i1-1] = wa[i-1];
+ wa[i1] = wa[i];
+ }
+ }
+ l1 = l2;
+ }
+}
+
+cfft_info cffti(uint16_t n)
+{+ cfft_info cfft;
+
+ cfft.n = n;
+ cfft.work = malloc(2*n*sizeof(real_t));
+ cfft.tab = malloc(2*n*sizeof(real_t));
+
+ cffti1(n, cfft.tab, cfft.ifac);
+
+ return cfft;
+}
+
+void cfftu(cfft_info cfft)
+{+ if (cfft.work) free(cfft.work);
+ if (cfft.tab) free(cfft.tab);
+}
\ No newline at end of file
--- /dev/null
+++ b/libfaad/cfft.h
@@ -1,0 +1,63 @@
+/*
+** FAAD - Freeware Advanced Audio Decoder
+** Copyright (C) 2002 M. Bakker
+**
+** 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.
+**
+** $Id: cfft.h,v 1.1 2002/07/14 19:11:11 menno Exp $
+**/
+
+#ifndef __CFFT_H__
+#define __CFFT_H__
+
+#ifdef __cplusplus
+extern "C" {+#endif
+
+typedef struct
+{+ real_t *work;
+ real_t *tab;
+ uint16_t ifac[15];
+ uint16_t n;
+} cfft_info;
+
+void cfftf(cfft_info cfft, real_t *c); // complex transform
+void cfftb(cfft_info cfft, real_t *c); // its inverse
+cfft_info cffti(uint16_t n); // initializer of the above routine pair
+void cfftu(cfft_info cfft);
+
+
+static void passf2(uint16_t ido, uint16_t l1, real_t *cc, real_t *ch,
+ real_t *wa1, int8_t isign);
+static void passf3(uint16_t ido, uint16_t l1, real_t *cc, real_t *ch,
+ real_t *wa1, real_t *wa2, int8_t isign);
+static void passf4(uint16_t ido, uint16_t l1, real_t *cc, real_t *ch,
+ real_t *wa1, real_t *wa2, real_t *wa3, int8_t isign);
+static void passf5(uint16_t ido, uint16_t l1, real_t *cc, real_t *ch,
+ real_t *wa1, real_t *wa2, real_t *wa3, real_t *wa4,
+ int8_t isign);
+static void passf(uint16_t *nac, uint16_t ido, uint16_t ip, uint16_t l1,
+ uint16_t idl1, real_t *cc, real_t *ch, real_t *wa,
+ int8_t isign);
+INLINE void cfftf1(uint16_t n, real_t *c, real_t *ch, real_t *wa,
+ uint16_t *ifac, int8_t isign);
+static void cffti1(uint16_t n, real_t *wa, uint16_t *ifac);
+
+
+#ifdef __cplusplus
+}
+#endif
+#endif
--- a/libfaad/common.h
+++ b/libfaad/common.h
@@ -16,7 +16,7 @@
** along with this program; if not, write to the Free Software
** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
**
-** $Id: common.h,v 1.13 2002/05/30 18:31:51 menno Exp $
+** $Id: common.h,v 1.14 2002/07/14 19:11:11 menno Exp $
**/
#ifndef __COMMON_H__
@@ -53,6 +53,9 @@
#endif
/* COMPILE TIME DEFINITIONS */
+
+/* use the somewhat faster, but a lot larger FFTW library */
+//#define USE_FFTW
/* use double precision */
/* #define USE_DOUBLE_PRECISION */
--- a/libfaad/libfaad.dsp
+++ b/libfaad/libfaad.dsp
@@ -135,308 +135,16 @@
SOURCE=.\codebook\hcb_sf.c
# End Source File
# End Group
-# Begin Group "fftw"
-
-# PROP Default_Filter ""
# Begin Source File
-SOURCE=.\fftw\config.c
+SOURCE=.\bits.c
# End Source File
# Begin Source File
-SOURCE=.\fftw\executor.c
+SOURCE=.\cfft.c
# End Source File
# Begin Source File
-SOURCE=.\fftw\fftwf77.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fftwnd.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_1.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_10.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_11.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_12.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_13.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_14.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_15.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_16.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_2.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_3.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_32.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_4.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_5.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_6.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_64.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_7.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_8.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fn_9.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_1.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_10.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_11.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_12.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_13.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_14.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_15.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_16.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_2.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_3.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_32.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_4.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_5.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_6.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_64.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_7.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_8.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fni_9.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftw_10.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftw_16.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftw_2.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftw_3.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftw_32.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftw_4.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftw_5.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftw_6.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftw_64.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftw_7.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftw_8.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftw_9.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftwi_10.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftwi_16.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftwi_2.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftwi_3.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftwi_32.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftwi_4.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftwi_5.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftwi_6.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftwi_64.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftwi_7.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftwi_8.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\ftwi_9.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\generic.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\malloc.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\planner.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\putils.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\rader.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\timer.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\twiddle.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\wisdom.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\wisdomio.c
-# End Source File
-# End Group
-# Begin Source File
-
-SOURCE=.\bits.c
-# End Source File
-# Begin Source File
-
SOURCE=.\common.c
# End Source File
# Begin Source File
@@ -541,6 +249,10 @@
# End Source File
# Begin Source File
+SOURCE=.\cfft.h
+# End Source File
+# Begin Source File
+
SOURCE=.\common.h
# End Source File
# Begin Source File
@@ -566,14 +278,6 @@
# Begin Source File
SOURCE=.\fftw\f77_func.h
-# End Source File
-# Begin Source File
-
-SOURCE=".\fftw\fftw-int.h"
-# End Source File
-# Begin Source File
-
-SOURCE=.\fftw\fftw.h
# End Source File
# Begin Source File
--- a/libfaad/mdct.c
+++ b/libfaad/mdct.c
@@ -16,7 +16,7 @@
** along with this program; if not, write to the Free Software
** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
**
-** $Id: mdct.c,v 1.10 2002/05/24 17:26:12 menno Exp $
+** $Id: mdct.c,v 1.11 2002/07/14 19:11:11 menno Exp $
**/
/*
@@ -42,13 +42,16 @@
#include <stdlib.h>
#include <assert.h>
+#ifdef USE_FFTW
/* uses fftw (http://www.fftw.org) for very fast arbitrary-n FFT and IFFT */
#include <fftw.h>
+#else
+#include "cfft.h"
+#endif
#include "mdct.h"
-
void faad_mdct_init(mdct_info *mdct, uint16_t N)
{uint16_t k;
@@ -57,8 +60,13 @@
mdct->N = N;
mdct->sincos = (faad_sincos*)malloc(N/4*sizeof(faad_sincos));
+#ifdef USE_FFTW
mdct->Z1 = (fftw_complex*)malloc(N/4*sizeof(fftw_complex));
mdct->Z2 = (fftw_complex*)malloc(N/4*sizeof(fftw_complex));
+#else
+ mdct->Z1 = (real_t*)malloc(N/2*sizeof(real_t));
+ mdct->Z2 = (faad_complex*)malloc(N/4*sizeof(faad_complex));
+#endif
for (k = 0; k < N/4; k++)
{@@ -67,18 +75,27 @@
mdct->sincos[k].cos = -cos(angle);
}
+#ifdef USE_FFTW
mdct->plan_backward = fftw_create_plan(N/4, FFTW_BACKWARD, FFTW_ESTIMATE);
#ifdef LTP_DEC
mdct->plan_forward = fftw_create_plan(N/4, FFTW_FORWARD, FFTW_ESTIMATE);
#endif
+#else
+ /* own implementation */
+ mdct->cfft = cffti(N/4);
+#endif
}
void faad_mdct_end(mdct_info *mdct)
{+#ifdef USE_FFTW
fftw_destroy_plan(mdct->plan_backward);
#ifdef LTP_DEC
fftw_destroy_plan(mdct->plan_forward);
#endif
+#else
+ cfftu(mdct->cfft);
+#endif
if (mdct->Z2) free(mdct->Z2);
if (mdct->Z1) free(mdct->Z1);
@@ -89,9 +106,15 @@
{uint16_t k;
+#ifdef USE_FFTW
fftw_complex *Z1 = mdct->Z1;
fftw_complex *Z2 = mdct->Z2;
+#else
+ real_t *Z1 = mdct->Z1;
+ faad_complex *Z2 = mdct->Z2;
+#endif
faad_sincos *sincos = mdct->sincos;
+ real_t fftdata[1024];
uint16_t N = mdct->N;
uint16_t N2 = N >> 1;
@@ -106,18 +129,32 @@
uint16_t n = k << 1;
real_t x0 = X_in[ n];
real_t x1 = X_in[N2 - 1 - n];
+#ifdef USE_FFTW
Z1[k].re = MUL(fac, MUL(x1, sincos[k].cos) - MUL(x0, sincos[k].sin));
Z1[k].im = MUL(fac, MUL(x0, sincos[k].cos) + MUL(x1, sincos[k].sin));
+#else
+ Z1[2*k] = MUL(fac, MUL(x1, sincos[k].cos) - MUL(x0, sincos[k].sin));
+ Z1[2*k+1] = MUL(fac, MUL(x0, sincos[k].cos) + MUL(x1, sincos[k].sin));
+#endif
}
/* complex IFFT */
+#ifdef USE_FFTW
fftw_one(mdct->plan_backward, Z1, Z2);
+#else
+ cfftb(mdct->cfft, Z1);
+#endif
/* post-IFFT complex multiplication */
for (k = 0; k < N4; k++)
{+#ifdef USE_FFTW
real_t zr = Z2[k].re;
real_t zi = Z2[k].im;
+#else
+ real_t zr = Z1[2*k];
+ real_t zi = Z1[2*k+1];
+#endif
Z2[k].re = MUL(zr, sincos[k].cos) - MUL(zi, sincos[k].sin);
Z2[k].im = MUL(zi, sincos[k].cos) + MUL(zr, sincos[k].sin);
}
@@ -142,8 +179,12 @@
{uint16_t k;
+#ifdef USE_FFTW
fftw_complex *Z1 = mdct->Z1;
fftw_complex *Z2 = mdct->Z2;
+#else
+ real_t *Z1 = mdct->Z1;
+#endif
faad_sincos *sincos = mdct->sincos;
uint16_t N = mdct->N;
@@ -159,25 +200,44 @@
real_t zr = X_in[N - N4 - 1 - n] + X_in[N - N4 + n];
real_t zi = X_in[ N4 + n] - X_in[ N4 - 1 - n];
+#ifdef USE_FFTW
Z1[k ].re = -MUL(zr, sincos[k ].cos) - MUL(zi, sincos[k ].sin);
Z1[k ].im = -MUL(zi, sincos[k ].cos) + MUL(zr, sincos[k ].sin);
+#else
+ Z1[k*2 ] = -MUL(zr, sincos[k ].cos) - MUL(zi, sincos[k ].sin);
+ Z1[k*2+1 ] = -MUL(zi, sincos[k ].cos) + MUL(zr, sincos[k ].sin);
+#endif
zr = X_in[ N2 - 1 - n] - X_in[ n];
zi = X_in[ N2 + n] + X_in[N - 1 - n];
+#ifdef USE_FFTW
Z1[k + N8].re = -MUL(zr, sincos[k + N8].cos) - MUL(zi, sincos[k + N8].sin);
Z1[k + N8].im = -MUL(zi, sincos[k + N8].cos) + MUL(zr, sincos[k + N8].sin);
+#else
+ Z1[k*2 + N8] = -MUL(zr, sincos[k + N8].cos) - MUL(zi, sincos[k + N8].sin);
+ Z1[k*2+1 + N8] = -MUL(zi, sincos[k + N8].cos) + MUL(zr, sincos[k + N8].sin);
+#endif
}
/* complex FFT */
+#ifdef USE_FFTW
fftw_one(mdct->plan_forward, Z1, Z2);
+#else
+ cfftf(mdct->cfft, Z1);
+#endif
/* post-FFT complex multiplication */
for (k = 0; k < N4; k++)
{uint16_t n = k << 1;
+#ifdef USE_FFTW
real_t zr = MUL(2.0, MUL(Z2[k].re, sincos[k].cos) + MUL(Z2[k].im, sincos[k].sin));
real_t zi = MUL(2.0, MUL(Z2[k].im, sincos[k].cos) - MUL(Z2[k].re, sincos[k].sin));
+#else
+ real_t zr = MUL(2.0, MUL(Z1[k*2], sincos[k].cos) + MUL(Z1[k*2+1], sincos[k].sin));
+ real_t zi = MUL(2.0, MUL(Z1[k*2+1], sincos[k].cos) - MUL(Z1[k*2], sincos[k].sin));
+#endif
X_out[ n] = -zr;
X_out[N2 - 1 - n] = zi;
--- a/libfaad/mdct.h
+++ b/libfaad/mdct.h
@@ -16,7 +16,7 @@
** along with this program; if not, write to the Free Software
** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
**
-** $Id: mdct.h,v 1.8 2002/04/23 21:08:26 menno Exp $
+** $Id: mdct.h,v 1.9 2002/07/14 19:11:12 menno Exp $
**/
#ifndef __MDCT_H__
@@ -26,7 +26,11 @@
extern "C" {#endif
+#ifdef USE_FFTW
#include <fftw.h>
+#else
+#include "cfft.h"
+#endif
typedef struct {real_t sin;
@@ -33,13 +37,26 @@
real_t cos;
} faad_sincos;
+#ifndef USE_FFTW
typedef struct {+ real_t re;
+ real_t im;
+} faad_complex;
+#endif
+
+typedef struct {+ faad_sincos *sincos;
+#ifdef USE_FFTW
fftw_complex *Z1;
fftw_complex *Z2;
- faad_sincos *sincos;
fftw_plan plan_backward;
#ifdef LTP_DEC
fftw_plan plan_forward;
+#endif
+#else
+ real_t *Z1;
+ faad_complex *Z2;
+ cfft_info cfft;
#endif
uint16_t N;
} mdct_info;