|
996
|
1 The program fft.i provides a basic complex Fast Fourier Transform,
|
|
|
2 along with a simple write-in-read-out wrapper. This program requires
|
|
|
3 the floating-point library in lib/floatlib.i, which should be appended
|
|
|
4 to fft.i before compiling.
|
|
|
5
|
|
|
6 The first number to enter specifies the direction of the transform.
|
|
|
7 Enter ONE to transform from time to frequency, or TWO for the inverse
|
|
|
8 transform. The second number to enter is the number of complex pairs.
|
|
|
9 This number must be a power of two, between TWO and THREE TWO SEVEN
|
|
|
10 SIX EIGHT inclusive. After this comes the actual data. Each entry must
|
|
|
11 consist of two floating-point numbers, the real value followed by the
|
|
|
12 imaginary value.
|
|
|
13
|
|
|
14 After the data has been entered, the FFT is performed on it. When it
|
|
|
15 is completed, the resulting data will be output. The data is preceded
|
|
|
16 by I or II to indicate time or frequency data (i.e., the opposite of
|
|
|
17 the number input at the start), and a second integer giving the
|
|
|
18 number of pairs being output.
|
|
|
19
|
|
|
20 Included are two files that produce simple curves when transformed.
|
|
|
21 delta.fft, the sum of two delta functions, transforms into a cosine
|
|
|
22 wave, and tophat.fft's output traces the sinc function.
|
|
|
23
|
|
|
24 INTERCAL brings new meaning to the "fast" in Fast Fourier Transform.
|
|
|
25 On Muppetlabs' Pentium, 133MHz, a transform of 2048 data points takes
|
|
|
26 over 5 minutes of CPU time.
|
|
|
27
|
|
|
28
|
|
|
29 The actual FFT routine is at line label (100), and is self-contained
|
|
|
30 other than its dependence on the standard and floating-point
|
|
|
31 libraries, and a small subroutine described below. Its three arguments
|
|
|
32 are .9, .10, and ;1. .9 must contain #1 or #2 to select the forward or
|
|
|
33 inverse transform, respectively. .10 must contain the length of the
|
|
|
34 input data, and ;1 the data itself. ;1 must have the dimensions .10 BY
|
|
|
35 #2. No error checking is done on the parameters within the routine.
|
|
|
36 (The wrapper routine sanity-checks the user's input.) On exit, .9 will
|
|
|
37 contain #2 if it was #1 originally and vice versa, and ;1 will contain
|
|
|
38 the transformed data.
|
|
|
39
|
|
|
40 The FFT algorithm was taken from "Numerical Recipes" (Press, et al;
|
|
|
41 Cambridge University Press, 1994 edition). For data centered on the
|
|
|
42 y-axis, it is usually best to cut the data in half at the axis and
|
|
|
43 switch their positions (i.e., enter the data for 0 to N/2 - 1 first,
|
|
|
44 followed by the data for -N/2 to -1) and then do the same with the
|
|
|
45 output.
|
|
|
46
|
|
|
47 The FFT routine uses a subroutine, labelled (153), as a wrapper around
|
|
|
48 the multiplication routine (5030). This subroutine is an example of
|
|
|
49 how to handle errors as indicated by the .5 variable.
|
|
|
50
|
|
|
51 Another noteworthy point is in the first part of the routine, where
|
|
|
52 the positions of the data are swapped with the bit-reversals of their
|
|
|
53 indices. Amazingly, bit-reversal is one of the few things that
|
|
|
54 INTERCAL can do more tersely than most other languages. It is nice to
|
|
|
55 see it actually finding a real-world application.
|
