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 Z-array is destroyed; the transform is in the Z-array. 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...M-1] 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 time-dependent signal as an integral over sinusoids of various frequencies. See standard texts such as A.V. Oppenheim and A.S. Willsky, Signals and Systems, Prentice-Hall, 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 Nyquist-freq. 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 .. (M-1)*DT, that is, at T[K] = K*DT, K = 0..(M-1) The signal magnitudes at these times are the small-g array, which are complex numbers, whose imaginary parts are zero, because the signal is real. g[0], g[1] .. g[M-1] The signal is represented by a sum of complex exponentials at the following frequencies (Hz): 0, DF, 2*DF .. (M-1)*DF, that is, at F[L] = L*DF, L = 0..(M-1) These frequencies are subdivided into the following categories: L = 0 DC term L = 1..(M/2-1) positive frequencies L = M/2 Nyquist frequency L = (M/2+1)..(M-1) "negative frequencies" The forward transform is defined: G[L] = (1/M)* sum K=0 to M-1 of g[K] * EXP(-2*PI*j*F[L]*T[K]) = (1/M)* sum K=0 to M-1 of g[K] * EXP(-2*PI*j*L*DF*K*DT) = (1/M)* sum K=0 to M-1 of g[K] * EXP(-2*PI*j*L*K/M) because DF*DT=1/M = (1/M)* sum K=0 to M-1 of g[K] * W^(-K*L), (1) because W = EXP(2*PI*j/M). The capital-G 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 M-1 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 M-1 (1/M)* sum K=0 to M-1 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 M-1 g[K] * sum L=0 to M-1 of W^(L*(I-K)) (4) It will now be proven that: sum L=0 to M-1 of W^(L*(I-K)) = M if I-K=0 or (5) = 0 otherwise First, if I-K=0, then W^0 = 1 and the summation gives a total of M. Second, if I-K = an integer other than zero, say N, then: W^(L*(I-K)) = 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 | -----x---------------------------x 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 = M-1, which is called the first negative frequency. Consider a complex sinusoid whose frequency is (M-1)*DF, and is a function of time K*DT: EXP (2*PI*j* (M-1)*DF * K*DT) Use DF*DT = 1/M = EXP (2*PI*j* (M-1)*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 negative-frequency number G[M-L] is the complex conjugate of the corresponding positive-frequency number G[L]. Proof: By eq. 1: G[L] = (1/M)* sum K=0 to M-1 of g[K] * W^(-K*L). By def. W: = (1/M)* sum K=0 to M-1 of g[K]*(COS(2*PI*K*L/M) - j*SIN(2*PI*K*L/M)) (3) Also: G[M-L] = (1/M)* sum K=0 to M-1 of g[K] * W^(-K*(M-L)) = (1/M)* sum K=0 to M-1 of g[K] * W^(-K*M) * W^(K*L) = (1/M)* sum K=0 to M-1 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 M-1 of g[K] * EXP (2*PI*j*K*L/M) =(1/M)* sum K=0 to M-1 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 50-fold! 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 A-W*B The numbers A and B are in designated memory locations. They are replaced IN THE SAME MEMORY LOCATIONS by A+W*B and A-W*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 G0-G1 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 G0-G2 W^0 G0-G2+j*(G1-G3) \ / \ / \/ \/ /\ /\ / \ / \ G1 G1+G3 W^1 G0+G2-G1-G3 \ / / \ W^0 / \ / \ / \ Fast and slow transforms agree. G3 G1-G3 G0-G2-j*(G1-G3) 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 G0-G4 W^0 G0-G4+j(G2-G6)\ / G0-G4+j(G2-G6)+W(G1-G5)+W^3(G3-G7)=g1 \/ \/ \ \ / / /\ /\ \ \ / / / \ / \ \ \ / / G2 G2+G6 W^2 G0+G4-G2-G6 \ W^0 / G0+G4-G2-G6+j(G1+G5-G3-G7) =g2 \ / / \ \ \ / \ / / W^0 / \ \ \/ \/ / / \ / \ \ /\ /\ / G6 G2-G6 G0-G4-j(G2-G6) \/ W^1 \/ G0-G4+W(G3-G7)-j(G2-G6)+W^3(G1-G5)=g3 \ /\ / \ /\ / \/ \/ \/ \/ /\ /\ /\ /\ G1 G1+G5 G1+G5+G3+G7 \/ W^2 \/ G0+G4+G2+G6-G1-G5-G3-G7 =g4 \ / \ / /\ / \ /\ W^0 \ / / \/ \/ \ / \ \ / / /\ /\ \ G5 G1-G5 W^0 G1-G5+j(G3-G7)/ W^3 \ G0-G4+J(G2-G6)-W(G1-G5)-W^3(G3-G7)=g5 \/ \ / / / \ \ /\ \/ / / \ \ / \ /\ / / \ \ G3 G3+G7 W^2 G1+G5-G3-G7 / \ G0+G4-G2-G6-j(G1+G5-G3-G7) =g6 \ / / \ / \ W^0 / \ / \ / \ / \ / \ G7 G3-G7 G1-G5-j(G3-G7) G0-G4-j(G2-G6)-W^3(G1-G5)-W(G3-G7)=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[I-1]; {TABULATE POWERS OF 2} LX:= E2[LAYER]; FLXPI2:= (6.28318530 * SIGN)/LX; FOR L:= 1 TO LAYER DO BEGIN NBLOCK:= E2[L-1]; LBLOCK:= LX DIV NBLOCK; LBHALF:= LBLOCK DIV 2; NW:= 0; FOR IBLOCK:=1 TO NBLOCK DO BEGIN LSTART:= LBLOCK*(IBLOCK-1); W[1]:= COS(FLXPI2*NW); {W=EXP(FLXPI2*NW j)} W[2]:= SIN(FLXPI2*NW); FOR I:=0 TO LBHALF-1 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[LAYER-I] THEN FLAG:= TRUE ELSE NW:= NW-E2[LAYER-I]; {MAKE ITH BIT=0} I:= I + 1 END; NW:= NW + E2[LAYER-II] {MAKE ITH BIT ONE} END; END; END;; NW:=0; FOR K:=0 TO LX-1 DO BEGIN IF NW+1-K > 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[LAYER-I] THEN FLAG:=TRUE ELSE NW:= NW - E2[LAYER-I]; I:= I +1; END; NW:= NW + E2[LAYER-II]; END; IF SIGN = -1 THEN {SCALE ONLY FORWARD TRANSFORM} BEGIN SCALE:= 1/E2[LAYER]; FOR I:=0 TO LX-1 DO BEGIN Z[I,1]:= Z[I,1]*SCALE; Z[I,2]:= Z[I,2]*SCALE; END; END; END;