Usenet 1994 October/usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso/unix/volume18/fft

Text File | 1989-03-12 | 25KB | 1,047 lines

Subject:  v18i020:  General Fast Foufier Transform package
Newsgroups: comp.sources.unix
Sender: sources
Approved: [e-mail]

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 ([e-mail]).
--------------------------------------------------------------------------------
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 ([e-mail]).
+ 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 ([e-mail])
+ 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