SUBROUTINE FFT (Z,M,IWK) ************************************************************************ C This FFT subroutine is taken from the NASA LaRC COMPUTER MANUAL - C VOL II Section E2.4 (1975). C C C FUNCTION - COMPUTE THE FAST FOURIER TRANSFORM, GIVEN A C COMPLEX VECTOR OF LENGTH EQUAL TO A POWER C OF TWO C USAGE - CALL FFT (Z,M,IWK) C PARAMETERS Z - COMPLEX VECTOR OF LENGTH N=2**M C WHICH CONTAINS ON INPUT THE C DATA TO BE TRANSFORMED. ON C OUTPUT,A CONTAINS THE FOURIER C COEFFICIENTS. C M - N = 2**M IS THE NUMBER OF DATA POINTS. C M= +N FFT WILL PERFORM FOURIER C TRANSFORM. C M= -N FFT WILL PERFORM INVERSE C TRANSFORM. C IWK - WORK AREA VECTOR OF LENGTH M+1. C PRECISION - SINGLE C LANGUAGE - FORTRAN C LATEST REVISION - APRIL 16, 1980 C----------------------------------------------------------------------- C DIMENSION IWK(*),Z(*),Z0(2),Z1(2),Z2(2),Z3(2) COMPLEX Z,ZA0,ZA1,ZA2,ZA3,AK2 EQUIVALENCE (ZA0,Z0(1)),(ZA1,Z1(1)),(ZA2,Z2(1)),(ZA3,Z3(1)) *, (A0,Z0(1)),(B0,Z0(2)),(A1,Z1(1)),(B *1,Z1(2)), (A2,Z2(1)),(B2,Z2(2)),(A3,Z3(1)),(B *3,Z3(2)) DATA SQ,SK,CK /.70710678118655,.38268343236509, * .92387953251129/ DATA TWOPI/6.2831853071796/,ZERO/0.0/,ONE/1.0/ C SQ=SQRT2/2,SK=SIN(PI/8),CK=COS(PI/8) C TWOPI=2*PI SIGN = 1.0 IF (M .LT. 0) SIGN = -1.0 M= IABS(M) MP = M+1 N = 2**M IF (SIGN .EQ. 1.0) GO TO 4 DO 2 I=1,N Z(I) = CONJG(Z(I)) 2 CONTINUE 4 IWK(1)=1 MM = (M/2)*2 KN = N+1 C INITIALIZE WORK VECTOR DO 5 I=2,MP IWK(I) = IWK(I-1)+IWK(I-1) 5 CONTINUE RAD=TWOPI/N MK = M - 4 KB = 1 IF (MM .EQ. M) GO TO 15 K2 = KN K0 = IWK(MM+1) + KB 10 K2 = K2 - 1 K0 = K0 - 1 AK2 = Z(K2) Z(K2) = Z(K0)- AK2 Z(K0) = Z(K0)+ AK2 IF (K0 .GT. KB) GO TO 10 15 C1 = ONE S1 = ZERO JJ = 0 K = MM - 1 J = 4 IF (K .GE. 1) GO TO 30 GO TO 9005 20 IF (IWK(J) .GT. JJ) GO TO 25 JJ = JJ - IWK(J) J = J-1 IF (IWK(J) .GT. JJ) GO TO 25 JJ = JJ - IWK(J) J = J - 1 K = K + 2 GO TO 20 25 JJ = IWK(J) + JJ J = 4 30 ISP = IWK(K) IF (JJ .EQ. 0) GO TO 40 C RESET TRIGONOMETRIC PARAMETERS C2 = JJ * ISP * RAD C1 = COS(C2) S1 = SIN(C2) 35 C2 = C1 * C1 - S1 * S1 S2 = C1 * (S1 + S1) C3 = C2 * C1 - S2 * S1 S3 = C2 * S1 + S2 * C1 40 JSP = ISP + KB C DETERMINE FOURIER COEFFICIENTS C IN GROUPS OF 4 DO 50 I=1,ISP K0 = JSP - I K1 = K0 + ISP K2 = K1 + ISP K3 = K2 + ISP ZA0 = Z(K0) ZA1 = Z(K1) ZA2 = Z(K2) ZA3 = Z(K3) IF (S1 .EQ. ZERO) GO TO 45 TEMP = A1 A1 = A1 * C1 - B1 * S1 B1 = TEMP * S1 + B1 * C1 TEMP = A2 A2 = A2 * C2 - B2 * S2 B2 = TEMP * S2 + B2 * C2 TEMP = A3 A3 = A3 * C3 - B3 * S3 B3 = TEMP * S3 + B3 * C3 45 TEMP = A0 + A2 A2 = A0 - A2 A0 = TEMP TEMP = A1 + A3 A3 = A1 - A3 A1 = TEMP TEMP = B0 + B2 B2 = B0 - B2 B0 = TEMP TEMP = B1 + B3 B3 = B1 - B3 B1 = TEMP Z(K0) = CMPLX(A0+A1,B0+B1) Z(K1) = CMPLX(A0-A1,B0-B1) Z(K2) = CMPLX(A2-B3,B2+A3) Z(K3) = CMPLX(A2+B3,B2-A3) 50 CONTINUE IF (K .LE. 1) GO TO 55 K = K - 2 GO TO 30 55 KB = K3 + ISP C CHECK FOR COMPLETION OF FINAL C ITERATION IF (KN .LE. KB) GO TO 9005 IF (J .NE. 1) GO TO 60 K = 3 J = MK GO TO 20 60 J = J - 1 C2 = C1 IF (J .NE. 2) GO TO 65 C1 = C1 * CK + S1 * SK S1 = S1 * CK - C2 * SK GO TO 35 65 C1 = (C1 - S1) * SQ S1 = (C2 + S1) * SQ GO TO 35 9005 CONTINUE IF (SIGN .EQ. 1.0) GO TO 75 XN = N DO 70 I=1,N Z(I)=CONJG(Z(I))/XN 70 CONTINUE 75 CALL QXZ136(Z,M,IWK) IF(SIGN.EQ.-1.0) M = -M RETURN END C --- SUBPROGRAM QXZ136 --- FORMERLY KNOWN AS ROUTINE FFRDR2 --- SUBROUTINE QXZ136 (Z,M,IWK) C-FFRDR2--------S-------LIBRARY 3--------------------------------------- C ************************************************************************ C C FUNCTION - THIS SUBROUTINE PERMUTES A COMPLEX DATA VECTOR C IN REVERSE BINARY ORDER TO NORMAL ORDER. THE C ROUTINE CAN ALSO BE USED TO PERMUTE A COM- C PLEX DATA VECTOR IN NORMAL ORDER TO REVERSE C BINARY ORDER SINCE THE PERMUTATION IS SYM- C METRIC. C USAGE - CALL QXZ136(Z,M,IWK) C PARAMETERS Z - COMPLEX VECTOR OF LENGTH N=2**M WHICH C CONTAINS ON INPUT THE DATA TO BE C PERMUTED. ON OUTPUT, Z CONTAINS THE C PERMUTED DATA VECTOR. C M - N=2**M IS THE NUMBER OF DATA POINTS. C IWK - WORK AREA VECTOR OF LENGTH M+1 C PRECISION - SINGLE C LANGUAGE - FORTRAN C LATEST REVISION - MARCH 16, 1973 C ----------------------------------------------------------------- C DIMENSION IWK(*),Z(*) COMPLEX Z,TEMP C IF(M .LE. 1) GO TO 9005 MP = M+1 JJ = 1 C INITIALIZE WORK VECTOR IWK(1) = 1 DO 5 I = 2,MP IWK(I) = IWK(I-1) * 2 5 CONTINUE N4 = IWK(MP-2) IF (M .GT. 2) N8 = IWK(MP-3) N2 = IWK(MP-1) LM = N2 NN = IWK(MP)+1 MP = MP-4 C DETERMINE INDICES AND SWITCH A*S J = 2 10 JK = JJ + N2 TEMP= Z(J) Z(J)=Z(JK) Z(JK)= TEMP J = J+1 IF (JJ .GT. N4) GO TO 15 JJ = JJ + N4 GO TO 35 15 JJ = JJ - N4 IF (JJ .GT. N8) GO TO 20 JJ = JJ + N8 GO TO 35 20 JJ = JJ - N8 K = MP 25 IF (IWK(K) .GE. JJ) GO TO 30 JJ = JJ - IWK(K) K = K - 1 GO TO 25 30 JJ = IWK(K) + JJ 35 IF (JJ .LE. J) GO TO 40 K = NN - J JK = NN - JJ TEMP= Z(J) Z(J) = Z(JJ) Z(JJ) = TEMP TEMP = Z(K) Z(K) = Z(JK) Z(JK)= TEMP 40 J = J + 1 C CYCLE REPEATED UNTIL LIMITING NUMBER C OF CHANGES IS ACHIEVED IF (J .LE. LM) GO TO 10 9005 RETURN END