uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 237 insertions, 0 deletions

diff --git a/libm/double/fftr.c b/libm/double/fftr.c
new file mode 100644
index 000000000..d4ce23463
--- /dev/null
+++ b/libm/double/fftr.c
@@ -0,0 +1,237 @@
+/*							fftr.c
+ *
+ *	FFT of Real Valued Sequence
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x[], sine[];
+ * int m;
+ *
+ * fftr( x, m, sine );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Computes the (complex valued) discrete Fourier transform of
+ * the real valued sequence x[].  The input sequence x[] contains
+ * n = 2**m samples.  The program fills array sine[k] with
+ * n/4 + 1 values of sin( 2 PI k / n ).
+ *
+ * Data format for complex valued output is real part followed
+ * by imaginary part.  The output is developed in the input
+ * array x[].
+ *
+ * The algorithm takes advantage of the fact that the FFT of an
+ * n point real sequence can be obtained from an n/2 point
+ * complex FFT.
+ *
+ * A radix 2 FFT algorithm is used.
+ *
+ * Execution time on an LSI-11/23 with floating point chip
+ * is 1.0 sec for n = 256.
+ *
+ *
+ *
+ * REFERENCE:
+ *
+ * E. Oran Brigham, The Fast Fourier Transform;
+ * Prentice-Hall, Inc., 1974
+ *
+ */
+
+
+#include <math.h>
+
+static short n0 = 0;
+static short n4 = 0;
+static short msav = 0;
+
+extern double PI;
+
+#ifdef ANSIPROT
+extern double sin ( double );
+static int bitrv(int, int);
+#else
+double sin();
+static int bitrv();
+#endif
+
+fftr( x, m0, sine )
+double x[];
+int m0;
+double sine[];
+{
+int th, nd, pth, nj, dth, m;
+int n, n2, j, k, l, r;
+double xr, xi, tr, ti, co, si;
+double a, b, c, d, bc, cs, bs, cc;
+double *p, *q;
+
+/* Array x assumed filled with real-valued data */
+/* m0 = log2(n0)			*/
+/* n0 is the number of real data samples */
+
+if( m0 != msav )
+	{
+	msav = m0;
+
+	/* Find n0 = 2**m0	*/
+	n0 = 1;
+	for( j=0; j<m0; j++ )
+		n0 <<= 1;
+
+	n4 = n0 >> 2;
+
+	/* Calculate array of sines */
+	xr = 2.0 * PI / n0;
+	for( j=0; j<=n4; j++ )
+		sine[j] = sin( j * xr );
+	}
+
+n = n0 >> 1;	/* doing half length transform */
+m = m0 - 1;
+
+
+/*							fftr.c	*/
+
+/*  Complex Fourier Transform of n Complex Data Points */
+
+/*	First, bit reverse the input data	*/
+
+for( k=0; k<n; k++ )
+	{
+	j = bitrv( k, m );
+	if( j > k )
+		{ /* executed approx. n/2 times */
+		p = &x[2*k];
+		tr = *p++;
+		ti = *p;
+		q = &x[2*j+1];
+		*p = *q;
+		*(--p) = *(--q);
+		*q++ = tr;
+		*q = ti;
+		}
+	}
+	
+/*							fftr.c	*/
+/*			Radix 2 Complex FFT			*/
+n2 = n/2;
+nj = 1;
+pth = 1;
+dth = 0;
+th = 0;
+
+for( l=0; l<m; l++ )
+	{	/* executed log2(n) times, total */
+	j = 0;
+	do
+		{	/* executed n-1 times, total */
+		r = th << 1;
+		si = sine[r];
+		co = sine[ n4 - r ];
+		if( j >= pth )
+			{
+			th -= dth;
+			co = -co;
+			}
+		else
+			th += dth;
+
+		nd = j;
+
+		do
+			{ /* executed n/2 log2(n) times, total */
+			r = (nd << 1) + (nj << 1);
+			p = &x[ r ];
+			xr = *p++;
+			xi = *p;
+			tr = xr * co + xi * si;
+			ti = xi * co - xr * si;
+			r = nd << 1;
+			q = &x[ r ];
+			xr = *q++;
+			xi = *q;
+			*p = xi - ti;
+			*(--p) = xr - tr;
+			*q = xi + ti;
+			*(--q) = xr + tr;
+			nd += nj << 1;
+			}
+		while( nd < n );
+		}
+	while( ++j < nj );
+
+	n2 >>= 1;
+	dth = n2;
+	pth = nj;
+	nj <<= 1;
+	}
+
+/*							fftr.c	*/
+
+/*	Special trick algorithm			*/
+/*	converts to spectrum of real series	*/
+
+/* Highest frequency term; add space to input array if wanted */
+/*
+x[2*n] = x[0] - x[1];
+x[2*n+1] = 0.0;
+*/
+
+/* Zero frequency term */
+x[0] = x[0] + x[1];
+x[1] = 0.0;
+n2 = n/2;
+
+for( j=1; j<=n2; j++ )
+	{	/* executed n/2 times */
+	si = sine[j];
+	co = sine[ n4 - j ];
+	p = &x[ 2*j ];
+	xr = *p++;
+	xi = *p;
+	q = &x[ 2*(n-j) ];
+	tr = *q++;
+	ti = *q;
+	a = xr + tr;
+	b = xi + ti;
+	c = xr - tr;
+	d = xi - ti;
+	bc = b * co;
+	cs = c * si;
+	bs = b * si;
+	cc = c * co;
+	*p = ( d - bs - cc )/2.0;
+	*(--p) = ( a + bc - cs )/2.0;
+	*q = -( d + bs + cc )/2.0;
+	*(--q) = ( a - bc + cs )/2.0;
+	}
+
+return(0);
+}
+
+/*							fftr.c	*/
+
+/*	Bit reverser	*/
+
+int bitrv( j, m )
+int j, m;
+{
+register int j1, ans;
+short k;
+
+ans = 0;
+j1 = j;
+
+for( k=0; k<m; k++ )
+	{
+	ans = (ans << 1) + (j1 & 1);
+	j1 >>= 1;
+	}
+
+return( ans );
+}