5 September 1988

FFT Program for Unix C compiler
Steve Sampson, Box 45668, Tinker AFB, OK 73145
sampson@killer	(TAC Trained)


   This Unix version of the FFT program uses a general purpose display routine,
and the number of samples is based on the amount of memory.  This version can
also be used in MS-DOS (if you change the file modes to "rb" and "wb").

The original Byte Magazine program (see references below) was designed for real
data only.  In my experiments I needed to preserve both real and imaginary
data.  If you feed the FFT real data only, then the output will be a mirror
image, and you can ignore the left side.  Two signal generators are included.
One generates sine waves (sine) and the other generates pulses (pulse).  Some
papers I found on the subject of FFTs are included at the end.  There are
several books devoted to the subject also.

For the Unix example try:

	sine 16 in
	1000
	3000

Which will sample the 1 Khz data every 333 microseconds (1 / 3 Khz).  Note: The
sample frequency should be greater than 2 times the input frequency (Nyquist
and all that...).

Then run fft:

	fft 16 in out

And you should see a display like so:

0	|=======				      (-1500.0 Hz)
1	|=====					      (-1312.5 Hz)
2	|====					      (-1125.0 Hz)
3	|====					       (-937.0 Hz)
4	|===					       (-750.0 Hz)
5	|===					       (-562.5 Hz)
6	|===					       (-375.0 Hz)
7	|===					       (-187.5 Hz)
8	|====		<-------   DC			(000.0 Hz)
9	|====		<-------   Fundamental		(187.5 Hz)
10	|======		<-------   Second Harmonic	(375.0 Hz)
11	|========					(562.5 Hz)
12	|==============					(750.0 Hz)
13	|========================================================
14	|============================		    (1125.0 Hz)	^
15	|===========				    (1312.5 Hz)	|
								|
				[13 - 8 (center)] * 187.5 = 937.0 Hz

The fundamental display frequency is:

	T  = Time Increment Between Samples
	N  = Number Of Samples
	Tp = N * T

	Then F = 1 / Tp

	In the example above, the time increment between samples is
	1 / 3000 or 333 microseconds.  N = 16, so Tp = 5333 microseconds
	and 1 / .005333 is 187.5 Hz.

	Therefore each filter is a multiple of 187.5 Hertz.  Filter 8 in this
	example is center, so that would be zero, 9 would be one, etc.

In this case the sample interval didn't work so good for the frequency and
shows the limitation of the Discrete Fourier Transform in representing a
continuous signal.  A better sample rate for 1000 Hz would be 4000 Hz,
in which case T = 250 ms, Tp = 4 ms, and F = 250 Hz.  1000 / 250 = 4.  The
power should all be in filter 12 (8 + 4) in this example.

Let's run it and see:

	sine 16 in
	1000
	4000

	fft 16 in out

0	|
1	|
2	|
3	|
4	|
5	|
6	|
7	|
8	|
9	|
10	|
11	|
12	|========================================================
13	|
14	|
15	|

Well what do you know...

The output data file can be used by other programs as needed.

By using negative frequencies in 'sine' you can generate opening targets:

	sine 16 in
	-1000
	3000

	fft 16 in out

Produces:

0	|=======
1	|===========
2	|============================
3	|=======================================================
4	|==============
5	|========
6	|======
7	|====
8	|====		<-------- Zero Hertz (DC)
9	|===
10	|===
11	|===
12	|===
13	|====
14	|====
15	|=====

You can see in these examples where weighting functions would be nice.

For an example of what happens when the imaginary data is not input
(ie, zeros put in) for a 1000 Hz frequency at 3000 Hz sample rate:

0	|===============
1	|==================
2	|===================================
3	|========================================================
4	|===========
5	|====
6	|==
7	|=				Delete this part
---------------------------------------------------------------------
8	|
9	|=
10	|==
11	|====
12	|===========
13	|=======================================================
14	|===================================
15	|==================

The left side is redundant and can be deleted.  This is what the original
Byte Magazine article did.

For generating pulses, a second program 'pulse' is provided.  It pre-loads
imaginary data with zeros.   For example:

	pulse 16 in
	.000006
	.0000008

Is a radar with a 6 microsecond pulse and 800 nanosecond range gates.

	fft 16 in out

Will produce:

0	|
1	|=======
2	|
3	|========
4	|
5	|============
6	|
7	|===================================
8	|========================================================
9	|===================================
10	|
11	|============
12	|
13	|========
14	|
15	|=======

The more filters you use, the prettier it looks.


FFT References
--------------

1.	Fast Fourier Transforms On Your Home Computer, William D. Stanley,
	Steven J. Peterson, BYTE Magazine, December 1978.  Basic idea comes
	from this program.

2.	8052 Microcomputer simplifies FFT Design, Arnold Rosenberg,
	Electronics, May 5, 1983.  Used a bit reverse table idea based on the
	routine in this program.

3.	A Fast Fourier Transform for the 8080, Robert D. Fusfeld,
	Dr. Dobbs, Number 44.  Gave me some ideas.

4.	A Guided Tour of the Fast Fourier Transform, G. D. Bergland,
	IEEE Spectrum, July 1969.  Math!

5.	FFT - Shortcut to Fourier Analysis, Richard Klahn, Richard R. Shively,
	Electronics, April 15 1968. Math!

/* EOF */