sampson@killer.DALLAS.TX.US.UUCP (Steve Sampson) (08/15/88)
Posting-number: Volume 4, Issue 28
Submitted-by: "Steve Sampson" <sampson@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@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@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
^^