PROCEDURE FFT (VAR Z: ARRAYCOMPLEX; LAYER, SIGN: INTEGER);
{PASCAL version of FFT. FFTEST.PAS is a demonstration program. See also
FFT.MAC and FFTEST.MAC.
TABLE OF CONTENTS
1. DIRECTIONS
2. SLOW FOURIER TRANSFORM
3. FAST FOURIER TRANSFORM
4. PASCAL VERSION OF FFT
1. DIRECTIONS
The main program must define:
TYPE COMPLEX = ARRAY [1..2] OF REAL;
TYPE ARRAYCOMPLEX = ARRAY [0..1023] OF COMPLEX; M = 2^10 POINTS
%INCLUDE 'FFT'; OR OTHER DIMENSION
The original Zarray is destroyed; the transform is in the Zarray.
Z[0] DC component
Z[ 1..M/2 ] positive freq terms: 1/T, 2/T ...
Z[ M/2 ] Nyquist freq term
Z[ M/2+1...M1] negative freq terms
A scale factor of 1/M is applied only to the forward transform.
2. SLOW FOURIER TRANSFORM
The Fourier transform is a method for representing a timedependent
signal as an integral over sinusoids of various frequencies. See
standard texts such as A.V. Oppenheim and A.S. Willsky, Signals and
Systems, PrenticeHall, Englewood Cliffs, NJ.
The following constants are defined:
j = SQRT(1)
LAYER = number of layers in the FFT
M = number of points = 2^LAYER (^ means exponentiation)
DT = sampling interval in seconds
T = M*DT = total sampling period
F = 1/DT = sampling freq. in Hz
DF = 1/T = freq. step size in Hz
NYQ = M/2 = number of the Nyquistfreq. point
FNYQ = Nyquist freq. = 0.5*F Hz
W = principle Mth root of unity = "twiddle factor"
= EXP (2*PI*j/M) = COS (2*PI/M) + j * SIN (2*PI/M)
Example: What are the 3 cube roots of 1?
Note that 1 = EXP (0) = EXP(2*PI*j) = EXP(4*PI*j). Thus the roots are:
EXP(0) = 1
EXP (2*PI*j/3) = COS (120 deg) + j*SIN (120 deg) = 0.5 + 0.866*j
EXP (4*PI*j/3) = COS (240) + j*SIN(240) = 0.5  0.866*J
In the complex plane, these lie on the unit circle spaced 120 deg apart.
Signals may be represented either in terms of sines and cosines, or
in terms of complex exponentials; the latter are used by the FFT
algorithm because the math is simpler.
The two methods are equivalent because of this theorem:
EXP(j*X) = COS(X) + j*SIN(X)
Theorems: From the above relationship, one can prove:
COS(X) = 0.5*(EXP(j*X) + EXP(j*X)) and
SIN(X) = 0.5*j*(EXP(j*X)  EXP(j*X))
The signal has been measured (with an A/D converter) at the times:
0, DT, 2*DT .. (M1)*DT, that is, at T[K] = K*DT, K = 0..(M1)
The signal magnitudes at these times are the smallg array, which are
complex numbers, whose imaginary parts are zero, because the signal is real.
g[0], g[1] .. g[M1]
The signal is represented by a sum of complex exponentials at the
following frequencies (Hz):
0, DF, 2*DF .. (M1)*DF, that is, at F[L] = L*DF, L = 0..(M1)
These frequencies are subdivided into the following categories:
L = 0 DC term
L = 1..(M/21) positive frequencies
L = M/2 Nyquist frequency
L = (M/2+1)..(M1) "negative frequencies"
The forward transform is defined:
G[L] = (1/M)* sum K=0 to M1 of g[K] * EXP(2*PI*j*F[L]*T[K])
= (1/M)* sum K=0 to M1 of g[K] * EXP(2*PI*j*L*DF*K*DT)
= (1/M)* sum K=0 to M1 of g[K] * EXP(2*PI*j*L*K/M) because DF*DT=1/M
= (1/M)* sum K=0 to M1 of g[K] * W^(K*L), (1)
because W = EXP(2*PI*j/M). The capitalG array contains the Fourier
coefficients. The minus sign is used for the forward transform, while a
plus sign will be used for the inverse transform below. The 2*PI is used
to convert Hz into radians/s. The 1/M is used for normalization.

Theorem: The inverse transform is:
g[I] = sum L=0 to M1 of G[L]*W^(I*L) (2)
The inverse transform is the same as the forward transform, except for the
omission of the 1/M normalizing factor, and the plus sign in the exponent.
Proof: Substitute (1) into (2):
g[I]=sum L=0 to M1 (1/M)* sum K=0 to M1 g[K] * W^(K*L) * W^(I*L) (3)
Combine the two exponents of W, and reverse the order of summation:
g[I] = (1/M) * sum K=0 to M1 g[K] * sum L=0 to M1 of W^(L*(IK)) (4)
It will now be proven that:
sum L=0 to M1 of W^(L*(IK)) = M if IK=0 or (5)
= 0 otherwise
First, if IK=0, then W^0 = 1 and the summation gives a total of M.
Second, if IK = an integer other than zero, say N, then:
W^(L*(IK)) = W^(L*N) = (W^N)^L
The result can be seen most easily by graphing the values of (W^N)^L
in the complex plane, for all the values of L in the summation. For example,
for M=8 and N=1:
j x L=2
L=3 x  x L=1


L=4 
xx L=0
1  1


L=5 x  x L=7
j x L=6
The x's fall on a unit circle. Obviously the vector sum of all 8 points
is zero. This nonrigorous graphical method can be used to prove the
theorem (5) for all values of M and N. Thus in eq. 4, the inner sum is
zero except when I=K, when the sum equals M. Thus:
g[I] = g[K] when I=K QED.

What is the meaning of a "negative frequency" term? Consider L = M1,
which is called the first negative frequency. Consider a complex sinusoid
whose frequency is (M1)*DF, and is a function of time K*DT:
EXP (2*PI*j* (M1)*DF * K*DT) Use DF*DT = 1/M
= EXP (2*PI*j* (M1)*K/M)
= EXP (2*PI*j*M*K/M) * EXP( 2*PI*j*K/M) Replace 1/M by DF*DT
= 1^K * EXP( 2*PI*j (DF) * K*DT)
This behaves like a complex sinusoid of frequency DF.

Theorem: For a real signal having all g[K] real, the negativefrequency
number G[ML] is the complex conjugate of the corresponding positivefrequency
number G[L].
Proof: By eq. 1: G[L] = (1/M)* sum K=0 to M1 of g[K] * W^(K*L). By def. W:
= (1/M)* sum K=0 to M1 of g[K]*(COS(2*PI*K*L/M)  j*SIN(2*PI*K*L/M)) (3)
Also: G[ML] = (1/M)* sum K=0 to M1 of g[K] * W^(K*(ML))
= (1/M)* sum K=0 to M1 of g[K] * W^(K*M) * W^(K*L)
= (1/M)* sum K=0 to M1 of g[K] * (W^M)^(K) * W^(K*L)
(Because W is the Mth root of unity, W^M = 1. Use the definition of W.)
= (1/M)* sum K=0 to M1 of g[K] * EXP (2*PI*j*K*L/M)
=(1/M)* sum K=0 to M1 of g[K] * (COS(2*PI*K*L/M) + j* SIN(2*PI*K*L/M)
This is the complex conjugate of (3) QED.
3. FAST FOURIER TRANSFORM
The FFT algorithm by J.W. Cooley and J.W. Tukey (An algorithm for the
machine computation of complex Fourier series, Math. Comput.
19(1965)297) saves computer time by avoiding the need to calculate the
powers of W. For an array of M=1024 points, for instance, the ordinary
Fourier transform (eq. 1) requires about M^2 ~= 1,000,000 complex
multiplications and additions, whereas the FFT algorithm requires only
2*M*log(M) to base 2 = 20,000, a savings of 50fold!
The FFT is done on a complex array Z with M = 2^LAYER points,
where LAYER is a small integer.
The FFT algorithm can be partially understood by comparing the slow
transform (1) with a manual execution of the PASCAL program. We write out
explicitly the FFTs for M = 2, 4 and 8. We do the inverse transform with
the + sign so that we do not have to carry the sign along.
Rules for a "butterfly" calculation: A A+W*B
\ /
W
/ \
B AW*B
The numbers A and B are in designated memory locations. They are
replaced IN THE SAME MEMORY LOCATIONS by A+W*B and AW*B, respectively.
The upper term has a + sign, while the lower has . The W represents a
specific power of W, as shown below.
M = 2. W = EXP(2*PI*j/2) = 1. Slow transform (eq. 2):
g0 = G0 + G1
g1 = G0 + W*G1 = G0  G1
Fast transform:
G0 G0+G1
\ /
W^0
/ \ Slow and fast transforms agree. This has 1
G1 G0G1 layer (1 butterfly).

M = 4. W = j. The slow transform produces polynomials in W.
g0 = G0 + G1 + G2 + G3
g1 = G0 + G1*W + G2*W^2 + G3*W^3 = G0 + j*G1 G2  j*G3
g2 = G0 + G1*W^2 + G2*W^4 + G3*W^6 = G0  G1 + G2  G3
g3 = G0 + G1*W^3 + G2*W^6 + G3*W^9 = G0  j*G1  G2 + j*G3
Fast transform:
G0 G0+G2 G0+G2+G1+G3
\ / \ /
W^0 \ /
/ \ \ /
G2 G0G2 W^0 G0G2+j*(G1G3)
\ / \ /
\/ \/
/\ /\
/ \ / \
G1 G1+G3 W^1 G0+G2G1G3
\ / / \
W^0 / \
/ \ / \ Fast and slow transforms agree.
G3 G1G3 G0G2j*(G1G3)
M=8; W = EXP (2*PI*j/8) = 0.707 + j*0.707; Slow transform:
g0 = G0 + G1 + G2 + G3 + G4 + G5 + G6 + G7
g1 = G0 + G1*W + G2*W^2 + G3*W^3 + G4*W^4 + G5*W^5 + G6*W^6 + G7*W^7
= G0 + G1*W + G2*j + G3*W^3  G4  G5*W  G6*j  G7*W^3
g2 = G0 + G1*W^2 + G2*W^4 + G3*W^6 + G4*W^8 + G5*W^10 + G6*W^12 + G7*W^14
= G0 + G1*j  G2  G3*j + G4 + G5*j  G6  G7*j
g3 = G0 + G1*W^3 + G2*W^6 + G3*W^9 + G4*W^12 + G5*W^15 + G6*W^18 + G7*W^21
= G0 + G1*W^3  G2*j + G3*W  G4  G5*W^3 + G6*j  G7*W
g4 = G0 + G1*W^4 + G2*W^8 + G3*W^12 + G4*W^16 + G5*W^20 + G6*W^24 + G7*W^28
= G0  G1 + G2  G3 + G4  G5 + G6  G7
g5 = G0 + G1*W^5 + G2*W^10 + G3*W^15 + G4*W^20 + G5*W^25 + G6*W^30 + G7*W^35
= G0  G1*W + G2*j  G3*W^3  G4 + G5*W  G6*j + G7*W^3
g6 = G0 + G1*W^6 + G2*W^12 + G3*W^18 + G4*W^24 + G5*W^30 + G6*W^36 + G7*W^42
= G0  G1*j  G2 + G3*j + G4  G5*j  G6 + G7*j
g7 = G0 + G1*W^7 + G2*W^14 + G3*W^21 + G4*W^28 + G5*W^35 + G6*W^42 + G7*W^49
= G0  G1*W^3  G2*j  G3*W  G4 + G5*W^3 + G6*j + G7*W
There are too many powers to compute. Note that only powers up to W^3 have
to be computed because W^4 = 1 and W^8 = 1. Fast transform:
G0 G0+G4 G0+G4+G2+G6 G0+G4+G2+G6+G1+G5+G3+G7 =g0
\ / \ / \ /
W^0 \ / \ /
/ \ \ / \ /
G4 G0G4 W^0 G0G4+j(G2G6)\ / G0G4+j(G2G6)+W(G1G5)+W^3(G3G7)=g1
\/ \/ \ \ / /
/\ /\ \ \ / /
/ \ / \ \ \ / /
G2 G2+G6 W^2 G0+G4G2G6 \ W^0 / G0+G4G2G6+j(G1+G5G3G7) =g2
\ / / \ \ \ / \ / /
W^0 / \ \ \/ \/ /
/ \ / \ \ /\ /\ /
G6 G2G6 G0G4j(G2G6) \/ W^1 \/ G0G4+W(G3G7)j(G2G6)+W^3(G1G5)=g3
\ /\ / \ /\ /
\/ \/ \/ \/
/\ /\ /\ /\
G1 G1+G5 G1+G5+G3+G7 \/ W^2 \/ G0+G4+G2+G6G1G5G3G7 =g4
\ / \ / /\ / \ /\
W^0 \ / / \/ \/ \
/ \ \ / / /\ /\ \
G5 G1G5 W^0 G1G5+j(G3G7)/ W^3 \ G0G4+J(G2G6)W(G1G5)W^3(G3G7)=g5
\/ \ / / / \ \
/\ \/ / / \ \
/ \ /\ / / \ \
G3 G3+G7 W^2 G1+G5G3G7 / \ G0+G4G2G6j(G1+G5G3G7) =g6
\ / / \ / \
W^0 / \ / \
/ \ / \ / \
G7 G3G7 G1G5j(G3G7) G0G4j(G2G6)W^3(G1G5)W(G3G7)=g7
The slow and fast transforms agree.
4. PASCAL VERSION OF FFT:}
VAR Q, W, HOLD: COMPLEX;
L,I,II,J,K,LX,NBLOCK,LBLOCK,LBHALF,NW,IBLOCK,LSTART: INTEGER;
FLAG: BOOLEAN;
SCALE, FLXPI2: REAL;
E2: ARRAY [0..15] OF INTEGER;
BEGIN
E2[0]:= 1;
FOR I:= 1 TO 15 DO E2[I]:= 2*E2[I1]; {TABULATE POWERS OF 2}
LX:= E2[LAYER];
FLXPI2:= (6.28318530 * SIGN)/LX;
FOR L:= 1 TO LAYER DO
BEGIN
NBLOCK:= E2[L1];
LBLOCK:= LX DIV NBLOCK;
LBHALF:= LBLOCK DIV 2;
NW:= 0;
FOR IBLOCK:=1 TO NBLOCK DO
BEGIN
LSTART:= LBLOCK*(IBLOCK1);
W[1]:= COS(FLXPI2*NW); {W=EXP(FLXPI2*NW j)}
W[2]:= SIN(FLXPI2*NW);
FOR I:=0 TO LBHALF1 DO
BEGIN
J:= I + LSTART;
K:= J + LBHALF;
Q[1]:= Z[K,1]*W[1]  Z[K,2]*W[2]; {COMPLEX MULT}
Q[2]:= Z[K,1]*W[2] + Z[K,2]*W[1];
Z[K,1]:= Z[J,1]  Q[1]; {COMPLEX SUBTRACT}
Z[K,2]:= Z[J,2]  Q[2];
Z[J,1]:= Z[J,1] + Q[1]; {COMPLEX ADD}
Z[J,2]:= Z[J,2] + Q[2]
END;
IF LAYER > 1 THEN
BEGIN
I:=2; {LOOP COUNTER}
FLAG:=FALSE;
WHILE ( (I <= LAYER) AND NOT FLAG) DO
BEGIN
II:=I; {IN INTEGER NW,IS ITH BIT FROM LEFT=0?}
IF NW < E2[LAYERI]
THEN FLAG:= TRUE
ELSE NW:= NWE2[LAYERI]; {MAKE ITH BIT=0}
I:= I + 1
END;
NW:= NW + E2[LAYERII] {MAKE ITH BIT ONE}
END;
END;
END;;
NW:=0;
FOR K:=0 TO LX1 DO
BEGIN
IF NW+1K > 0 THEN
BEGIN
HOLD:= Z[NW]; {SWITCH COMPLEX NUMBERS}
Z[NW]:= Z[K];
Z[K]:= HOLD
END;
I:=1;
FLAG:= FALSE;
WHILE ((I <= LAYER) AND NOT FLAG) DO
BEGIN
II:= I;
IF NW < E2[LAYERI] THEN FLAG:=TRUE
ELSE NW:= NW  E2[LAYERI];
I:= I +1;
END;
NW:= NW + E2[LAYERII];
END;
IF SIGN = 1 THEN {SCALE ONLY FORWARD TRANSFORM}
BEGIN
SCALE:= 1/E2[LAYER];
FOR I:=0 TO LX1 DO
BEGIN
Z[I,1]:= Z[I,1]*SCALE;
Z[I,2]:= Z[I,2]*SCALE;
END;
END;
END;