c It runs and computes correctly although it is not in any c sense sophisticated. However, it is free, and useful if c you are in a hurry and don't want to think things out for yourself. c c Arthur Wouk (wouk@brl-vgr) c c ******************* c Fast Fourier Transform c ******************* c SUBROUTINE FFT (X, N, K) C FFT COMPUTES THE (FAST) FOURIER TRANSFORM OF THE VECTOR X C (A COMPLEX ARRAY OF DIMENSION N). SOURCE: Ferziger; Numerical C methods for engineering applications. C C X = DATA TO BE TRANSFORMED; ON RETURN IT CONTAINS THE TRANSFORM. C N = SIZE OF VECTOR. MUST BE A POWER OF 2 (<32769). C K = 1 FOR FORWARD TRANSFORM. C K = -1 FOR INVERSE TRANSFORM. C IMPLICIT INTEGER (A-Z) INTEGER SBY2,S REAL GAIN, PI2, ANG, RE, IM COMPLEX X(N), XTEMP, T, U(16), V, W LOGICAL NEW DATA PI2,GAIN,NO,KO /6.283185307, 1., 0, 0/ C C TEST FIRST CALL? C NEW = ( NO .NE. N) IF ( .NOT. NEW) GO TO 2 C C IF FIRST CALL COMPUTE LOG2 (N). C L2N = 0 NO = 1 1 L2N = L2N + 1 NO = NO + NO IF (NO .LT. N) GO TO 1 GAIN = 1./N ANG = PI2*GAIN RE = COS (ANG) IM = SIN (ANG) C C COMPUTE COMPLEX EXPONENTIALS IF NOT FIRST CALL C 2 IF (.NOT. NEW .AND. K*KO .GE. 1) GO TO 4 U(1) = CMPLX (RE, -SIGN(IM, FLOAT(K))) DO 3 I = 2,L2N U(I) = U(I-1)*U(I-1) 3 CONTINUE K0 = K C C MAIN LOOP C 4 SBY2 = N DO 7 STAGE = 1,L2N V = U(STAGE) W = (1., 0.) S = SBY2 SBY2 = S/2 DO 6 L = 1,SBY2 DO 5 I = 1,N,S P = I + L- 1 Q = P + SBY2 T = X(P) + X(Q) X(Q) = ( X(P) - X(Q))*W X(P) =T 5 CONTINUE W = W*V 6 CONTINUE 7 CONTINUE C C REORDER THE ELEMENTS BY BIT REVERSAL C DO 9 I = 1,N INDEX = I-1 JNDEX = 0 DO 8 J = 1,L2N JNDEX = JNDEX+JNDEX ITEMP = INDEX/2 IF (ITEMP+ITEMP .NE. INDEX) JNDEX = JNDEX + 1 INDEX = ITEMP 8 CONTINUE J = JNDEX + 1 IF (J .LT. I) GO TO 9 XTEMP = X(J) X(J) = X(I) X(I) = XTEMP 9 CONTINUE C C FORWARD TRANSFORM DONE C IF (K .GT. 0) RETURN C C INVERSE TRANSFORM C DO 10 I = 1,N X(I) = X(I)*GAIN 10 CONTINUE RETURN END