v04i028: An FFT Program (Unix and MS-DOS)

Steve Sampson sampson at killer.DALLAS.TX.US.UUCP
Mon Aug 15 19:32:05 AEST 1988


Posting-number: Volume 4, Issue 28
Submitted-by: "Steve Sampson" <sampson at killer.DALLAS.TX.US.UUCP>
Archive-name: fft.port

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#		"End of shell archive."
# Contents:  readme.doc fft.c gen_sine.c gen_pulse.c
# Wrapped by sampson at killer on Mon Aug 15 04:27:11 1988
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f readme.doc -a "${1}" != "-c" ; then 
  echo shar: Will not over-write existing file \"readme.doc\"
else
echo shar: Extracting \"readme.doc\" \(5214 characters\)
sed "s/^X//" >readme.doc <<'END_OF_readme.doc'
XReadme.doc  15 August 1988
X
X+----------------------------------------------------------------------------+
X| Thanks for the feedback on the original posting, here's an archive with all|
X| the known bugs fixed, and some improvements. [srs]			     |
X|									     |
X| This has compiled and run on Turbo-C V1.5, MS-C V5.0 and Unix SV3.1	     |
X|									     |
X| If your wondering what all the castings about, I had a hell of a time :-)  |
X| It works I ain't touching it (smile).  But yes, it's probably overkill.    |
X| The last killer bug was that Unix seems to truncate on an (int) cast, and  |
X| Turbo-C and MS-C seem to round.  At least that's what I decided, so there. |
X+----------------------------------------------------------------------------+
X
XI originally saw an FFT program in Byte Magazine many years ago.  I wrote
Xa version for BASIC that worked pretty good.  Then I thought I'd translate
Xit into C.  These programs are the result.
X
XNeeds a graphic interface.  I'm interested in any better graphic versions.
X
XThe original Byte Magazine program was designed for real data only.  In my
Xexperiments I needed to preserve both real and imaginary data.  If you feed
Xthe FFT real data only, then the output will be a mirror image, and you can
Xignore the left side.  Two signal generators are included.  One generates
Xsine waves (gen_sine) and the other generates pulses (gen_pulse).
X
XFor an example try:
X
X	gen_sine 16 in
X	1000
X	3000
X
XWhich will sample the 1 Khz data every 333 microseconds (1 / 3 Khz).
XNote: The sample frequency should be greater than 2 times the input
Xfrequency (Nyquist and all that...).
X
XThen run fft like so:
X
X	fft 16 in out
X
XAnd you should see a display like so:
X
X0	|=======				      (-1500.0 Hz)
X1	|=====					      (-1312.5 Hz)
X2	|====					      (-1125.0 Hz)
X3	|====					       (-937.0 Hz)
X4	|===					       (-750.0 Hz)
X5	|===					       (-562.5 Hz)
X6	|===					       (-375.0 Hz)
X7	|===					       (-187.5 Hz)
X8	|====		<-------   DC			(000.0 Hz)
X9	|====		<-------   Fundamental		(187.5 Hz)
X10	|======		<-------   Second Harmonic	(375.0 Hz)
X11	|========					(562.5 Hz)
X12	|==============					(750.0 Hz)
X13	|========================================================
X14	|============================		    (1125.0 Hz)	^
X15	|===========				    (1312.5 Hz)	|
X								|
X				[13 - 8 (center)] * 187.5 = 937.0 Hz
X
XThe fundamental display frequency is:
X
X	T  = Time Increment Between Samples
X	N  = Number Of Samples
X	Tp = N * T
X
X	Then F = 1 / Tp
X
X	In the example above, the time increment between samples is
X	1 / 3000 or 333 microseconds.  N = 16, so Tp = 5333 microseconds
X	and 1 / .005333 is 187.5 Hz.
X
X	Therefore each filter is a multiple of 187.5 Hertz.  Filter 8 in this
X	example is center, so that would be zero, 9 would be one, etc.
X
XIn this case the sample interval didn't work so good for the frequency and
Xshows the limitation of the Discrete Fourier Transform in representing a
Xcontinuous signal.  A better sample rate for 1000 Hz would be 4000 Hz,
Xin which case T = 250 ms, Tp = 4 ms, and F = 250 Hz.  1000 / 250 = 4.  The
Xpower should all be in filter 12 (8 + 4) in this case.
X
XLet's run it and see:
X
X	gen_sine 16 in
X	1000
X	4000
X
X	fft 16 in out
X
X0	|
X1	|
X2	|
X3	|
X4	|
X5	|
X6	|
X7	|
X8	|
X9	|
X10	|
X11	|
X12	|========================================================
X13	|
X14	|
X15	|
X
XWell what do you know...
X
XThe output data file data can be used by other programs as needed.
X
XBy using negative frequencies in gen_sine you can generate opening targets:
X
X	gen_sine 16 in
X	-1000
X	3000
X	fft 16 in out
X
XProduces:
X
X0	|=======
X1	|===========
X2	|============================
X3	|=======================================================
X4	|==============
X5	|========
X6	|======
X7	|====
X8	|====		<-------- Zero Hertz (DC)
X9	|===
X10	|===
X11	|===
X12	|===
X13	|====
X14	|====
X15	|=====
X
XYou can see in these examples where weighting functions would be nice.
X
XFor an example of what happens when the imaginary data is not input
X(ie, zeros put in) for a 1000 Hz frequency at 3000 Hz sample rate:
X
X0	|===============
X1	|==================
X2	|===================================
X3	|========================================================
X4	|===========
X5	|====
X6	|==
X7	|=				Trash this part
X---------------------------------------------------------------------
X8	|
X9	|=
X10	|==
X11	|====
X12	|===========
X13	|=======================================================
X14	|===================================
X15	|==================
X
XThe left side is redundant and can be deleted.  This is what the original
XByte Magazine article did (December 1978 Issue).
X
XFor generating pulses, a second program gen_pulse is provided.  It pre-loads
Ximaginary data with zeros.   For example:
X
X	gen_pulse 16 in
X	.000006
X	.0000008
X
XIs a radar with a 6 microsecond pulse and 800 nanosecond range gates.
X
X	fft 16 in out
X
XWill produce:
X
X0	|
X1	|=======
X2	|
X3	|========
X4	|
X5	|============
X6	|
X7	|===================================
X8	|========================================================
X9	|===================================
X10	|
X11	|============
X12	|
X13	|========
X14	|
X15	|=======
X
XThe more filters you use, the prettier it looks.  Good luck, have fun with it.
XSteve Sampson, CompuServe: 75136,626  Unix: sampson at killer.dallas.tx.us
END_OF_readme.doc
if test 5214 -ne `wc -c <readme.doc`; then
    echo shar: \"readme.doc\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f fft.c -a "${1}" != "-c" ; then 
  echo shar: Will not over-write existing file \"fft.c\"
