rsalz@uunet.uu.net (Rich Salz) (03/14/89)
Submitted-by: Peter Valkenburg <cs.vu.nl!valke>
Posting-number: Volume 18, Issue 20
Archive-name: fft
The packages included below perform fast fourier transformations on an
arbitrary number of real or complex samples. It uses a generalized version
of the well-known Cooley-Tukey algorithm, meaning that it will also work
for a number of samples that is not a power of 2.
Since I didn't bother to squeeze out the very last machine cycle, you are
not advised to use this in real-time applications. However, the recursive
algorithm used is very neat.
Enjoy.
Peter Valkenburg (valke@cs.vu.nl).
--------------------------------------------------------------------------------
WARNING: This package shows serious deficiencies if used in SDI-systems or
AEGIS-alikes. So all of you defense-people: HANDS-OFF!
--------------------------------------------------------------------------------
----------cut here--------------------------------------------------------------
#!/bin/sh
: This is a shar archive. Extract with sh, not csh.
: This archive ends with exit, so do not worry about trailing junk.
: --------------------------- cut here --------------------------
PATH=/bin:/usr/bin:/usr/ucb
if test -f 'fft'
then echo Removing 'fft'
rm 'fft'
fi
if test -d 'fft'
then :
else echo Creating 'fft'
mkdir 'fft'
fi
chmod 'u=rwx,g=rx,o=rx' 'fft'
echo Extracting 'fft/README'
sed 's/^X//' > 'fft/README' << '+ END-OF-FILE ''fft/README'
XThis directory contains the packages realfft(3) and fft(3) that do Cooley-Tukey
Xfast Fourier transform on an arbitrary number of real or complex samples,
Xrespectively.
X
XThe package was tested on 4.[12] bsd & Sun Release 3.5, but it should work on
Xany UNIX system featuring sin(3M) and malloc(3).
X
XContents:
X complex.h - include file (for users of fft(3) routines).
X fft - manual and source of fft(3).
X realfft - manual and source of realfft(3).
X example - example of how to use realfft(3).
X
XSee the README's in fft/, realfft/ and example/ to find out how to compile and
Xuse things.
+ END-OF-FILE fft/README
chmod 'u=rw,g=r,o=r' 'fft/README'
set `wc -c 'fft/README'`
count=$1
case $count in
585) :;;
*) echo 'Bad character count in ''fft/README' >&2
echo 'Count should be 585' >&2
esac
if test -f 'fft/complex'
then echo Removing 'fft/complex'
rm 'fft/complex'
fi
if test -d 'fft/complex'
then :
else echo Creating 'fft/complex'
mkdir 'fft/complex'
fi
chmod 'u=rwx,g=rx,o=rx' 'fft/complex'
echo Extracting 'fft/complex/Makefile'
sed 's/^X//' > 'fft/complex/Makefile' << '+ END-OF-FILE ''fft/complex/Makefile'
XLTYPE=A18
XTARGET=fft.a
XCFLAGS=-O
XPREFLAGS=
XIDIRS=-I../
XLLIBS=
X
XLINT=lint -uhbaxpc
XCTAGS=ctags
XPRINT=@pr -t
X
XCFILES=\
X fourier.c\
X ft.c\
X w.c
XHFILES=\
X /usr/include/math.h\
X ../complex.h\
X w.h
XOBJECTS=\
X fourier.o\
X ft.o\
X w.o
X
X.SUFFIXES: .i
X
X$(TARGET): $(OBJECTS)
X ar rv $@ $?
X ranlib $@
X
Xlint:
X $(LINT) $(PREFLAGS) $(IDIRS) $(CFILES) $(LLIBS) -lc
X
Xtags: $(HFILES) $(CFILES)
X $(CTAGS) $(HFILES) $(CFILES)
X
Xprint: fft.man
X $(PRINT) fft.man tech complex.h w.h $(CFILES)
X
Xfft.man: fft.3
X @nroff -man fft.3 > fft.man
X
X.c.o:
X $(CC) $(CFLAGS) -c $(IDIRS) $<
X
X.c.i:
X $(CC) $(CFLAGS) -P $(IDIRS) $<
X
X.c.s:
X $(CC) $(CFLAGS) -S $(IDIRS) $<
X
Xfourier.o:\
X ../complex.h\
X w.h
X
Xft.o:\
X ../complex.h\
X w.h
X
Xw.o:\
X ../complex.h\
X w.h\
X /usr/include/math.h\
+ END-OF-FILE fft/complex/Makefile
chmod 'u=rw,g=r,o=r' 'fft/complex/Makefile'
set `wc -c 'fft/complex/Makefile'`
count=$1
case $count in
741) :;;
*) echo 'Bad character count in ''fft/complex/Makefile' >&2
echo 'Count should be 741' >&2
esac
echo Extracting 'fft/complex/README'
sed 's/^X//' > 'fft/complex/README' << '+ END-OF-FILE ''fft/complex/README'
XThis fft(3) package does Cooley-Tukey fast Fourier transform on an arbitrary
Xnumber of complex samples.
X
XHow to make the stuff:
X make - create library "fft.a"
X make fft.man - nroff manual page "fft.3" into "fft.man"
X make print - print source, "tech" and "fft.man"
X
XFile "tech" contains a short description of the functions and variables in the
Ximplementation.
X
XPrograms using fft(3), should include "../complex.h" and be loaded (ld(1)) with
Xlibraries "fft.a" and "/usr/lib/libm.a". The package also uses malloc(3) of
Xthe standard C-library.
+ END-OF-FILE fft/complex/README
chmod 'u=rw,g=r,o=r' 'fft/complex/README'
set `wc -c 'fft/complex/README'`
count=$1
case $count in
544) :;;
*) echo 'Bad character count in ''fft/complex/README' >&2
echo 'Count should be 544' >&2
esac
echo Extracting 'fft/complex/fft.3'
sed 's/^X//' > 'fft/complex/fft.3' << '+ END-OF-FILE ''fft/complex/fft.3'
X.TH FFT 3
X.SH NAME
Xfft, rft \- forward and reverse complex fourier transform
X.SH SYNOPSIS
X.nf
X#include "complex.h"
X
Xfft (in, n, out)
XCOMPLEX *in;
Xunsigned n;
XCOMPLEX *out;
X
Xrft (in, n, out)
XCOMPLEX *in;
Xunsigned n;
XCOMPLEX *out;
X
Xc_re (c)
XCOMPLEX c;
X
Xc_im (c)
XCOMPLEX c;
X.fi
X.SH DESCRIPTION
X.I
XFft
Xand
X.I rft
Xperform, respectively, forward and reverse discrete
Xfast fourier transform on the
X.I n
X(an arbitrary number) complex
Xsamples of array
X.IR in .
XThe result is placed in
X.IR out .
X.br
XThe functions are a recursive implementation of the Cooley-Tukey algorithm.
XBoth use O
X.RI ( n
X*
X.RI ( p1
X+ .. +
X.IR pk ))
Xoperations, where
X.I pi
Xis the
X.IR i -th
Xfactor in the
Xprime-decomposition of size
X.I k
Xof
X.IR n .
X.br
XBoth functions compute a sine/cosine table internally.
XThis table is not recomputed on successive calls with the same
X.IR n .
X
X.I C_re
Xand
X.I c_im
Xare C preprocessor defines that yield the real and imaginary
Xparts of
X.IR c ,
Xrespectively.
XThey are used to assign and inspect arrays
X.I in
Xand
X.IR out .
X.SH FILES
Xfft.a \- library containing fft and rft.
X.br
X/usr/lib/libm.a \- library used by fft.a.
X.SH DIAGNOSTICS
X.I Fft
Xand
X.I rft
Xreturn -1 if allocating space for the internal table fails, else 0.
X.SH BUGS
XThe original contents of
X.I in
Xare destroyed.
X
XThe transform isn't optimized for interesting special cases of
X.IR n ,
Xe.g.\&
X.I n
Xis a power of 2, although it will work in O
X.RI ( n
X* 2log
X.RI ( n )).
X
XThe error for a forward-reverse sequence is about 10e-13 for
X.I n
X= 1024.
X.SH AUTHOR
XPeter Valkenburg (valke@cs.vu.nl).
+ END-OF-FILE fft/complex/fft.3
chmod 'u=rw,g=r,o=r' 'fft/complex/fft.3'
set `wc -c 'fft/complex/fft.3'`
count=$1
case $count in
1548) :;;
*) echo 'Bad character count in ''fft/complex/fft.3' >&2
echo 'Count should be 1548' >&2
esac
echo Extracting 'fft/complex/fourier.c'
sed 's/^X//' > 'fft/complex/fourier.c' << '+ END-OF-FILE ''fft/complex/fourier.c'
X/*
X * "fourier.c", Pjotr '87.
X */
X
X#include <complex.h>
X#include "w.h"
X
X/*
X * Recursive (reverse) complex fast Fourier transform on the n
X * complex samples of array in, with the Cooley-Tukey method.
X * The result is placed in out. The number of samples, n, is arbitrary.
X * The algorithm costs O (n * (r1 + .. + rk)), where k is the number
X * of factors in the prime-decomposition of n (also the maximum
X * depth of the recursion), and ri is the i-th primefactor.
X */
XFourier (in, n, out)
XCOMPLEX *in;
Xunsigned n;
XCOMPLEX *out;
X{
X unsigned r;
X unsigned radix ();
X
X if ((r = radix (n)) < n)
X split (in, r, n / r, out);
X join (in, n / r, n, out);
X}
X
X/*
X * Give smallest possible radix for n samples.
X * Determines (in a rude way) the smallest primefactor of n.
X */
Xstatic unsigned radix (n)
Xunsigned n;
X{
X unsigned r;
X
X if (n < 2)
X return 1;
X
X for (r = 2; r < n; r++)
X if (n % r == 0)
X break;
X return r;
X}
X
X/*
X * Split array in of r * m samples in r parts of each m samples,
X * such that in [i] goes to out [(i % r) * m + (i / r)].
X * Then call for each part of out Fourier, so the r recursively
X * transformed parts will go back to in.
X */
Xstatic split (in, r, m, out)
XCOMPLEX *in;
Xregister unsigned r, m;
XCOMPLEX *out;
X{
X register unsigned k, s, i, j;
X
X for (k = 0, j = 0; k < r; k++)
X for (s = 0, i = k; s < m; s++, i += r, j++)
X out [j] = in [i];
X
X for (k = 0; k < r; k++, out += m, in += m)
X Fourier (out, m, in);
X}
X
X/*
X * Sum the n / m parts of each m samples of in to n samples in out.
X * r - 1
X * Out [j] becomes sum in [j % m] * W (j * k). Here in is the k-th
X * k = 0 k n k
X * part of in (indices k * m ... (k + 1) * m - 1), and r is the radix.
X * For k = 0, a complex multiplication with W (0) is avoided.
X */
Xstatic join (in, m, n, out)
XCOMPLEX *in;
Xregister unsigned m, n;
XCOMPLEX *out;
X{
X register unsigned i, j, jk, s;
X
X for (s = 0; s < m; s++)
X for (j = s; j < n; j += m) {
X out [j] = in [s];
X for (i = s + m, jk = j; i < n; i += m, jk += j)
X c_add_mul (out [j], in [i], W (n, jk));
X }
X}
+ END-OF-FILE fft/complex/fourier.c
chmod 'u=rw,g=r,o=r' 'fft/complex/fourier.c'
set `wc -c 'fft/complex/fourier.c'`
count=$1
case $count in
2046) :;;
*) echo 'Bad character count in ''fft/complex/fourier.c' >&2
echo 'Count should be 2046' >&2
esac
echo Extracting 'fft/complex/ft.c'
sed 's/^X//' > 'fft/complex/ft.c' << '+ END-OF-FILE ''fft/complex/ft.c'
X/*
X * "ft.c", Pjotr '87.
X */
X
X#include <complex.h>
X#include "w.h"
X
X/*
X * Forward Fast Fourier Transform on the n samples of complex array in.
X * The result is placed in out. The number of samples, n, is arbitrary.
X * The W-factors are calculated in advance.
X */
Xint fft (in, n, out)
XCOMPLEX *in;
Xunsigned n;
XCOMPLEX *out;
X{
X unsigned i;
X
X for (i = 0; i < n; i++)
X c_conj (in [i]);
X
X if (W_init (n) == -1)
X return -1;
X
X Fourier (in, n, out);
X
X for (i = 0; i < n; i++) {
X c_conj (out [i]);
X c_realdiv (out [i], n);
X }
X
X return 0;
X}
X
X/*
X * Reverse Fast Fourier Transform on the n complex samples of array in.
X * The result is placed in out. The number of samples, n, is arbitrary.
X * The W-factors are calculated in advance.
X */
Xrft (in, n, out)
XCOMPLEX *in;
Xunsigned n;
XCOMPLEX *out;
X{
X if (W_init (n) == -1)
X return -1;
X
X Fourier (in, n, out);
X
X return 0;
X}
+ END-OF-FILE fft/complex/ft.c
chmod 'u=rw,g=r,o=r' 'fft/complex/ft.c'
set `wc -c 'fft/complex/ft.c'`
count=$1
case $count in
865) :;;
*) echo 'Bad character count in ''fft/complex/ft.c' >&2
echo 'Count should be 865' >&2
esac
echo Extracting 'fft/complex/tech'
sed 's/^X//' > 'fft/complex/tech' << '+ END-OF-FILE ''fft/complex/tech'
X Short technical description of functions in the fft(3) package
X
X
X"ft.c"
XThe entry-points:
X fft - Forward Complex Fast Fourier Transform
X rft - Reverse Complex Fast Fourier Transform
X
X"fourier.c"
XRecursive implementation of the Cooley-Tukey algorithm:
X Fourier - top level call
X radix - determine radix for a number of samples
X split - split samples in groups, and recursively call Fourier
X join - join (add) groups of samples into a new group
X
X"complex.h"
XManipulation of complex numbers:
X COMPLEX - type for complex numbers
X c_re - real part of complex number
X c_im - imaginary part of complex number
X c_add_mul - multiply and add complex numbers
X c_conj - convert a complex number into its conjugate
X c_realdiv - divide a complex by a real number
X
X"w.h"
XW-factors:
X W - give previously calculated W-factor
X
X"w.c"
XComputation of W-factors:
X W_factors - array of W-factors
X Nfactors - number of factors in W_factors
X W_init - prepare W-factors array
+ END-OF-FILE fft/complex/tech
chmod 'u=rw,g=r,o=r' 'fft/complex/tech'
set `wc -c 'fft/complex/tech'`
count=$1
case $count in
951) :;;
*) echo 'Bad character count in ''fft/complex/tech' >&2
echo 'Count should be 951' >&2
esac
echo Extracting 'fft/complex/w.c'
sed 's/^X//' > 'fft/complex/w.c' << '+ END-OF-FILE ''fft/complex/w.c'
X/*
X * "w.c", Pjotr '87.
X */
X
X#include <complex.h>
X#include "w.h"
X#include <math.h>
X
XCOMPLEX *W_factors = 0; /* array of W-factors */
Xunsigned Nfactors = 0; /* number of entries in W-factors */
X
X/*
X * W_init puts Wn ^ k (= e ^ (2pi * i * k / n)) in W_factors [k], 0 <= k < n.
X * If n is equal to Nfactors then nothing is done, so the same W_factors
X * array can used for several transforms of the same number of samples.
X * Notice the explicit calculation of sines and cosines, an iterative approach
X * introduces substantial errors.
X */
Xint W_init (n)
Xunsigned n;
X{
X char *malloc ();
X# define pi 3.1415926535897932384626434
X unsigned k;
X
X if (n == Nfactors)
X return 0;
X if (Nfactors != 0 && W_factors != 0)
X free ((char *) W_factors);
X if ((Nfactors = n) == 0)
X return 0;
X if ((W_factors = (COMPLEX *) malloc (n * sizeof (COMPLEX))) == 0)
X return -1;
X
X for (k = 0; k < n; k++) {
X c_re (W_factors [k]) = cos (2 * pi * k / n);
X c_im (W_factors [k]) = sin (2 * pi * k / n);
X }
X
X return 0;
X}
+ END-OF-FILE fft/complex/w.c
chmod 'u=rw,g=r,o=r' 'fft/complex/w.c'
set `wc -c 'fft/complex/w.c'`
count=$1
case $count in
996) :;;
*) echo 'Bad character count in ''fft/complex/w.c' >&2
echo 'Count should be 996' >&2
esac
echo Extracting 'fft/complex/w.h'
sed 's/^X//' > 'fft/complex/w.h' << '+ END-OF-FILE ''fft/complex/w.h'
X/*
X * "w.h", Pjotr '87.
X */
X
Xextern COMPLEX *W_factors;
Xextern unsigned Nfactors;
X
X/*
X * W gives the (already computed) Wn ^ k (= e ^ (2pi * i * k / n)).
X * Notice that the powerseries of Wn has period Nfactors.
X */
X#define W(n, k) (W_factors [((k) * (Nfactors / (n))) % Nfactors])
+ END-OF-FILE fft/complex/w.h
chmod 'u=rw,g=r,o=r' 'fft/complex/w.h'
set `wc -c 'fft/complex/w.h'`
count=$1
case $count in
283) :;;
*) echo 'Bad character count in ''fft/complex/w.h' >&2
echo 'Count should be 283' >&2
esac
echo Extracting 'fft/complex.h'
sed 's/^X//' > 'fft/complex.h' << '+ END-OF-FILE ''fft/complex.h'
X/*
X * "complex.h", Pjotr '87.
X */
X
Xtypedef struct {
X double re, im;
X } COMPLEX;
X#define c_re(c) ((c).re)
X#define c_im(c) ((c).im)
X
X/*
X * C_add_mul adds product of c1 and c2 to c.
X */
X#define c_add_mul(c, c1, c2) { COMPLEX C1, C2; C1 = (c1); C2 = (c2); \
X c_re (c) += C1.re * C2.re - C1.im * C2.im; \
X c_im (c) += C1.re * C2.im + C1.im * C2.re; }
X
X/*
X * C_conj substitutes c by its complex conjugate.
X */
X#define c_conj(c) { c_im (c) = -c_im (c); }
X
X/*
X * C_realdiv divides complex c by real.
X */
X#define c_realdiv(c, real) { c_re (c) /= (real); c_im (c) /= (real); }
+ END-OF-FILE fft/complex.h
chmod 'u=rw,g=r,o=r' 'fft/complex.h'
set `wc -c 'fft/complex.h'`
count=$1
case $count in
583) :;;
*) echo 'Bad character count in ''fft/complex.h' >&2
echo 'Count should be 583' >&2
esac
if test -f 'fft/example'
then echo Removing 'fft/example'
rm 'fft/example'
fi
if test -d 'fft/example'
then :
else echo Creating 'fft/example'
mkdir 'fft/example'
fi
chmod 'u=rwx,g=rx,o=rx' 'fft/example'
echo Extracting 'fft/example/README'
sed 's/^X//' > 'fft/example/README' << '+ END-OF-FILE ''fft/example/README'
XAn example of realfft(3) usage is in example.c. It contains two fourier
Xtransforms on real data.
X
XCompile with:
X cc example.c ../real/realfft.a ../complex/fft.a -lm
+ END-OF-FILE fft/example/README
chmod 'u=rw,g=r,o=r' 'fft/example/README'
set `wc -c 'fft/example/README'`
count=$1
case $count in
166) :;;
*) echo 'Bad character count in ''fft/example/README' >&2
echo 'Count should be 166' >&2
esac
echo Extracting 'fft/example/example.c'
sed 's/^X//' > 'fft/example/example.c' << '+ END-OF-FILE ''fft/example/example.c'
X/*
X * Test for realfft(3).
X */
X
X#define N 8
X#define pi 3.1415926535897932384626434
X
Xdouble sin (), cos ();
Xchar *malloc ();
X
Xmain ()
X{
X unsigned i, j;
X
X double in [N], out [N];
X
X printf ("Example #1:\n");
X for (i = 0; i < N; i++)
X in [i] = i;
X printsamp (in, N);
X
X realfft (in, N, out);
X printf ("After a fast fft\n");
X printampl (out, N);
X
X /* check */
X printf ("A reverse slow ft gives:\n");
X srft (out, N, in);
X printsamp (in, N);
X
X printf ("And the reverse fast ft yields:\n");
X realrft (out, N, in);
X printsamp (in, N);
X
X printf ("\n\nExample #2\n");
X for (i = 0; i < N; i++) {
X in [i] = 0;
X for (j = 0; j <= N / 2; j++)
X in [i] += cos (2 * pi * i * j / N) +
X sin (2 * pi * i * j / N);
X }
X printsamp (in, N);
X
X realfft (in, N, out);
X printf ("After a forward fast ft:\n");
X printampl (out, N);
X
X /* check */
X printf ("A reverse slow ft yields:\n");
X srft (out, N, in);
X printsamp (in, N);
X
X printf ("And a reverse fast ft gives:\n");
X realrft (out, N, in);
X printsamp (in, N);
X}
X
Xprintampl (ampl, n)
Xdouble *ampl;
Xunsigned n;
X{
X unsigned i;
X
X printf ("Amplitudes\n");
X
X if (n == 0)
X return;
X
X printf ("%f (dc)\n", ampl [0]);
X for (i = 1; i < (n + 1) / 2; i++)
X printf ("%f, %f (%u-th harmonic)\n", ampl [2 * i - 1],
X ampl [2 * i], i);
X if (n % 2 == 0)
X printf ("%f (Nyquist)\n", ampl [n - 1]);
X
X printf ("\n");
X}
X
Xprintsamp (samp, n)
Xdouble *samp;
Xunsigned n;
X{
X unsigned i;
X printf ("Samples\n");
X
X for (i = 0; i < n; i++)
X printf ("%f\n", samp [i]);
X
X printf ("\n");
X}
X
X/*
X * Slow reverse fourier transform. In [0] contains dc, in [n - 1] Nyquist.
X * This is just a gimmick to compare with realfft(3).
X */
Xsrft (in, n, out)
Xdouble *in;
Xunsigned n;
Xdouble *out;
X{
X unsigned i, j;
X
X for (i = 0; i < n; i++) {
X out [i] = in [0]; /* dc */
X for (j = 1; j < (n + 1) / 2; j++) /* j-th harmonic */
X out [i] += in [2 * j - 1] * cos (2 * pi * i * j / n) +
X in [2 * j] * sin (2 * pi * i * j / n);
X if (n % 2 == 0) /* Nyquist */
X out [i] += in [n - 1] * cos (2 * pi * i * j / n);
X }
X}
+ END-OF-FILE fft/example/example.c
chmod 'u=rw,g=r,o=r' 'fft/example/example.c'
set `wc -c 'fft/example/example.c'`
count=$1
case $count in
2026) :;;
*) echo 'Bad character count in ''fft/example/example.c' >&2
echo 'Count should be 2026' >&2
esac
if test -f 'fft/real'
then echo Removing 'fft/real'
rm 'fft/real'
fi
if test -d 'fft/real'
then :
else echo Creating 'fft/real'
mkdir 'fft/real'
fi
chmod 'u=rwx,g=rx,o=rx' 'fft/real'
echo Extracting 'fft/real/Makefile'
sed 's/^X//' > 'fft/real/Makefile' << '+ END-OF-FILE ''fft/real/Makefile'
XLTYPE=A18
XTARGET=realfft.a
XCFLAGS=-O
XPREFLAGS=
XIDIRS=-I../
XLLIBS=
X
XLINT=lint -uhbaxc
XCTAGS=ctags
XPRINT=@pr -t
X
XCFILES=\
X realfft.c
XHFILES=\
X ../complex.h
XOBJECTS=\
X realfft.o
X
X.SUFFIXES: .i
X
X$(TARGET): $(OBJECTS)
X ar rv $@ $?
X ranlib $@
X
Xlint:
X $(LINT) $(PREFLAGS) $(IDIRS) $(CFILES) $(LLIBS) -lc
X
Xtags: $(HFILES) $(CFILES)
X $(CTAGS) $(HFILES) $(CFILES)
X
Xprint: realfft.man
X $(PRINT) realfft.man $(CFILES)
X
Xrealfft.man: realfft.3
X @nroff -man realfft.3 > realfft.man
X
X.c.o:
X $(CC) $(CFLAGS) -c $(IDIRS) $<
X
X.c.i:
X $(CC) $(CFLAGS) -P $(IDIRS) $<
X
X.c.s:
X $(CC) $(CFLAGS) -S $(IDIRS) $<
X
Xrealfft.o:\
X ../complex.h
+ END-OF-FILE fft/real/Makefile
chmod 'u=rw,g=r,o=r' 'fft/real/Makefile'
set `wc -c 'fft/real/Makefile'`
count=$1
case $count in
611) :;;
*) echo 'Bad character count in ''fft/real/Makefile' >&2
echo 'Count should be 611' >&2
esac
echo Extracting 'fft/real/README'
sed 's/^X//' > 'fft/real/README' << '+ END-OF-FILE ''fft/real/README'
XThis realfft(3) package does Cooley-Tukey fast Fourier transform on an arbitrary
Xnumber of real samples. It uses fft(3) for the actual complex transform.
X
XHow to make the stuff:
X make - create library "realfft.a"
X make fft.man - nroff manual page "realfft.3" into "realfft.man"
X make print - print source and "realfft.man"
X
XPrograms using realfft(3), should be loaded (ld(1)) with libraries "realfft.a",
X"../complex/fft.a" and "/usr/lib/libm.a". The package also uses malloc(3) of
Xthe standard C-library.
+ END-OF-FILE fft/real/README
chmod 'u=rw,g=r,o=r' 'fft/real/README'
set `wc -c 'fft/real/README'`
count=$1
case $count in
508) :;;
*) echo 'Bad character count in ''fft/real/README' >&2
echo 'Count should be 508' >&2
esac
echo Extracting 'fft/real/realfft.3'
sed 's/^X//' > 'fft/real/realfft.3' << '+ END-OF-FILE ''fft/real/realfft.3'
X.TH REALFFT 3
X.SH NAME
Xrealfft, realrft \- forward and reverse real fourier transform
X.SH SYNOPSIS
X.nf
Xrealfft (in, n, out)
Xdouble *in;
Xunsigned n;
Xdouble *out;
X
Xrealrft (in, n, out)
Xdouble *in;
Xunsigned n;
Xdouble *out;
X.fi
X.SH DESCRIPTION
X.I Realfft
Xand
X.I realrft
Xperform, respectively, forward and reverse discrete
Xfast fourier transform on the
X.I n
X(an arbitrary number) reals
Xof array
X.IR in .
XThe result (also
X.I n
Xreals) is placed in array
X.IR out .
XThe original contents of
X.I in
Xare not disturbed.
X
XThe format of the
X.I out
Xarray of
X.I realfft
Xand the
X.I in
Xarray of
X.I realrft
Xis as follows:
XThe cosine component of the dc frequency is under index 0
Xand the
X.IR i -th
Xharmonic's cosine and sine are under, respectively, index
X2 *
X.I i
X- 1 and 2 *
X.IR i .
XIf
X.I n
Xis even then the cosine component of the
XNyquist frequency is under index
X.I n
X- 1.
XNote that the dc and Nyquist sine components need not be passed, since they
Xare always 0.
X
XThe actual transform is done by
X.IR fft (3)
Xin complex space.
X.SH "SEE ALSO"
Xfft(3)
X.SH FILES
Xrealfft.a \- contains realfft and realrft.
X.br
Xfft.a \- library used by realfft.a.
X.br
X/usr/lib/libm.a \- library used by fft.a.
X.SH DIAGNOSTICS
X.I Realfft
Xand
X.I realrft
Xreturn -1 if routines in
X.IR fft (3)
Xfail, else 0.
X.SH BUGS
XThe error for a forward-reverse sequence is about 1e-13 for
X.I n
X= 1024.
X.SH AUTHOR
XPeter Valkenburg (valke@cs.vu.nl)
+ END-OF-FILE fft/real/realfft.3
chmod 'u=rw,g=r,o=r' 'fft/real/realfft.3'
set `wc -c 'fft/real/realfft.3'`
count=$1
case $count in
1391) :;;
*) echo 'Bad character count in ''fft/real/realfft.3' >&2
echo 'Count should be 1391' >&2
esac
echo Extracting 'fft/real/realfft.c'
sed 's/^X//' > 'fft/real/realfft.c' << '+ END-OF-FILE ''fft/real/realfft.c'
X/*
X * "realfft.c", Pjotr '87
X */
X
X/*
X * Bevat funkties realfft en realrft die resp. forward en reverse fast fourier
X * transform op reele samples doen. Gebruikt pakket fft(3).
X */
X
X#include <complex.h>
X
Xchar *malloc ();
X
X/*
X * Reele forward fast fourier transform van n samples van in naar
X * amplitudes van out.
X * De cosinus komponent van de dc komt in out [0], dan volgen in
X * out [2 * i - 1] en out [2 * i] steeds resp. de cosinus en sinus
X * komponenten van de i-de harmonische. Bij een even aantal samples
X * bevat out [n - 1] de cosinus komponent van de Nyquist frequentie.
X * Extraatje: Na afloop is in onaangetast.
X */
Xrealfft (in, n, out)
Xdouble *in;
Xunsigned n;
Xdouble *out;
X{
X COMPLEX *c_in, *c_out;
X unsigned i;
X
X if (n == 0 ||
X (c_in = (COMPLEX *) malloc (n * sizeof (COMPLEX))) == 0 ||
X (c_out = (COMPLEX *) malloc (n * sizeof (COMPLEX))) == 0)
X return;
X
X for (i = 0; i < n; i++) {
X c_re (c_in [i]) = in [i];
X c_im (c_in [i]) = 0;
X }
X
X fft (c_in, n, c_out);
X
X out [0] = c_re (c_out [0]); /* cos van dc */
X for (i = 1; i < (n + 1) / 2; i++) { /* cos/sin i-de harmonische */
X out [2 * i - 1] = c_re (c_out [i]) * 2;
X out [2 * i] = c_im (c_out [i]) * -2;
X }
X if (n % 2 == 0) /* cos van Nyquist */
X out [n - 1] = c_re (c_out [n / 2]);
X
X free ((char *) c_in);
X free ((char *) c_out);
X}
X
X/*
X * Reele reverse fast fourier transform van amplitudes van in naar
X * n samples van out.
X * De cosinus komponent van de dc staat in in [0], dan volgen in
X * in [2 * i - 1] en in [2 * i] steeds resp. de cosinus en sinus
X * komponenten van de i-de harmonische. Bij een even aantal samples
X * bevat in [n - 1] de cosinus komponent van de Nyquist frequentie.
X * Extraatje: Na afloop is in onaangetast.
X */
Xrealrft (in, n, out)
Xdouble *in;
Xunsigned n;
Xdouble *out;
X{
X COMPLEX *c_in, *c_out;
X unsigned i;
X
X if (n == 0 ||
X (c_in = (COMPLEX *) malloc (n * sizeof (COMPLEX))) == 0 ||
X (c_out = (COMPLEX *) malloc (n * sizeof (COMPLEX))) == 0)
X return;
X
X c_re (c_in [0]) = in [0]; /* dc */
X c_im (c_in [0]) = 0;
X for (i = 1; i < (n + 1) / 2; i++) { /* geconj. symm. harmonischen */
X c_re (c_in [i]) = in [2 * i - 1] / 2;
X c_im (c_in [i]) = in [2 * i] / -2;
X c_re (c_in [n - i]) = in [2 * i - 1] / 2;
X c_im (c_in [n - i]) = in [2 * i] / 2;
X }
X if (n % 2 == 0) { /* Nyquist */
X c_re (c_in [n / 2]) = in [n - 1];
X c_im (c_in [n / 2]) = 0;
X }
X
X rft (c_in, n, c_out);
X
X for (i = 0; i < n; i++)
X out [i] = c_re (c_out [i]);
X
X free ((char *) c_in);
X free ((char *) c_out);
X}
+ END-OF-FILE fft/real/realfft.c
chmod 'u=rw,g=r,o=r' 'fft/real/realfft.c'
set `wc -c 'fft/real/realfft.c'`
count=$1
case $count in
2503) :;;
*) echo 'Bad character count in ''fft/real/realfft.c' >&2
echo 'Count should be 2503' >&2
esac
exit 0
--
Please send comp.sources.unix-related mail to rsalz@uunet.uu.net.