fft.c.sh

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

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