dtynan@sultra.UUCP (Der Tynan) (10/01/88)
Hi!
I received enough responses (and enough bounced mail!), that I'm hereby
posting my Fast Fourier Transform (FFT) program. Do with it what you
want. I don't really care. However, if you really *must* send huge
amounts of money to ease your conscience, my postal address is in the
maps :-) Just kidding. It's free. A word of caution, however. This
program was NOT written to be public domain. It was written for a client,
to check the fidelity of a real-time FFT system. For this reason, it
focuses more on high fidelity, than it does on speed, or portability. I'd
appreciate peoples comments on how they got on with this thing. What little
inline documentation there is, is poor - like I said, it was written to
serve one unique purpose (which it did!). The TWIDDLE factors are computed
at the beginning of the program, to save some time. Also, the maximum
transform size is extremely limited (on this system it only goes to 8K).
The reason for this, is that it uses Singletons algorithm, which is an
out-of-place algorithm (twice as much memory). Also, each data sample
(complex) is 'double', plus the storage for the TWIDDLE factors...
If you don't have complex data (and don't know what a HILBERT transform is),
then just zero out the imaginary component (this causes aliasing).
Have fun!
- Der
-----cut here-----cut here-----cut here-----cut here-----
#! /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 to create the files:
# compute.c
# This archive created: Thu Sep 29 13:36:12 1988
export PATH; PATH=/bin:$PATH
echo shar: extracting "'compute.c'" '(5372 characters)'
if test -f 'compute.c'
then
echo shar: will not over-write existing file "'compute.c'"
else
sed 's/^X//' << \SHAR_EOF > 'compute.c'
X#ifndef lint
Xstatic char *sccsid = "@(#)compute.c 1.0 4/10/87";
X#endif
X
X/*
X * This program is Public Domain. It may be freely distributed, both as
X * source- and object-code. No warranties, either implied or expressed
X * are provided with this code. It is purely on an "AS-IS" basis.
X * Inclusion of this source (or related object) in any commercial product
X * is permitted. This notice must remain intact, however.
X *
X * Copyright (C) 1987, 1988, Tynan Computers. All rights reserved.
X */
X
X#include <stdio.h>
X#include <math.h>
X
X#define MAXITER 13
X#define MAXXFORM 8192
X#define OUTPUT "fft.out"
X
Xstruct sample {
X double s_real; /* REAL value */
X double s_imag; /* IMAGINARY value */
X};
X
Xint nsample; /* No. of samples */
Xint niter; /* No. of iterations */
Xint amask; /* sin() address mask */
Xint pageno;
Xdouble xsin(), xcos(); /* Local functions */
Xdouble *twp; /* Pointer to twiddle */
Xstruct sample samplea[MAXXFORM], *dsrc; /* Actual dataset (a) */
Xstruct sample sampleb[MAXXFORM], *ddst; /* Actual dataset (b) */
X
X/*
X * ***** FFT Computation Program *****
X *
X * Written by: Dermot Tynan, 10-Apr-87 (dtynan@zorba.uucp)
X *
X * This program uses Singletons Algorithm, Radix-2, Out-of-Place, Decimation-
X * in-Time.
X *
X * Compile with; cc -O -o compute compute.c -lm
X *
X */
Xmain(argc, argv)
Xchar *argv[];
X{
X register i, j;
X register struct sample *sp;
X double twd;
X double xm;
X FILE *fp;
X
X /*
X * After each iteration, the source and destination banks are
X * swapped. Hence, two pointers are used, which alternate
X * between banks of data.
X */
X dsrc = samplea;
X ddst = sampleb;
X printf("FFT Computation Program - Written by Dermot Tynan.\n\n");
X if (argc != 2) {
X fprintf(stderr, "usage: compute filename\n");
X exit(2);
X }
X if ((fp = fopen(argv[1], "r")) == NULL) {
X perror(argv[1]);
X exit(1);
X }
X printf("Input phase...\n");
X for (i = 0, sp = dsrc; i < MAXXFORM; i++, sp++) {
X if (feof(fp))
X break;
X fscanf(fp, "%lf %lf\n", &sp->s_real, &sp->s_imag);
X }
X fclose(fp);
X nsample = i;
X for (i = 0, niter = 1; i < MAXITER; i++) {
X if (niter == nsample)
X break;
X if (niter > nsample) {
X fprintf(stderr, "?Invalid # of samples.\n");
X exit(1);
X }
X niter <<= 1;
X }
X niter = i;
X printf("[%d samples - %d iterations]\n", nsample, niter);
X /*
X * As an optimization, instead of constantly re-referencing the
X * sin() and cosine(), we generate a table with all the values
X * we will use. For example, a 64pt FFT only uses 64 discrete
X * sine angles. Some of these, however, are accessed up to six
X * times. A table is thus built, holding all such values.
X */
X printf("Twiddle generation phase...\n");
X if ((twp = (double *)malloc(nsample * sizeof(double))) == NULL) {
X fprintf(stderr, "?Error - insufficient memory.\n");
X exit(1);
X }
X for (i = 0; i < nsample; i++) {
X twd = (2.0 * M_PI * (double )i)/(double )nsample;
X twp[i] = sin(twd);
X }
X /*
X * Open the output file now, just in case there's a problem.
X */
X if ((fp = fopen(OUTPUT, "w")) == NULL) {
X perror(OUTPUT);
X exit(1);
X }
X printf("Transform phase...\n");
X amask = 0;
X for (i = 0; i < niter; i++) {
X printf("--- iteration #%d\n", i + 1);
X compute(); /* compute an iter. */
X amask <<= 1; /* Open up the mask */
X amask |= 1;
X sp = dsrc; /* Swap banks */
X dsrc = ddst;
X ddst = sp;
X }
X printf("Output phase...\n");
X for (i = 0; i < nsample; i++) {
X /*
X * Do bit-reversal, and output.
X */
X sp = &dsrc[bitrev(i)];
X fprintf(fp, "%26.15lf\t%26.15lf\n", sp->s_real, sp->s_imag);
X }
X printf("Done!\n");
X fclose(fp);
X exit(0);
X}
X
X/*
X * This routine performs one iteration of an FFT. The data is read from
X * dsrc, through the butterfly computation, and written to ddst (this is
X * singletons out-of-place algorithm).
X */
Xcompute()
X{
X register n, k, th;
X register struct sample *sp;
X double ar, ai;
X double br, bi;
X double wr, wi;
X
X k = nsample / 2; /* k = N/2 */
X for (n = 0; n < k; n++) {
X /*
X * Read x(n) and x(n + N/2)
X * x(n) = A(n) = (ar, ai).
X * x(n + N/2) = B(n) = (br, bi).
X */
X sp = &dsrc[n];
X ar = sp->s_real;
X ai = sp->s_imag;
X sp = &dsrc[n+k];
X br = sp->s_real;
X bi = sp->s_imag;
X /*
X * Calculate TWIDDLE factors.
X */
X th = bitrev(n & amask) >> 1;
X wr = xcos(th);
X wi = xsin(th);
X /*
X * Do butterfly calculation.
X * A'(n) = A + BWk
X * B'(n) = A - BWk
X */
X sp = &ddst[n << 1];
X sp->s_real = ar + br * wr + bi * wi;
X sp->s_imag = ai + bi * wr - br * wi;
X sp++;
X sp->s_real = ar - br * wr - bi * wi;
X sp->s_imag = ai - bi * wr + br * wi;
X }
X}
X
X/*
X * These two functions replicate the actual math library routines. They
X * use an internal memory image, instead of computing the actual values.
X */
Xdouble
Xxcos(th)
Xregister th;
X{
X th += (nsample / 4);
X return(xsin(th));
X}
X
Xdouble
Xxsin(th)
Xregister th;
X{
X register i, j;
X th &= (nsample - 1);
X return(twp[th]);
X}
X
X/*
X * In the FFT world, bit reversals are a fact of life. For example, assume
X * a seven bit word: 0x65. In binary, this is 1100101. When bit reversed,
X * the result is 1010011 (see the difference?). Giving a result of 0x53.
X * This is further complicated by the fact that the no. of bits is a function
X * of the transform size.
X */
Xbitrev(n)
Xregister n;
X{
X register i, r;
X
X for (i = r = 0; i < niter; i++, r <<= 1, n >>= 1)
X if (n & 1)
X r |= 1;
X r >>= 1; /* We shifted once too often */
X return(r);
X}
SHAR_EOF
if test 5372 -ne "`wc -c < 'compute.c'`"
then
echo shar: error transmitting "'compute.c'" '(should have been 5372 characters)'
fi
fi # end of overwriting check
# End of shell archive
exit 0
--
Reply: dtynan@sultra.UUCP (Der Tynan @ Tynan Computers)
{mips,pyramid}!sultra!dtynan
Cast a cold eye on life, on death. Horseman, pass by... [WBY]