v04i009: An FFT program
Steve Sampson
sampson at killer.DALLAS.TX.US
Mon Aug 1 11:44:07 AEST 1988
Posting-number: Volume 4, Issue 9
Submitted-by: "Steve Sampson" <sampson at killer.DALLAS.TX.US>
Archive-name: dos-fft
[MS-DOS C of some kind -- you may have to patch "fopen"'s under Unix. At
least. No binaries were provided. ++bsa]
#! /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 gen.c fft.c
# Wrapped by sampson at killer on Sun Jul 31 20:40:12 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\" \(3814 characters\)
sed "s/^X//" >readme.doc <<'END_OF_readme.doc'
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. I don't do windows though...
X
XNeeds a graphic interface, but the time escapes me. Please upload any better
Xgraphic 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.
X
XFor an example try:
X
X gen 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.exe 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 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 file data can be used by other programs as needed.
X
XBy using negative frequencies in gen.exe you can generate opening targets:
X
X gen 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.
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
XGood luck, have fun with it,
XSteve Sampson, CompuServe: 75136,626 Unix: sampson at killer.dallas.tx.us
END_OF_readme.doc
if test 3814 -ne `wc -c <readme.doc`; then
echo shar: \"readme.doc\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f gen.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"gen.c\"
else
echo shar: Extracting \"gen.c\" \(1292 characters\)
sed "s/^X//" >gen.c <<'END_OF_gen.c'
X/*
X * gen.c
X *
X * C Version 1.0 by Steve Sampson, Public Domain
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 samples output
X */
X
X#include <stdio.h>
X#include <alloc.h>
X#include <math.h>
X
X#define PI2 ((double)2.0 * M_PI)
X
Xmain(argc, argv)
Xint argc;
Xchar *argv[];
X{
X FILE *fp;
X double sample, freq, time, *real, *imag;
X int loop, samples;
X
X if (argc != 3) {
X printf("Usage: gen samples output_file\n");
X printf("Where samples is a power of 2\n");
X exit(-1);
X }
X
X if ((fp = fopen(argv[2], "wb")) == (FILE *)NULL) {
X printf("Unable to create write file\n");
X exit(-1);
X }
X
X samples = abs(atoi(argv[1]));
X
X real = (double *)malloc(samples * sizeof(double));
X imag = (double *)malloc(samples * sizeof(double));
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(PI2 * freq * time);
X imag[loop] = -cos(PI2 * freq * time);
X time += sample;
X }
X
X fwrite(real, sizeof(double), samples, fp);
X fwrite(imag, sizeof(double), samples, fp);
X
X fclose(fp);
X putchar('\n');
X}
X
X/* EOF */
END_OF_gen.c
if test 1292 -ne `wc -c <gen.c`; then
echo shar: \"gen.c\" 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\" \(4203 characters\)
sed "s/^X//" >fft.c <<'END_OF_fft.c'
X/*
X * fft.c
X *
X * C Version 1.0 by Steve Sampson, Public Domain
X *
X * This program is based on the work by W. D. Stanley
X * and S. J. Peterson, Old Dominion University.
X *
X * This program produces a Frequency Domain display
X * from 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_data output_data
X * Where 'samples' is a power of two
X *
X * Array Version for Turbo C 1.5
X */
X
X/* Includes */
X
X#include <stdlib.h>
X#include <stdio.h>
X#include <math.h>
X#include <alloc.h>
X
X/* Defines */
X
X#define TWO_PI ((double)2.0 * M_PI)
X
X/* Globals */
X
Xint samples, power;
Xdouble *real, *imag, max;
XFILE *fpi, *fpo;
X
X/* Prototypes and forward declarations */
X
Xvoid fft(void), max_amp(void), display(void);
Xint permute(int);
Xdouble magnitude(int);
X
X/* The program */
X
Xmain(argc, argv)
Xint argc;
Xchar *argv[];
X{
X int n;
X
X if (argc != 4) {
Xerr1: fprintf(stderr, "Usage: fft samples input output\n");
X fprintf(stderr, "Where samples is a power of 2\n");
X exit(1);
X }
X
X samples = abs(atoi(argv[1]));
X power = log10((double)samples) / log10((double)2.0);
X
X if ((real = (double *)malloc(samples * sizeof(double))) == NULL)
Xerr2: fprintf(stderr, "Out of memory\n");
X
X if ((imag = (double *)malloc(samples * sizeof(double))) == NULL)
X goto err2;
X
X if ((fpo = fopen(argv[3], "wb")) == (FILE *)NULL)
X goto err1;
X
X if ((fpi = fopen(argv[2], "rb")) != (FILE *)NULL) {
X fread(real, sizeof(double), samples, fpi);
X fread(imag, sizeof(double), samples, fpi);
X fclose(fpi);
X }
X else
X goto err1;
X
X fft();
X max_amp();
X display();
X
X fwrite(real, sizeof(double), samples, fpo);
X fwrite(imag, sizeof(double), samples, fpo);
X
X fclose(fpo);
X}
X
X
Xvoid fft()
X{
X unsigned i1, i2, i3, i4, y;
X int loop, loop1, loop2;
X double a1, a2, b1, b2, z1, z2, v;
X
X /* Scale the data */
X
X for (loop = 0; loop < samples; loop++) {
X real[loop] /= (double)samples;
X imag[loop] /= (double)samples;
X }
X
X i1 = samples >> 1;
X i2 = 1;
X v = TWO_PI * ((double)1.0 / (double)samples);
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 / i1);
X z1 = cos(v * y);
X z2 = -sin(v * y);
X
X for (loop2 = i3; loop2 < i4; loop2++) {
X a1 = real[loop2];
X a2 = imag[loop2];
X
X b1 = z1*real[loop2+i1] - z2*imag[loop2+i1];
X b2 = z2*real[loop2+i1] + 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 * Find maximum amplitude
X */
X
Xvoid max_amp()
X{
X double mag;
X int loop;
X
X max = (double)0.0;
X for (loop = 0; loop < samples; loop++) {
X if ((mag = magnitude(loop)) > max)
X max = mag;
X }
X}
X
X/*
X * Display the frequency domain.
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
Xvoid display()
X{
X int c, n, x, loop;
X
X n = samples / 2;
X
X for (loop = n; loop < samples; loop++) {
X x = (int)(magnitude(loop) * (double)56.0 / max);
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 / max);
X printf("%d\t|", loop + n);
X c = 0;
X while (++c <= x)
X putchar('=');
X
X putchar('\n');
X }
X}
X
X/*
X * Calculate Power Magnitude
X */
X
Xdouble magnitude(n)
Xint n;
X{
X n = permute(n);
X return (sqrt(real[n] * real[n] + imag[n] * imag[n]));
X}
X
X/*
X * Bit reverse the number
X *
X * Change 11100000b to 00000111b or vice-versa
X */
X
Xint permute(index)
Xint index;
X{
X int n1, result, loop;
X
X n1 = samples;
X result = 0;
X
X for (loop = 0; loop < power; loop++) {
X n1 >>= 1; /* n1 / 2.0 */
X if (index < n1)
X continue;
X
X result += (int) pow((double)2.0, (double)loop);
X index -= n1;
X }
X
X return result;
X}
X
X/* EOF */
END_OF_fft.c
if test 4203 -ne `wc -c <fft.c`; then
echo shar: \"fft.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