aeusemrs@csun.UUCP (Mike Stump) (07/29/87)
All right, here it is!
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
# fft.doc
# fft.c
# This archive created: Sun Jul 26 11:32:35 1987
export PATH; PATH=/bin:$PATH
if test -f 'fft.doc'
then
echo shar: will not over-write existing file "'fft.doc'"
else
cat << \SHAR_EOF > 'fft.doc'
Documentation for fft.c
This program is a direct translation from FORTRAN IV of E. O. Brigham's
program on the fast Fourier transform, The Fast Fourier Transform. I suggest
reading the book, to get the details of the method.
In translating this program from FORTRAN, I retained the bit swapping
function described by Brigham. Using bit operators, I assume that a more
efficient manner of bit swapping can be performed.
I used pointers to the data to be transformed to allow the routine to
be more general. Thus fft() needs to know the number of data points when it is
called which is only limited by memory size, which can get used up quickly
since each point is represented by two 8-byte numbers. A solution might be to
hold the data on disk and use disk access functions, if you run out of memory
for the data. This method should be avoided for small data sets, since the
calculation time increases dramatically with more data points to transform.
With disk overhead, calculation time will slow down considerably. To save
memory space you can redefine the variable types from double to single
precision at the expense of precision.
I have compared the performance of the fft to a regular Fourier trans-
form, and the fft wins out especially when there are more than 512 data points.
The restriction of having the data be a multiple of 2 is usually not prohibi-
tive since the data usually are generated by a computer in base 2, eg. radio
astronomy (my original purpose for writing this program), image analysis, etc.
Enclosed is some sample data and its transform. Use this to debug the
program if necessary ( I hate to receive a data analysis program with no data
to test the program ). These data have been checked for correctness by using
a standard FFT from a mainframe computer math library.
Finally if you have any questions or suggestions, please contact me:
Jim Pisano
P.O. Box 3134
Charlottesville, VA 22903
Or contact aeusemrs@csun.uucp (Mike Stump).
INPUT OUTPUT
index real imaginary real imaginary
_____________________________________________________________________________
1 -28.66 0.0 -00.11 0.0
2 -32.43 0.0 -29.98644 -38.211
3 -8.42 0.0 -62.32944 -23.91025
4 40.22 0.0 -58.32872 -10.02253
5 1.14 0.0 -58.46804 -7.78748
6 -25.05 0.0 -68.11068 -5.28995
7 12.14 0.0 -86.80929 -11.77746
8 27.82 0.0 -137.68 -0.06883
9 2.4 0.0 -33.03304 202.24
10 -19.68 0.0 82.43765 72.43995
11 7.55 0.0 25.90461 18.1517
12 17.22 0.0 8.26912 6.555
13 -0.58 0.0 -4.33195 3.91253
14 -12.5 0.0 -8.54575 -11.5672
15 0.89 0.0 -20.32586 -3.86289
16 10.45 0.0 -8.05521 10.74816
17 -1.16 0.0 1.77 0.0
18 -7.06 0.0 -8.0552 -10.74816
19 -0.23 0.0 -20.32586 3.86289
20 4.59 0.0 -8.54575 11.5672
21 1.45 0.0 -4.33195 -3.91253
22 -3.4 0.0 8.26913 -6.55501
23 3.03 0.0 25.9046 -18.15169
24 3.05 0.0 82.43763 -72.43996
25 3.67 0.0 -33.02996 -202.24
26 -0.51 0.0 -137.68 0.06884
27 -4.35 0.0 -86.80931 11.77746
28 -1.45 0.0 -68.11066 5.28995
29 5.64 0.0 -58.46805 7.78747
30 -0.96 0.0 -58.32872 10.02253
31 -4.46 0.0 -62.32945 23.91205
32 -1.25 0.0 -29.98643 38.21101
SHAR_EOF
fi # end of overwriting check
if test -f 'fft.c'
then
echo shar: will not over-write existing file "'fft.c'"
else
cat << \SHAR_EOF > 'fft.c'
/*
Introduced into UseNet By: aeusemrs@csun.uucp (Mike Stump) Sun Jul 26 11:26:31 PDT 1987
HEADER:
TITLE: Fast Fourier Transform;
DATE: 05/18/1985;
DESCRIPTION: "Performs fast fourier transform using method described
by E. O. Brigham. For details of the method, refer
to Brigham's book. THE FAST FOURIER TRANSFORM";
KEYWORDS: Fourier, transform;
FILENAME: FFT.C;
WARNINGS:
"This program is self-contained, all that is needed is a manner of getting
the data into the array real_data (& imag_data, if applicable). The
transformed data will reside in these two arrays upon return with the
original data being destroyed."
AUTHORS: Jim Pisano;
COMPILERS: DeSmet;
REFERENCES: AUTHORS: E. O. Brigham;
TITLE: "THE FAST FOURIER TRANSFORM";
CITATION: "";
ENDREF
*/
/* file name fft.c
* program name fft() ... Fast Fourier Transform
*
* Perform fast fourier transform using method described by E. O. Brigham.
* For details of the method, refer to Brigham's book
*
* Translated to C from FORTRAN by
*
* Jim Pisano
* P.O. Box 3134
* University Station
* Charlottesville, VA 22903
*
* This program is in the public domain & may be used by anyone for commercial
* or non-commercial purposes.
*
* real_data ... ptr. to real part of data to be transformed
* imag_data ... ptr. to imag " " " " " "
* inv ..... Switch to flag normal or inverse transform
* n_pts ... Number of real data points
* nu ...... logarithm in base 2 of n_pts e.g. nu = 5 if n_pts = 32.
*
* This program is self-contained, all that is needed is a manner of getting
* the data into the array real_data (& imag_data, if applicable). The
* transformed data will reside in these two arrays upon return with the
* original data being destroyed.
*/
#include <stdio.h>
#include <math.h> /* declare math functions to use */
#define PI 3.1419527
fft( real_data, imag_data, n_pts, nu, inv )
double *real_data, *imag_data;
int n_pts, nu, inv;
{
int n2, j, j1, l, i, ib, k, k1, k2;
int sgn;
double tr, ti, arg, nu1; /* intermediate values in calcs. */
double c, s; /* cosine & sine components of Fourier trans. */
n2 = n_pts / 2;
nu1 = nu - 1.0;
k = 0;
/*
* sign change for inverse transform
*/
sgn = inv ? -1 : 1;
/*
* Calculate the componets of the Fourier series of the function
*/
for( l = 0; l != nu; l++ )
{
do
{
for( i = 0; i != n2; i++ )
{
j = k / ( pow( 2.0, nu1 ) );
ib = bit_swap( j, nu );
arg = 2.0 * PI * ib / n_pts;
c = cos( arg );
s = sgn * sin( arg );
k1 = k;
k2 = k1 + n2;
tr = *(real_data+k2) * c + *(imag_data+k2) * s;
ti = *(imag_data+k2) * c - *(real_data+k2) * s;
*(real_data+k2) = *(real_data+k1) - tr;
*(imag_data+k2) = *(imag_data+k1) - ti;
*(real_data+k1) = *(real_data+k1) + tr;
*(imag_data+k1) = *(imag_data+k1) + ti;
k++;
}
k += n2;
} while( k < n_pts - 1);
k = 0;
nu1 -= 1.0;
n2 /= 2;
}
for( k = 0; k != n_pts; k++ )
{
ib = bit_swap( k, nu );
if( ib > k)
{
swap( (real_data+k), (real_data+ib) );
swap( (imag_data+k), (imag_data+ib) );
}
}
/*
* If calculating the inverse transform, must divide the data by the number of
* data points.
*/
if( inv )
for( k = 0; k != n_pts; k++)
{
*(real_data+k) /= n_pts;
*(imag_data+k) /= n_pts;
}
}
/*
* Bit swaping routine in which the bit pattern of the integer i is reordered.
* See Brigham's book for details
*/
bit_swap( i, nu )
int i, nu;
{
int ib, i1, i2;
ib = 0;
for( i1 = 0; i1 != nu; i1++ )
{
i2 = i / 2;
ib = ib * 2 + i - 2 * i2;
i = i2;
}
return( ib );
}
/*
* Simple exchange routine where *x1 & *x2 are swapped
*/
swap( x1, x2 )
double *x1, *x2;
{
int *temp_iy;
double *temp_x;
*temp_x = *x1;
*x1 = *x2;
*x2 = *temp_x;
}
SHAR_EOF
fi # end of overwriting check
# End of shell archive
exit 0
--
Mike Stump, Cal State Univ, Northridge Comp Sci Department
uucp: {sdcrdcf, ihnp4, hplabs, ttidca, psivax, csustan}!csun!aeusemrs