Index: /dsfilters/temp/libmpcenc/src/fft4g.c
===================================================================
--- /dsfilters/temp/libmpcenc/src/fft4g.c	(revision 396)
+++ /dsfilters/temp/libmpcenc/src/fft4g.c	(revision 396)
@@ -0,0 +1,525 @@
+/*
+ * Musepack audio compression
+ * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * This library 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+
+
+/* F U N C T I O N S */
+static          void  makewt       ( const int nw, int* ip, float* w );
+static          void  makect       ( const int nc, int* ip, float* c );
+
+# define cftmdl(n,l,a,w)   cftmdl_i386  ( n, l, a, w )
+
+
+/* -------- child routines -------- */
+static mpc_inline void
+bitrv2 ( const int n, int* ip, float* a )
+{
+	int    j, j1, k, k1, l, m, m2;
+	float  xr, xi, yr, yi;
+
+	ip[0] = 0;
+	l     = n;
+	m     = 1;
+	while ( (m << 3) < l ) {
+		l >>= 1;
+		for ( j = 0; j < m; j++ ) {
+			ip[m + j] = ip[j] + l;
+		}
+		m <<= 1;
+	}
+	m2 = 2 * m;
+	if ( (m << 3) == l ) {
+		for ( k = 0; k < m; k++ ) {
+			for ( j = 0; j < k; j++ ) {
+				j1        = 2 * j + ip[k];
+				k1        = 2 * k + ip[j];
+				xr        = a[j1];
+				xi        = a[j1 + 1];
+				yr        = a[k1];
+				yi        = a[k1 + 1];
+				a[j1]     = yr;
+				a[j1 + 1] = yi;
+				a[k1]     = xr;
+				a[k1 + 1] = xi;
+				j1       += m2;
+				k1       += 2 * m2;
+				xr        = a[j1];
+				xi        = a[j1 + 1];
+				yr        = a[k1];
+				yi        = a[k1 + 1];
+				a[j1]     = yr;
+				a[j1 + 1] = yi;
+				a[k1]     = xr;
+				a[k1 + 1] = xi;
+				j1       += m2;
+				k1       -= m2;
+				xr        = a[j1];
+				xi        = a[j1 + 1];
+				yr        = a[k1];
+				yi        = a[k1 + 1];
+				a[j1]     = yr;
+				a[j1 + 1] = yi;
+				a[k1]     = xr;
+				a[k1 + 1] = xi;
+				j1       += m2;
+				k1       += 2 * m2;
+				xr        = a[j1];
+				xi        = a[j1 + 1];
+				yr        = a[k1];
+				yi        = a[k1 + 1];
+				a[j1]     = yr;
+				a[j1 + 1] = yi;
+				a[k1]     = xr;
+				a[k1 + 1] = xi;
+			}
+			j1        = 2 * k + m2 + ip[k];
+			k1        = j1 + m2;
+			xr        = a[j1];
+			xi        = a[j1 + 1];
+			yr        = a[k1];
+			yi        = a[k1 + 1];
+			a[j1]     = yr;
+			a[j1 + 1] = yi;
+			a[k1]     = xr;
+			a[k1 + 1] = xi;
+		}
+	} else {
+		for ( k = 1; k < m; k++ ) {
+			for ( j = 0; j < k; j++ ) {
+				j1        = 2 * j + ip[k];
+				k1        = 2 * k + ip[j];
+				xr        = a[j1];
+				xi        = a[j1 + 1];
+				yr        = a[k1];
+				yi        = a[k1 + 1];
+				a[j1]     = yr;
+				a[j1 + 1] = yi;
+				a[k1]     = xr;
+				a[k1 + 1] = xi;
+				j1       += m2;
+				k1       += m2;
+				xr        = a[j1];
+				xi        = a[j1 + 1];
+				yr        = a[k1];
+				yi        = a[k1 + 1];
+				a[j1]     = yr;
+				a[j1 + 1] = yi;
+				a[k1]     = xr;
+				a[k1 + 1] = xi;
+			}
+		}
+	}
+	return;
+}
+
+static mpc_inline void
+cft1st ( const int n, float* a, float* w )
+{
+	int    j, k1;
+	float  wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
+	float  x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
+
+	x0r   = a[ 0] + a[ 2];
+	x0i   = a[ 1] + a[ 3];
+	x1r   = a[ 0] - a[ 2];
+	x1i   = a[ 1] - a[ 3];
+	x2r   = a[ 4] + a[ 6];
+	x2i   = a[ 5] + a[ 7];
+	x3r   = a[ 4] - a[ 6];
+	x3i   = a[ 5] - a[ 7];
+	a[ 0] = x0r + x2r;
+	a[ 1] = x0i + x2i;
+	a[ 4] = x0r - x2r;
+	a[ 5] = x0i - x2i;
+	a[ 2] = x1r - x3i;
+	a[ 3] = x1i + x3r;
+	a[ 6] = x1r + x3i;
+	a[ 7] = x1i - x3r;
+	wk1r  = w[ 2];
+	x0r   = a[ 8] + a[10];
+	x0i   = a[ 9] + a[11];
+	x1r   = a[ 8] - a[10];
+	x1i   = a[ 9] - a[11];
+	x2r   = a[12] + a[14];
+	x2i   = a[13] + a[15];
+	x3r   = a[12] - a[14];
+	x3i   = a[13] - a[15];
+	a[ 8] = x0r + x2r;
+	a[ 9] = x0i + x2i;
+	a[12] = x2i - x0i;
+	a[13] = x0r - x2r;
+	x0r   = x1r - x3i;
+	x0i   = x1i + x3r;
+	a[10] = wk1r * (x0r - x0i);
+	a[11] = wk1r * (x0r + x0i);
+	x0r   = x3i + x1r;
+	x0i   = x3r - x1i;
+	a[14] = wk1r * (x0i - x0r);
+	a[15] = wk1r * (x0i + x0r);
+
+	k1 = 0;
+	j  = 16;
+	do {
+		k1       += 2;
+		wk2r      = w[k1];
+		wk2i      = w[k1 + 1];
+		wk1r      = w[2*k1];
+		wk1i      = w[2*k1 + 1];
+		wk3r      = wk1r - 2 * wk2i * wk1i;
+		wk3i      = 2 * wk2i * wk1r - wk1i;
+		x0r       = a[j]     + a[j + 2];
+		x0i       = a[j + 1] + a[j + 3];
+		x1r       = a[j]     - a[j + 2];
+		x1i       = a[j + 1] - a[j + 3];
+		x2r       = a[j + 4] + a[j + 6];
+		x2i       = a[j + 5] + a[j + 7];
+		x3r       = a[j + 4] - a[j + 6];
+		x3i       = a[j + 5] - a[j + 7];
+		a[j]      = x0r + x2r;
+		a[j + 1]  = x0i + x2i;
+		x0r      -= x2r;
+		x0i      -= x2i;
+		a[j + 4]  = wk2r * x0r - wk2i * x0i;
+		a[j + 5]  = wk2r * x0i + wk2i * x0r;
+		x0r       = x1r - x3i;
+		x0i       = x1i + x3r;
+		a[j + 2]  = wk1r * x0r - wk1i * x0i;
+		a[j + 3]  = wk1r * x0i + wk1i * x0r;
+		x0r       = x1r + x3i;
+		x0i       = x1i - x3r;
+		a[j + 6]  = wk3r * x0r - wk3i * x0i;
+		a[j + 7]  = wk3r * x0i + wk3i * x0r;
+		wk1r      = w[2*k1 + 2];
+		wk1i      = w[2*k1 + 3];
+		wk3r      = wk1r - 2 * wk2r * wk1i;
+		wk3i      = 2 * wk2r * wk1r - wk1i;
+		x0r       = a[j +  8] + a[j + 10];
+		x0i       = a[j +  9] + a[j + 11];
+		x1r       = a[j +  8] - a[j + 10];
+		x1i       = a[j +  9] - a[j + 11];
+		x2r       = a[j + 12] + a[j + 14];
+		x2i       = a[j + 13] + a[j + 15];
+		x3r       = a[j + 12] - a[j + 14];
+		x3i       = a[j + 13] - a[j + 15];
+		a[j + 8]  = x0r + x2r;
+		a[j + 9]  = x0i + x2i;
+		x0r      -= x2r;
+		x0i      -= x2i;
+		a[j + 12] = -wk2i * x0r - wk2r * x0i;
+		a[j + 13] = -wk2i * x0i + wk2r * x0r;
+		x0r       = x1r - x3i;
+		x0i       = x1i + x3r;
+		a[j + 10] = wk1r * x0r - wk1i * x0i;
+		a[j + 11] = wk1r * x0i + wk1i * x0r;
+		x0r       = x1r + x3i;
+		x0i       = x1i - x3r;
+		a[j + 14] = wk3r * x0r - wk3i * x0i;
+		a[j + 15] = wk3r * x0i + wk3i * x0r;
+	} while ( j += 16, j < n );
+	return;
+}
+
+static mpc_inline void
+cftmdl_i386 ( const int n, const int l, float* a, float* w )
+{
+	int    j, j1, j2, j3, k, k1, m, m2;
+	float  wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
+	float  x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
+
+	m = l << 2;
+
+	for ( j = 0; j < l; j += 2 ) {
+		j1        = j  + l;
+		j2        = j1 + l;
+		j3        = j2 + l;
+		x0r       = a[j]      + a[j1];
+		x0i       = a[j + 1]  + a[j1 + 1];
+		x1r       = a[j]      - a[j1];
+		x1i       = a[j + 1]  - a[j1 + 1];
+		x2r       = a[j2]     + a[j3];
+		x2i       = a[j2 + 1] + a[j3 + 1];
+		x3r       = a[j2]     - a[j3];
+		x3i       = a[j2 + 1] - a[j3 + 1];
+		a[j]      = x0r + x2r;
+		a[j + 1]  = x0i + x2i;
+		a[j2]     = x0r - x2r;
+		a[j2 + 1] = x0i - x2i;
+		a[j1]     = x1r - x3i;
+		a[j1 + 1] = x1i + x3r;
+		a[j3]     = x1r + x3i;
+		a[j3 + 1] = x1i - x3r;
+	}
+
+	wk1r = w[2];
+	for ( j = m; j < l + m; j += 2 ) {
+		j1        = j  + l;
+		j2        = j1 + l;
+		j3        = j2 + l;
+		x0r       = a[j]      + a[j1];
+		x0i       = a[j + 1]  + a[j1 + 1];
+		x1r       = a[j]      - a[j1];
+		x1i       = a[j + 1]  - a[j1 + 1];
+		x2r       = a[j2]     + a[j3];
+		x2i       = a[j2 + 1] + a[j3 + 1];
+		x3r       = a[j2]     - a[j3];
+		x3i       = a[j2 + 1] - a[j3 + 1];
+		a[j]      = x0r + x2r;
+		a[j + 1]  = x0i + x2i;
+		a[j2]     = x2i - x0i;
+		a[j2 + 1] = x0r - x2r;
+		x0r       = x1r - x3i;
+		x0i       = x1i + x3r;
+		a[j1]     = wk1r * (x0r - x0i);
+		a[j1 + 1] = wk1r * (x0r + x0i);
+		x0r       = x3i + x1r;
+		x0i       = x3r - x1i;
+		a[j3]     = wk1r * (x0i - x0r);
+		a[j3 + 1] = wk1r * (x0i + x0r);
+	}
+
+	k1 = 0;
+	m2 = 2 * m;
+	for ( k = m2; k < n; k += m2 ) {
+		k1  += 2;
+		wk2r = w[k1];
+		wk2i = w[k1 + 1];
+		wk1r = w[2*k1];
+		wk1i = w[2*k1 + 1];
+		wk3r = wk1r - 2 * wk2i * wk1i;
+		wk3i = 2 * wk2i * wk1r - wk1i;
+		j    = k;
+		do {
+			j1        = j  + l;
+			j2        = j1 + l;
+			j3        = j2 + l;
+			x0r       = a[j]      + a[j1];
+			x0i       = a[j + 1]  + a[j1 + 1];
+			x1r       = a[j]      - a[j1];
+			x1i       = a[j + 1]  - a[j1 + 1];
+			x2r       = a[j2]     + a[j3];
+			x2i       = a[j2 + 1] + a[j3 + 1];
+			x3r       = a[j2]     - a[j3];
+			x3i       = a[j2 + 1] - a[j3 + 1];
+			a[j]      = x0r + x2r;
+			a[j + 1]  = x0i + x2i;
+			x0r      -= x2r;
+			x0i      -= x2i;
+			a[j2]     = wk2r * x0r - wk2i * x0i;
+			a[j2 + 1] = wk2r * x0i + wk2i * x0r;
+			x0r       = x1r - x3i;
+			x0i       = x1i + x3r;
+			a[j1]     = wk1r * x0r - wk1i * x0i;
+			a[j1 + 1] = wk1r * x0i + wk1i * x0r;
+			x0r       = x1r + x3i;
+			x0i       = x1i - x3r;
+			a[j3]     = wk3r * x0r - wk3i * x0i;
+			a[j3 + 1] = wk3r * x0i + wk3i * x0r;
+		} while ( j += 2, j < l + k );
+
+		wk1r = w[2*k1 + 2];
+		wk1i = w[2*k1 + 3];
+		wk3r = wk1r - 2 * wk2r * wk1i;
+		wk3i = 2 * wk2r * wk1r - wk1i;
+		j    = k + m;
+		do {
+			j1        = j  + l;
+			j2        = j1 + l;
+			j3        = j2 + l;
+			x0r       = a[j]      + a[j1];
+			x0i       = a[j + 1]  + a[j1 + 1];
+			x1r       = a[j]      - a[j1];
+			x1i       = a[j + 1]  - a[j1 + 1];
+			x2r       = a[j2]     + a[j3];
+			x2i       = a[j2 + 1] + a[j3 + 1];
+			x3r       = a[j2]     - a[j3];
+			x3i       = a[j2 + 1] - a[j3 + 1];
+			a[j]      = x0r + x2r;
+			a[j + 1]  = x0i + x2i;
+			x0r      -= x2r;
+			x0i      -= x2i;
+			a[j2]     = -wk2i * x0r - wk2r * x0i;
+			a[j2 + 1] = -wk2i * x0i + wk2r * x0r;
+			x0r       = x1r - x3i;
+			x0i       = x1i + x3r;
+			a[j1]     = wk1r * x0r - wk1i * x0i;
+			a[j1 + 1] = wk1r * x0i + wk1i * x0r;
+			x0r       = x1r + x3i;
+			x0i       = x1i - x3r;
+			a[j3]     = wk3r * x0r - wk3i * x0i;
+			a[j3 + 1] = wk3r * x0i + wk3i * x0r;
+		} while ( j += 2, j < l+k+m );
+	}
+	return;
+}
+
+static mpc_inline void
+cftfsub ( const int n, float* a, float* w )
+{
+	int    j, j1, j2, j3, l;
+	float  x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
+
+	l = 2;
+	if ( n > 8 ) {
+		cft1st ( n, a, w );
+		l = 8;
+		while ( (l << 2) < n ) {
+			cftmdl ( n, l, a, w );
+			l <<= 2;
+		}
+	}
+	if ( (l << 2) == n ) {
+		j = 0;
+		do {
+			j1        = j  + l;
+			j2        = j1 + l;
+			j3        = j2 + l;
+			x0r       = a[j]      + a[j1];
+			x0i       = a[j + 1]  + a[j1 + 1];
+			x1r       = a[j]      - a[j1];
+			x1i       = a[j + 1]  - a[j1 + 1];
+			x2r       = a[j2]     + a[j3];
+			x2i       = a[j2 + 1] + a[j3 + 1];
+			x3r       = a[j2]     - a[j3];
+			x3i       = a[j2 + 1] - a[j3 + 1];
+			a[j]      = x0r + x2r;
+			a[j + 1]  = x0i + x2i;
+			a[j2]     = x0r - x2r;
+			a[j2 + 1] = x0i - x2i;
+			a[j1]     = x1r - x3i;
+			a[j1 + 1] = x1i + x3r;
+			a[j3]     = x1r + x3i;
+			a[j3 + 1] = x1i - x3r;
+		} while ( j += 2, j < l );
+	} else {
+		j = 0;
+		do {
+			j1        = j + l;
+			x0r       = a[j]     - a[j1];
+			x0i       = a[j + 1] - a[j1 + 1];
+			a[j]     += a[j1];
+			a[j + 1] += a[j1 + 1];
+			a[j1]     = x0r;
+			a[j1 + 1] = x0i;
+		} while ( j += 2, j < l );
+	}
+	return;
+}
+
+
+static mpc_inline void
+rftfsub ( const int n, float* a, int nc, float* c )
+{
+	int    j, k, kk, ks, m;
+	float  wkr, wki, xr, xi, yr, yi;
+
+	m  = n >> 1;
+	ks = 2 * nc / m;
+	kk = ks;
+	j  = 2;
+	k  = n;
+	do {
+		k        -= 2;
+		nc       -= ks;
+		wkr       = 0.5f - c[nc];
+		wki       = c[kk];
+		xr        = a[j]     - a[k];
+		xi        = a[j + 1] + a[k + 1];
+		yr        = wkr * xr - wki * xi;
+		yi        = wkr * xi + wki * xr;
+		a[j]     -= yr;
+		a[j + 1] -= yi;
+		a[k]     += yr;
+		a[k + 1] -= yi;
+		kk       += ks;
+	} while ( j += 2, j < m );
+	return;
+}
+
+/* -------- initializing routines -------- */
+static void
+makewt ( const int nw, int* ip, float* w )
+{
+    int     j;
+    int     nwh;
+    float   x;
+    float   y;
+    double  delta;
+
+    ip[0] = nw;
+    ip[1] = 1;
+    if ( nw > 2 ) {
+        nwh        = nw >> 1;
+        delta      = (M_PI/4) / nwh;
+        w[0]       = 1.;
+        w[1]       = 0.;
+        w[nwh]     = COS (delta * nwh);
+        w[nwh + 1] = w[nwh];
+        if ( nwh > 2 ) {
+            for ( j = 2; j < nwh; j += 2 ) {
+                x             = COS (delta * j);
+                y             = SIN (delta * j);
+                w[j]          = x;
+                w[j + 1]      = y;
+                w[nw - j]     = y;
+                w[nw - j + 1] = x;
+            }
+            bitrv2 ( nw, ip + 2, w );
+        }
+    }
+    return;
+}
+
+
+static void
+makect ( const int nc, int* ip, float* c )
+{
+    int     j;
+    int     nch;
+    double  delta;
+
+    ip[1] = nc;
+    if ( nc > 1 ) {
+        nch    = nc >> 1;
+        delta  = (M_PI/4) / nch;
+        c[0]   = COS (delta * nch);
+        c[nch] = 0.5f * c[0];
+        for ( j = 1; j < nch; j++ ) {
+            c[j]      = 0.5f * COS (delta * j);
+            c[nc - j] = 0.5f * SIN (delta * j);
+        }
+    }
+    return;
+}
+
+void mpc_psy_model_generate_fft_tables(const int n, int *ip, float *w)
+{
+    int  nw;
+    int  nc;
+
+    nw = n >> 2;
+    makewt ( nw, ip, w );
+
+    nc = n >> 2;
+    makect ( nc, ip, w + nw );
+}
+
+
+
+
+