ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Radix 4 complex FFT c c Usage: CALL FASTF( XREAL, XIMAG, ISIZE, ITYPE, IFAULT ) c c Arguments: c XREAL - array containing real parts of transform sequence (i/o) c XIMAG - array containing imaginary parts of transform sequence (i/o) c ISIZE - Size of transforms -- ISIZE = 4*2**M (i) c ITYPE - +1 to denote forward transform c -1 to denote backward transform (i) c IFAULT - 1 if error in arguments, 0 otherwise (o) c c Forward transform computes c X(k) = sum_{j=0}^{ISIZE-1} x(j)*exp(-2ijk*pi/ISIZE) c c Backward transform computes c x(j) = (1/ISIZE) * sum_{k=0}^{ISIZE-1} X(k)*exp(2ijk*pi/ISIZE) c c Notice that these transforms are normalized in such a way that a call c to the forward transform followed by a call to the backward transform c will result in the original sequence. c c This is a modified (slightly) version of the original algorithm c referred to below. I believe that the original author was D. Monro, c Imperial College, London. (I have not checked out the reference, but c please let me know if you have.) c c Modifications - c Alg. 83.1 -- Entry error checking was rewritten c Alg. 83.2 -- Normalization was switched to backward transform c Alg. 83.3 -- Unchanged c c Steve Kifowit, kifowit@math.niu.edu, July 1997 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE FASTF( XREAL, XIMAG, ISIZE, ITYPE, IFAULT ) C C ALGORITHM AS 83.1 APPL. STATIST. (1975) VOL.24, P.153 C C RADIX 4 COMPLEX DISCRETE FAST FOURIER TRANSFORM WITH C UNSCRAMBLING OF THE TRANSFORMED ARRAYS C REAL XREAL(*), XIMAG(*) C C CHECK FOR VALID TRANSFORM SIZE - UP TO 2 ** MAX2 C DATA MAX2 /20/ C IFAULT = 1 IF ( ISIZE .LT. 4 ) RETURN II = 4 IPOW = 2 1 IF ( ( II - ISIZE ) .NE. 0 ) THEN II = II * 2 IPOW = IPOW + 1 IF ( IPOW .GT. MAX2 ) RETURN GOTO 1 ENDIF IF ( IABS( ITYPE ) .NE. 1 ) RETURN C C IF WE REACH THIS POINT, THERE ARE NO ENTRY ERRORS C IFAULT = 0 C C CALL FASTG (ALGORITHM AS 83.2) TO PERFORM THE TRANSFORM C CALL FASTG( XREAL, XIMAG, ISIZE, ITYPE ) C C CALL SCRAM (ALGORITHM AS 83.3) C TO UNSCRAMBLE THE RESULTS C CALL SCRAM( XREAL, XIMAG, ISIZE, IPOW ) RETURN c ... End of subroutine FASTF ... END C SUBROUTINE FASTG( XREAL, XIMAG, N, ITYPE ) C C ALGORITHM AS 83.2 APPL. STATIST. (1975) VOL.24, P.153 C C RADIX 4 COMPLEX DISCRETE FAST FOURIER TRANSFORM WITHOUT C UNSCRAMBLING, SUITABLE FOR CONVOLUTIONS OR OTHER C APPLICATIONS WHICH DO NOT REQUIRE UNSCRAMBLING. C SUBROUTINE FASTF USES THIS ROUTINE FOR TRANSFORMATION C AND ALSO PROVIDES UNSCRAMBLING C REAL XREAL(*), XIMAG(*), BCOS, BSIN, CW1, CW2, CW3, PI, $ SW1, SW2, SW3, TEMPR, X1, X2, X3, XS0, XS1, XS2, XS3, $ Y1, Y2, Y3, YS0, YS1, YS2, YS3, Z, ZERO, HALF, ONE, $ ONE5, TWO, FOUR, ZATAN, ZFLOAT, ZSIN C DATA ZERO, HALF, ONE, ONE5, TWO, FOUR $ /0.0, 0.5, 1.0, 1.5, 2.0, 4.0/ C ZATAN(Z) = ATAN(Z) ZFLOAT(K) = FLOAT(K) ZSIN(Z) = SIN(Z) C PI = FOUR * ZATAN(ONE) IFACA = N / 4 IF ( ITYPE .GT. 0 ) GOTO 5 C C IF THIS IS TO BE AN INVERSE TRANSFORM, CONJUGATE THE DATA C DO 4, K = 1, N XIMAG(K) = -XIMAG(K) 4 CONTINUE 5 IFCAB = IFACA * 4 C C DO THE TRANSFORMS REQUIRED BY THIS STAGE C Z = PI / ZFLOAT(IFCAB) BCOS = -TWO * ZSIN(Z) ** 2 BSIN = ZSIN(TWO * Z) CW1 = ONE SW1 = ZERO DO 10, LITLA = 1, IFACA DO 8, I0 = LITLA, N, IFCAB C C THIS IS THE MAIN CALCULATION OF RADIX 4 TRANSFORMS C I1 = I0 + IFACA I2 = I1 + IFACA I3 = I2 + IFACA XS0 = XREAL(I0) + XREAL(I2) XS1 = XREAL(I0) - XREAL(I2) YS0 = XIMAG(I0) + XIMAG(I2) YS1 = XIMAG(I0) - XIMAG(I2) XS2 = XREAL(I1) + XREAL(I3) XS3 = XREAL(I1) - XREAL(I3) YS2 = XIMAG(I1) + XIMAG(I3) YS3 = XIMAG(I1) - XIMAG(I3) XREAL(I0) = XS0 + XS2 XIMAG(I0) = YS0 + YS2 X1 = XS1 + YS3 Y1 = YS1 - XS3 X2 = XS0 - XS2 Y2 = YS0 - YS2 X3 = XS1 - YS3 Y3 = YS1 + XS3 IF ( LITLA .GT. 1 ) GOTO 7 XREAL(I2) = X1 XIMAG(I2) = Y1 XREAL(I1) = X2 XIMAG(I1) = Y2 XREAL(I3) = X3 XIMAG(I3) = Y3 GOTO 8 C C MULTIPLY BY TWIDDLE FACTORS IF REQUIRED C 7 XREAL(I2) = X1 * CW1 + Y1 * SW1 XIMAG(I2) = Y1 * CW1 - X1 * SW1 XREAL(I1) = X2 * CW2 + Y2 * SW2 XIMAG(I1) = Y2 * CW2 - X2 * SW2 XREAL(I3) = X3 * CW3 + Y3 * SW3 XIMAG(I3) = Y3 * CW3 - X3 * SW3 8 CONTINUE IF ( LITLA .EQ. IFACA ) GOTO 10 C C CALCULATE A NEW SET OF TWIDDLE FACTORS C Z = CW1 * BCOS - SW1 * BSIN + CW1 SW1 = BCOS * SW1 + BSIN * CW1 + SW1 TEMPR = ONE5 - HALF * ( Z * Z + SW1 * SW1 ) CW1 = Z * TEMPR SW1 = SW1 * TEMPR CW2 = CW1 * CW1 - SW1 * SW1 SW2 = TWO * CW1 * SW1 CW3 = CW1 * CW2 - SW1 * SW2 SW3 = CW1 * SW2 + CW2 * SW1 10 CONTINUE IF ( IFACA .LE. 1 ) GOTO 14 C C SET UP THE TRANSFORM SPLIT FOR THE NEXT STAGE C IFACA = IFACA / 4 IF ( IFACA .GT. 0 ) GOTO 5 C C THIS IS THE CALCULATION OF A RADIX TWO STAGE C DO 13, K = 1, N, 2 TEMPR = XREAL(K) + XREAL(K + 1) XREAL(K + 1) = XREAL(K) - XREAL(K + 1) XREAL(K) = TEMPR TEMPR = XIMAG(K) + XIMAG(K + 1) XIMAG(K + 1) = XIMAG(K) - XIMAG(K + 1) XIMAG(K) = TEMPR 13 CONTINUE 14 IF ( ITYPE .GT. 0 ) GOTO 17 C C IF THIS WAS AN INVERSE TRANSFORM, CONJUGATE AND SCALE C THE RESULT C Z = ONE / ZFLOAT(N) DO 16, K = 1, N XIMAG(K) = -XIMAG(K) * Z XREAL(K) = XREAL(K) * Z 16 CONTINUE C 17 RETURN c ... End of subroutine FASTG ... END C SUBROUTINE SCRAM( XREAL, XIMAG, N, IPOW ) C C ALGORITHM AS 83.3 APPL. STATIST. (1975) VOL.24, P.153 C C SUBROUTINE FOR UNSCRAMBLING FFT DATA. C REAL XREAL(*), XIMAG(*), TEMPR INTEGER L(19) EQUIVALENCE (L1, L(1)), (L2, L(2)), (L3, L(3)), $ (L4, L(4)), (L5, L(5)), (L6, L(6)), (L7, L(7)), $ (L8, L(8)), (L9, L(9)), (L10, L(10)), (L11, L(11)), $ (L12, L(12)), (L13, L(13)), (L14, L(14)), (L15, L(15)), $ (L16, L(16)), (L17, L(17)), (L18, L(18)), (L19, L(19)) C II = 1 ITOP = 2 ** ( IPOW - 1 ) I = 20 - IPOW DO 5, K = 1, I L(K) = II 5 CONTINUE L0 = II I = I + 1 DO 6, K = I, 19 II = II * 2 L(K) = II 6 CONTINUE II = 0 DO 9 J1 = 1, L1, L0 DO 9 J2 = J1, L2, L1 DO 9 J3 = J2, L3, L2 DO 9 J4 = J3, L4, L3 DO 9 J5 = J4, L5, L4 DO 9 J6 = J5, L6, L5 DO 9 J7 = J6, L7, L6 DO 9 J8 = J7, L8, L7 DO 9 J9 = J8, L9, L8 DO 9 J10 = J9, L10, L9 DO 9 J11 = J10, L11, L10 DO 9 J12 = J11, L12, L11 DO 9 J13 = J12, L13, L12 DO 9 J14 = J13, L14, L13 DO 9 J15 = J14, L15, L14 DO 9 J16 = J15, L16, L15 DO 9 J17 = J16, L17, L16 DO 9 J18 = J17, L18, L17 DO 9 J19 = J18, L19, L18 J20 = J19 DO 9 I = 1, 2 II = II + 1 IF (II .GE. J20) GOTO 8 C C J20 IS THE BIT-REVERSE OF II C PAIRWISE INTERCHANGE C TEMPR = XREAL(II) XREAL(II) = XREAL(J20) XREAL(J20) = TEMPR TEMPR = XIMAG(II) XIMAG(II) = XIMAG(J20) XIMAG(J20) = TEMPR 8 J20 = J20 + ITOP 9 CONTINUE RETURN c ... End of subroutine SCRAM ... END