fft.c.sh

UNIX系统之下的快速傅立叶变换包

#!/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'`