else
echo shar: Extracting \"fft.c\" \(7978 characters\)
sed "s/^X//" >fft.c <<'END_OF_fft.c'
X/*
X *	fft.c
X *
X *	C Version 1.6 by Steve Sampson, Public Domain, 15 August 1988
X *
X *	This program is based on the work by W. D. Stanley and S. J. Peterson,
X *	Old Dominion University.  It produces a Frequency Domain display from
X *	the Time Domain data input using the Fast Fourier Transform.
X *
X *	The Real data is generated by the in-phase (I) channel and the
X *	Imaginary data is produced by the quadrature-phase (Q) channel of
X *	a Doppler Radar receiver.  The middle filter is zero Hz.  Closing
X *	targets are displayed to the right, and Opening targets to the left.
X *
X *	Note: With Imaginary data set to zero the output is a mirror image.
X *
X *	Usage:	fft samples input output
X *		1. Samples is a power of two.
X *		2. Input is (samples * sizeof(double)) characters.
X *		3. Output is (samples * sizeof(double)) characters.
X *		4. Standard error is help or debugging messages.
X *
X *	Define SPEED to produce a program that pre-computes the sin/cos tables,
X *		this takes more memory however.
X *
X *	Define SIMULATE to substitute a quicky magnitude routine, Has a max
X *		12% error.
X *
X *	Define SYSV to create a Unix(tm) System V operating system version,
X *		else it defaults to MS-DOS(tm) and Turbo-C(tm) V1.5
X *		and MS-C V5.0.
X *
X *	See also: readme.doc, gen_pulse.c, and gen_sine.c.
X */
X
X/* Includes */
X
X#include <stdio.h>
X#include <malloc.h>	/* Rename Turbo-C alloc.h to malloc.h */
X#include <math.h>
X
X
X/* Defines */
X
X#define	TWO_PI	((double)2.0 * 3.14159265358979324)
X#define	Chunk	((unsigned)(samples * sizeof(double)))
X
X#ifdef	SYSV
X#define	RMODE	"r"
X#define	WMODE	"w"
X#else
X#define	RMODE	"rb"
X#define	WMODE	"wb"
X#endif
X
X
X/* Globals, Forward declarations */
X
Xstatic int	samples, power, permute();
Xstatic double	*real, *imag, maxn, magnitude();
Xstatic void	fft(), max_amp(), display(), err_report();
X
X#ifdef SPEED
Xstatic double	*sine, *cosine;
Xstatic void	build_trig();
X#endif
X
Xstatic FILE	*fpi, *fpo;
X
X
X/* The program */
X
Xmain(argc, argv)
Xint	argc;
Xchar	*argv[];
X{
X	if (argc != 4)
X		err_report(0);
X
X	samples = abs(atoi(argv[1]));
X	power   = (int)(log10((double)samples) / log10((double)2.0));
X
X	/* Allocate memory for the variable arrays */
X
X	if ((real = (double *)malloc(Chunk)) == (double *)NULL)
X		err_report(1);
X
X	if ((imag = (double *)malloc(Chunk)) == (double *)NULL)
X		err_report(2);
X
X#ifdef SPEED
X	if ((sine = (double *)malloc(Chunk)) == (double *)NULL)
X		err_report(3);
X
X	if ((cosine = (double *)malloc(Chunk)) == (double *)NULL)
X		err_report(4);
X#endif
X	/* open the data files */
X
X	if ((fpi = fopen(argv[2], RMODE)) == (FILE *)NULL)
X		err_report(5);
X
X	if ((fpo = fopen(argv[3], WMODE)) == (FILE *)NULL)
X		err_report(6);
X
X	/* read in the data from the input file */
X
X	fread((char *)real, sizeof(double), (unsigned)samples, fpi);
X	fread((char *)imag, sizeof(double), (unsigned)samples, fpi);
X	fclose(fpi);
X
X#ifdef SPEED
X	build_trig();
X#endif
X	fft();
X	display();
X
X	/* write the frequency domain data to the standard output */
X
X	fwrite((char *)real, sizeof(double), (unsigned)samples, fpo);
X	fwrite((char *)imag, sizeof(double), (unsigned)samples, fpo);
X	fclose(fpo);
X
X	free((char *)real);
X	free((char *)imag);
X
X#ifdef SPEED
X	free((char *)sine);
X	free((char *)cosine);
X#endif
X
X	exit(0);
X}
X
X
Xstatic void err_report(n)
Xint	n;
X{
X	switch (n)  {
X	case 0:
X		fprintf(stderr, "Usage: fft samples\n");
X		fprintf(stderr, "Where samples is a power of two\n");
X		break;
X	case 1:
X		fprintf(stderr, "fft: Out of memory getting real space\n");
X		break;
X	case 2:
X		fprintf(stderr, "fft: Out of memory getting imag space\n");
X		free((char *)real);
X		break;
X#ifdef SPEED
X	case 3:
X		fprintf(stderr, "fft: Out of memory getting sine space\n");
X		free((char *)real);
X		free((char *)imag);
X		break;
X	case 4:
X		fprintf(stderr, "fft: Out of memory getting cosine space\n");
X		free((char *)real);
X		free((char *)imag);
X		free((char *)sine);
X		break;
X#endif
X	case 5:
X		fprintf(stderr,"fft: Unable to open data input file\n");
X		free((char *)real);
X		free((char *)imag);
X#ifdef SPEED
X		free((char *)sine);
X		free((char *)cosine);
X#endif
X		break;
X	case 6:
X		fprintf(stderr,"fft: Unable to open data output file\n");
X		fclose(fpi);
X		free((char *)real);
X		free((char *)imag);
X#ifdef SPEED
X		free((char *)sine);
X		free((char *)cosine);
X#endif
X	}
X
X	exit(1);
X}
X
X
Xstatic void fft()
X{
X	register int	loop, loop1, loop2;
X	unsigned	i1;			/* going to right shift this */
X	int		i2, i3, i4, y;
X	double		a1, a2, b1, b2, z1, z2;
X#ifndef SPEED
X	double		v;
X
X	v  = TWO_PI * ((double)1.0 / (double)samples);
X#endif
X	i1 = samples / 2;
X	i2 = 1;
X
X	/* perform the butterflies */
X
X	for (loop = 0; loop < power; loop++)  {
X		i3 = 0;
X		i4 = i1;
X
X		for (loop1 = 0; loop1 < i2; loop1++)  {
X			y = permute(i3 / (int)i1);
X#ifdef SPEED
X			z1 =  cosine[y];
X			z2 = -sine[y];
X#else
X			z1 =  cos(v * (double)y);
X			z2 = -sin(v * (double)y);
X#endif
X			for (loop2 = i3; loop2 < i4; loop2++)  {
X
X				/* Don't depend on complex precedences */
X
X				a1 = real[loop2];
X				a2 = imag[loop2];
X
X				b1  = z1 * real[loop2 + i1];
X				b1 -= z2 * imag[loop2 + i1];
X
X				b2  = z2 * real[loop2 + i1];
X				b2 += z1 * imag[loop2 + i1];
X
X				real[loop2] = a1 + b1;
X				imag[loop2] = a2 + b2;
X
X				real[loop2 + i1] = a1 - b1;
X				imag[loop2 + i1] = a2 - b2;
X			}
X
X			i3 += (i1 << 1);
X			i4 += (i1 << 1);
X		}
X
X		i1 >>= 1;
X		i2 <<= 1;
X	}
X}
X
X/*
X *	Display the frequency domain.
X *
X *	The filters are aranged so that DC is in the middle filter.
X *	Thus -Doppler is on the left, +Doppler on the right.
X */
X
Xstatic void display()
X{
X	register int	c, loop;
X	int		n, x;
X
X	max_amp();
X	n = samples / 2;
X
X	for (loop = n; loop < samples; loop++)  {
X		x = (int)(magnitude(loop) * (double)56.0 / maxn);
X		printf("%d\t|", loop - n);
X		c = 0;
X		while (++c <= x)
X			putchar('=');
X
X		putchar('\n');
X	}
X
X	for (loop = 0; loop < n; loop++)  {
X		x = (int)(magnitude(loop) * (double)56.0 / maxn);
X		printf("%d\t|", loop + n);
X		c = 0;
X		while (++c <= x)
X			putchar('=');
X
X		putchar('\n');
X	}
X}
X
X/*
X *	Find maximum amplitude
X */
X
Xstatic void max_amp()
X{
X	register int	loop;
X	double		mag;
X
X	maxn = (double)0.0;
X	for (loop = 0; loop < samples; loop++)  {
X		if ((mag = magnitude(loop)) > maxn)
X			maxn = mag;
X	}
X}
X
X/*
X *	Calculate Power Magnitude
X */
X
X#ifdef SIMULATE
X/*
X *	Simulate sqrt(i^2 + q^2).  Do this by adding half the smaller to
X *	the larger.  See also excellent description by Bob Leedom in the
X *	June 1979 Byte Magazine Technical Forum P. 188.  I originally saw
X *	The algorithm used in both an FAA Moving Target Detector (MTD) and
X *	a Westinghouse radar Envelope Detector.  Both in Hardware of course,
X *	where the binary shifting worked very fast.  Don't need high accuracy
X *	to tell you the storm is gonna remove you from Kansas (eh - Toto)...
X */
X
Xstatic double magnitude(n)
Xint	n;
X{
X	double	i, q;
X
X	n = permute(n);
X
X	i = fabs(real[n]);	/* mask the sign bit */
X	q = fabs(imag[n]);
X
X	if (i > q)
X		return (i + q / (double)2.0);	/* shift right 1 and add */
X	else
X		return (q + i / (double)2.0);
X}
X#else
X
X/*
X *	Do it the normal way
X */
X
Xstatic double magnitude(n)
Xint	n;
X{
X	double	i, q;
X
X	n = permute(n);
X	i = real[n] * real[n];
X	q = imag[n] * imag[n];
X
X	return sqrt(i + q);
X}
X#endif
X
X/*
X *	Bit reverse the number
X *
X *	Change 11100000b to 00000111b or vice-versa
X */
X
Xstatic int permute(index)
Xint	index;
X{
X	register int	loop;
X	unsigned	n1;
X	int		result;
X
X	n1 = samples;
X	result = 0;
X
X	for (loop = 0; loop < power; loop++)  {
X		n1 >>= 1;
X		if (index < n1)
X			continue;
X
X		/* Unix has a rounding error that MS-DOS compilers don't */
X		/* so go ahead and add some bits (.05) */
X
X		result += (int)((double).05 + pow((double)2.0, (double)loop));
X		index -= n1;
X	}
X
X	return result;
X}
X
X#ifdef SPEED
X/*
X *	Pre-compute the sine and cosine tables
X */
X
Xstatic void build_trig()
X{
X	register int	loop;
X	double		angle, increment;
X
X	angle     = (double)0.0;
X	increment = TWO_PI * ((double)1.0 / (double)samples);
X
X	for (loop = 0; loop < samples; loop++)  {
X		sine[loop]   = sin(angle);
X		cosine[loop] = cos(angle);
X		angle       += increment;
X	}
X}
X#endif
X
X/* EOF */
END_OF_fft.c
if test 7978 -ne `wc -c <fft.c`; then
    echo shar: \"fft.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f gen_sine.c -a "${1}" != "-c" ; then 
  echo shar: Will not over-write existing file \"gen_sine.c\"
