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