diff options
Diffstat (limited to 'libm/double/fftr.c')
| -rw-r--r-- | libm/double/fftr.c | 237 |
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 ); -} |