summaryrefslogtreecommitdiff
path: root/libm/double/fftr.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/fftr.c')
-rw-r--r--libm/double/fftr.c237
1 files changed, 0 insertions, 237 deletions
diff --git a/libm/double/fftr.c b/libm/double/fftr.c
deleted file mode 100644
index d4ce23463..000000000
--- a/libm/double/fftr.c
+++ /dev/null
@@ -1,237 +0,0 @@
-/* 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 );
-}