else
echo shar: Extracting \"gen_sine.c\" \(1964 characters\)
sed "s/^X//" >gen_sine.c <<'END_OF_gen_sine.c'
X/*
X *	gen_sine.c
X *
X *	C Version 1.6 by Steve Sampson, Public Domain, 15 August 1988
X *
X *	This program is used to generate time domain sinewave data
X *	for fft.c.  If you want an opening target - negate the test frequency
X *
X *	Usage: gen_sine samples output
X */
X
X#include <stdio.h>
X#include <malloc.h>
X#include <math.h>
X
X#define	TWO_PI	((double)2.0 * 3.14159265358979324)
X#define Chunk	((unsigned)(samples * sizeof(double)))
X
X#ifdef	SYSV
X#define	WMODE	"w"
X#else
X#define	WMODE	"wb"
X#endif
X
Xstatic double	sample, freq, time, *real, *imag;
Xstatic int	loop, samples;
Xstatic void	err_report();
Xstatic FILE	*fp;
X
X
Xmain(argc, argv)
Xint	argc;
Xchar	*argv[];
X{
X	if (argc != 3)
X		err_report(0);
X
X	samples = abs(atoi(argv[1]));
X
X	if ((real = (double *)malloc(Chunk)) == (double *)NULL)
X		err_report(1);
X
X	if ((imag = (double *)malloc(Chunk)) == (double *)NULL)
X		err_report(2);
X
X	printf("Input The Test Frequency (Hz) ? ");
X	scanf("%lf", &freq);
X	printf("Input The Sampling Frequency (Hz) ? ");
X	scanf("%lf", &sample);
X	sample = (double)1.0 / sample;
X
X	time = (double)0.0;
X	for (loop = 0; loop < samples; loop++)  {
X		real[loop] =  sin(TWO_PI * freq * time);
X		imag[loop] = -cos(TWO_PI * freq * time);
X		time += sample;
X	}
X
X	if ((fp = fopen(argv[2], WMODE)) == (FILE *)NULL)
X		err_report(3);
X
X	fwrite((char *)real, sizeof(double), (unsigned)samples, fp);
X	fwrite((char *)imag, sizeof(double), (unsigned)samples, fp);
X	fclose(fp);
X
X	putchar('\n');
X
X	free((char *)real);
X	free((char *)imag);
X
X	exit(0);
X}
X
X
Xstatic void err_report(n)
Xint	n;
X{
X	switch (n)  {
X	case 0:
X		fprintf(stderr,"Usage: gen_sine samples output_file\n");
X		fprintf(stderr,"Where samples is a power of two\n");
X		break;
X	case 1:
X		fprintf(stderr,"gen_sine: Out of memory getting real space\n");
X		break;
X	case 2:
X		fprintf(stderr,"gen_sine: Out of memory getting imag space\n");
X		free((char *)real);
X		break;
X	case 3:
X		fprintf(stderr,"gen_sine: Unable to create write file\n");
X	}
X
X	exit(1);
X}
X
X/* EOF */
END_OF_gen_sine.c
if test 1964 -ne `wc -c <gen_sine.c`; then
    echo shar: \"gen_sine.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f gen_pulse.c -a "${1}" != "-c" ; then 
  echo shar: Will not over-write existing file \"gen_pulse.c\"
else
echo shar: Extracting \"gen_pulse.c\" \(1827 characters\)
sed "s/^X//" >gen_pulse.c <<'END_OF_gen_pulse.c'
X/*
X *	gen_pulse.c
X *
X *	C Version 1.6 by Steve Sampson, Public Domain, 15 August 1988
X *
X *	This program is used to generate time domain pulse data	for fft.c.
X *
X *	Usage: gen_pulse samples output
X */
X
X#include <stdio.h>
X#include <malloc.h>
X
X#define	Chunk	((unsigned)(samples * sizeof(double)))
X
X#ifdef	SYSV
X#define	WMODE	"w"
X#else
X#define	WMODE	"wb"
X#endif
X
Xstatic double	sample, pw, time, *real, *imag;
Xstatic int	loop, samples;
Xstatic void	err_report();
Xstatic FILE	*fp;
X
X
Xmain(argc, argv)
Xint	argc;
Xchar	*argv[];
X{
X	if (argc != 3)
X		err_report(0);
X
X	samples = abs(atoi(argv[1]));
X
X	if ((real = (double *)malloc(Chunk)) == (double *)NULL)
X		err_report(1);
X
X	if ((imag = (double *)malloc(Chunk)) == (double *)NULL)
X		err_report(2);
X
X	printf("Input The Pulse Width (Seconds) ? ");
X	scanf("%lf", &pw);
X	printf("Input The Sampling Time (Seconds) ? ");
X	scanf("%lf", &sample);
X
X	time = (double)0.0;
X	for (loop = 0; loop < samples; loop++)  {
X		if (time < pw)
X			real[loop] = (double)1.0;
X		else
X			real[loop] = (double)0.0;
X
X		imag[loop] = (double)0.0;
X		time += sample;
X	}
X
X	if ((fp = fopen(argv[2], WMODE)) == (FILE *)NULL)
X		err_report(3);
X
X	fwrite((char *)real, sizeof(double), (unsigned)samples, fp);
X	fwrite((char *)imag, sizeof(double), (unsigned)samples, fp);
X	fclose(fp);
X
X	putchar('\n');
X
X	free((char *)real);
X	free((char *)imag);
X
X	exit(0);
X}
X
X
Xstatic void err_report(n)
Xint	n;
X{
X	switch (n)  {
X	case 0:
X		fprintf(stderr,"Usage: gen_pulse samples output_file\n");
X		fprintf(stderr,"Where samples is a power of two\n");
X		break;
X	case 1:
X		fprintf(stderr,"gen_pulse: Out of memory getting real space\n");
X		break;
X	case 2:
X		fprintf(stderr,"gen_pulse: Out of memory getting imag space\n");
X		free((char *)real);
X		break;
X	case 3:
X		fprintf(stderr,"gen_pulse: Unable to create write file\n");
X	}
X
X	exit(1);
X}
X
X/* EOF */
END_OF_gen_pulse.c
if test 1827 -ne `wc -c <gen_pulse.c`; then
    echo shar: \"gen_pulse.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
echo shar: End of shell archive.
exit 0

^^



More information about the Comp.sources.misc mailing list