FFT Subroutines for Sequences With Special Properties
Return to Main Software Page
C
C-----------------------------------------------------------------------
C MAIN PROGRAM: TEST PROGRAM FOR FFT SUBROUTINES
C AUTHOR: L R RABINER
C BELL LABORATORIES, MURRAY HILL, NEW JERSEY, 07974
C INPUT: RANDOMLY CHOSEN SEQUENCES TO TEST FFT SUBROUTINES
C FOR SEQUENCES WITH SPECIAL PROPERTIES
C N IS THE FFT LENGTH (N MUST BE A POWER OF 2)
C 2<= N <= 4096
C-----------------------------------------------------------------------
C
DIMENSION X(4098), Y(4098)
C
C DEFINE I/0 DEVICE CODES
C INPUT: INPUT TO THIS PROGRAM IS USER-INTERACTIVE
C THAT IS - A QUESTION IS WRITTEN ON THE USER
C TERMINAL (IOUT1) AD THE USER TYPES IN THE ANSWER.
C
C OUTPUT: ALL OUTPUT IS WRITTEN ON THE STANDARD
C OUTPUT UNIT (IOUT2).
C
IND = I1MACH(1)
IOUT1 = I1MACH(4)
IOUT2 = I1MACH(2)
C
C READ IN ANALYSIS SIZE FOR FFT
C
10 WRITE (IOUT1,9999)
9999 FORMAT (30H FFT SIZE(2.LE.N.LE.4096)(I4)=)
READ (IND,9998) N
9998 FORMAT (I4)
IF (N.EQ.0) STOP
DO 20 I=1,12
ITEST = 2**I
IF (N.EQ.ITEST) GO TO 30
20 CONTINUE
WRITE (IOUT1,9997)
9997 FORMAT (45H N IS NOT A POWER OF 2 IN THE RANGE 2 TO 4096)
GO TO 10
30 WRITE (IOUT2,9996) N
9996 FORMAT (11H TESTING N=, I5, 17H RANDOM SEQUENCES)
WRITE (IOUT2,9992)
NP2 = N + 2
NO2 = N/2
NO2P1 = NO2 + 1
NO4 = N/4
NO4P1 = NO4 + 1
C
C CREATE SYMMETRICAL SEQUENCE OF SIZE N
C
DO 40 I=2,NO2
X(I) = UNI(0) - 0.5
IND1 = NP2 - I
X(IND1) = X(I)
40 CONTINUE
X(1) = UNI(0) - 0.5
X(NO2P1) = UNI(0) - 0.5
DO 50 I=1,NO2P1
Y(I) = X(I)
50 CONTINUE
WRITE (IOUT2,9995)
9995 FORMAT (28H ORIGINAL SYMMETRIC SEQUENCE)
WRITE (IOUT2,9993) (X(I),I=1,N)
WRITE (IOUT2,9992)
C
C COMPUTE TRUE FFT OF N POINT SEQUENCE
C
CALL FAST(X, N)
WRITE (IOUT2,9994) N
9994 FORMAT (1H , I4, 32H POINT FFT OF SYMMETRIC SEQUENCE)
WRITE (IOUT2,9993) (X(I),I=1,NP2)
9993 FORMAT (1H , 5E13.5)
WRITE (IOUT2,9992)
9992 FORMAT (1H /1H )
C
C USE SUBROUTINE FFTSYM TO OBTAIN DFT FROM NO2 POINT FFT
C
DO 60 I=1,NO2P1
X(I) = Y(I)
60 CONTINUE
CALL FFTSYM(X, N, Y)
WRITE (IOUT2,9991)
9991 FORMAT (17H OUTPUT OF FFTSYM)
WRITE (IOUT2,9993) (X(I),I=1,NO2P1)
WRITE (IOUT2,9992)
C
C USE SUBROUTINE IFTSYM TO OBTAIN ORIGINAL SEQUENCE FROM NO2 POINT DFT
C
CALL IFTSYM(X, N, Y)
WRITE (IOUT2,9990)
9990 FORMAT (17H OUTPUT OF IFTSYM)
WRITE (IOUT2,9993) (X(I),I=1,NO2P1)
WRITE (IOUT2,9992)
C
C CREATE ANTISYMMETRIC N POINT SEQUENCE
C
DO 70 I=2,NO2
X(I) = UNI(0) - 0.5
IND1 = NP2 - I
X(IND1) = -X(I)
70 CONTINUE
X(1) = 0.
X(NO2P1) = 0.
DO 80 I=1,NO2P1
Y(I) = X(I)
80 CONTINUE
WRITE (IOUT2,9989)
9989 FORMAT (32H ORIGINAL ANTISYMMETRIC SEQUENCE)
WRITE (IOUT2,9993) (X(I),I=1,N)
WRITE (IOUT2,9992)
C
C OBTAIN N POINT DFT OF ANTISYMMETRIC SEQUENCE
C
CALL FAST(X, N)
WRITE (IOUT2,9988) N
9988 FORMAT (1H , I4, 36H POINT FFT OF ANTISYMMETRIC SEQUENCE)
WRITE (IOUT2,9993) (X(I),I=1,NP2)
WRITE (IOUT2,9992)
C
C USE SUBROUTINE FFTASM TO OBTAIN DFT FROM NO2 POINT FFT
C
DO 90 I=1,NO2
X(I) = Y(I)
90 CONTINUE
CALL FFTASM(X, N, Y)
WRITE (IOUT2,9987)
9987 FORMAT (17H OUTPUT OF FFTASM)
WRITE (IOUT2,9993) (X(I),I=1,NO2P1)
WRITE (IOUT2,9992)
C
C USE SUBROUTINE IFTASM TO OBTAIN ORIGINAL SEQUENCE FROM NO2 POINT DFT
C
CALL IFTASM(X, N, Y)
WRITE (IOUT2,9986)
9986 FORMAT (17H OUTPUT OF IFTASM)
WRITE (IOUT2,9993) (X(I),I=1,NO2)
WRITE (IOUT2,9992)
C
C CREATE SEQUENCE WITH ONLY ODD HARMONICS--BEGIN IN FREQUENCY DOMAIN
C
DO 100 I=1,NP2,2
X(I) = 0.
X(I+1) = 0.
IF (MOD(I,4).EQ.1) GO TO 100
X(I) = UNI(0) - 0.5
X(I+1) = UNI(0) - 0.5
IF (N.EQ.2) X(I+1) = 0.
100 CONTINUE
WRITE (IOUT2,9985) N
9985 FORMAT (1H , I4, 35H POINT FFT OF ODD HARMONIC SEQUENCE)
WRITE (IOUT2,9993) (X(I),I=1,NP2)
WRITE (IOUT2,9992)
C
C TRANSFORM BACK TO TIME SEQUENCE
C
CALL FSST(X, N)
WRITE (IOUT2,9984)
9984 FORMAT (31H ORIGINAL ODD HARMONIC SEQUENCE)
WRITE (IOUT2,9993) (X(I),I=1,N)
WRITE (IOUT2,9992)
C
C USE SUBROUTINE FFTOHM TO OBTAIN DFT FROM NO2 POINT FFT
C
CALL FFTOHM(X, N)
WRITE (IOUT2,9983)
9983 FORMAT (17H OUTPUT OF FFTOHM)
WRITE (IOUT2,9993) (X(I),I=1,NO2)
WRITE (IOUT2,9992)
C
C USE SUBROUTINE IFTOHM TO OBTAIN ORIGINAL SEQUENCE FROM NO2 POINT DFT
C
CALL IFTOHM(X, N)
WRITE (IOUT2,9982)
9982 FORMAT (17H OUTPUT OF IFTOHM)
WRITE (IOUT2,9993) (X(I),I=1,NO2)
WRITE (IOUT2,9992)
C
C CREATE SEQUENCE WITH ONLY REAL VALUED ODD HARMONICS
C
DO 110 I=1,NP2,2
X(I) = 0.
X(I+1) = 0.
IF (MOD(I,4).EQ.1) GO TO 110
X(I) = UNI(0) - 0.5
110 CONTINUE
WRITE (IOUT2,9981) N
9981 FORMAT (1H , I4, 45H POINT FFT OF ODD HARMONIC, SYMMETRIC SEQUENC,
* 1HE)
WRITE (IOUT2,9993) (X(I),I=1,NP2)
WRITE (IOUT2,9992)
C
C TRANSFORM BACK TO TIME SEQUENCE
C
CALL FSST(X, N)
WRITE (IOUT2,9980)
9980 FORMAT (42H ORIGINAL ODD HARMONIC, SYMMETRIC SEQUENCE)
WRITE (IOUT2,9993) (X(I),I=1,N)
WRITE (IOUT2,9992)
C
C USE SUBROUTINE FFTSOH TO OBTAIN DFT FROM NO4 POINT FFT
C
CALL FFTSOH(X, N, Y)
WRITE (IOUT2,9979)
9979 FORMAT (17H OUTPUT OF FFTSOH)
WRITE (IOUT2,9993) (X(I),I=1,NO4)
WRITE (IOUT2,9992)
C
C USE SUBROUTINE IFTSOH TO OBTAIN ORIGINAL SEQUENCE FROM NO4 POINT DFT
C
CALL IFTSOH(X, N, Y)
WRITE (IOUT2,9978)
9978 FORMAT (17H OUTPUT OF IFTSOH)
WRITE (IOUT2,9993) (X(I),I=1,NO4)
WRITE (IOUT2,9992)
C
C CREATE SEQUENCE WITH ONLY IMAGINARY VALUED ODD HARMONICS--BEGIN
C IN FREQUENCY DOMAIN
C
DO 120 I=1,NP2,2
X(I) = 0.
X(I+1) = 0.
IF (MOD(I,4).EQ.1) GO TO 120
X(I+1) = UNI(0) - 0.5
120 CONTINUE
WRITE (IOUT2,9977) N
9977 FORMAT (1H , I4, 41H POINT FFT OF ODD HARMONIC, ANTISYMMETRIC,
* 9H SEQUENCE)
WRITE (IOUT2,9993) (X(I),I=1,NP2)
WRITE (IOUT2,9992)
C
C TRANSFORM BACK TO TIME SEQUENCE
C
CALL FSST(X, N)
WRITE (IOUT2,9976)
9976 FORMAT (46H ORIGINAL ODD HARMONIC, ANTISYMMETRIC SEQUENCE)
WRITE (IOUT2,9993) (X(I),I=1,N)
WRITE (IOUT2,9992)
C
C USE SUBROUTINE FFTAOH TO OBTAIN DFT FROM NO4 POINT FFT
C
CALL FFTAOH(X, N, Y)
WRITE (IOUT2,9975)
9975 FORMAT (17H OUTPUT OF FFTAOH)
WRITE (IOUT2,9993) (X(I),I=1,NO4)
WRITE (IOUT2,9992)
C
C USE SUBROUTINE IFTAOH TO OBTAIN ORIGINAL SEQUENCE FROM N/4 POINT DFT
C
CALL IFTAOH(X, N, Y)
WRITE (IOUT2,9974)
9974 FORMAT (17H OUTPUT OF IFTAOH)
WRITE (IOUT2,9993) (X(I),I=1,NO4P1)
WRITE (IOUT2,9992)
C
C BEGIN A NEW PAGE
C
WRITE (IOUT2,9973)
9973 FORMAT (1H1)
GO TO 10
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FFTSYM
C COMPUTE DFT FOR REAL, SYMMETRIC, N-POINT SEQUENCE X(M) USING
C N/2-POINT FFT
C SYMMETRIC SEQUENCE MEANS X(M)=X(N-M), M=1,...,N/2-1
C NOTE: INDEX M IS SEQUENCE INDEX--NOT FORTRAN INDEX
C-----------------------------------------------------------------------
C
SUBROUTINE FFTSYM(X, N, Y)
DIMENSION X(1), Y(1)
C
C X = REAL ARRAY WHICH ON INPUT CONTAINS THE N/2+1 POINTS OF THE
C INPUT SEQUENCE (SYMMETRICAL)
C ON OUTPUT X CONTAINS THE N/2+1 REAL POINTS OF THE TRANSFORM OF
C THE INPUT--I.E. THE ZERO VALUED IMAGINARY PARTS ARE NOT RETURNED
C N = TRUE SIZE OF INPUT
C Y = SCRATCH ARRAY OF SIZE N/2+2
C
C
C FOR N = 2, COMPUTE DFT DIRECTLY
C
IF (N.GT.2) GO TO 10
T = X(1) + X(2)
X(2) = X(1) - X(2)
X(1) = T
RETURN
10 TWOPI = 8.*ATAN(1.0)
C
C FIRST COMPUTE B0 TERM, WHERE B0=SUM OF ODD VALUES OF X(M)
C
NO2 = N/2
NO4 = N/4
NIND = NO2 + 1
B0 = 0.
DO 20 I=2,NIND,2
B0 = B0 + X(I)
20 CONTINUE
B0 = B0*2.
C
C FOR N = 4 SKIP RECURSION LOOP
C
IF (N.EQ.4) GO TO 40
C
C FORM NEW SEQUENCE, Y(M)=X(2*M)+(X(2*M+1)-X(2*M-1))
C
DO 30 I=2,NO4
IND = 2*I
T1 = X(IND) - X(IND-2)
Y(I) = X(IND-1) + T1
IND1 = NO2 + 2 - I
Y(IND1) = X(IND-1) - T1
30 CONTINUE
40 Y(1) = X(1)
Y(NO4+1) = X(NO2+1)
C
C TAKE N/2 POINT (REAL) FFT OF Y
C
CALL FAST(Y, NO2)
C
C FORM ORIGINAL DFT BY UNSCRAMBLING Y(K)
C USE RECURSION TO GIVE SIN(TPN*I) MULTIPLIER
C
TPN = TWOPI/FLOAT(N)
COSI = 2.*COS(TPN)
SINI = 2.*SIN(TPN)
COSD = COSI/2.
SIND = SINI/2.
NIND = NO4 + 1
DO 50 I=2,NIND
IND = 2*I
BK = Y(IND)/SINI
AK = Y(IND-1)
X(I) = AK + BK
NIND1 = N/2 + 2 - I
X(NIND1) = AK - BK
TEMP = COSI*COSD - SINI*SIND
SINI = COSI*SIND + SINI*COSD
COSI = TEMP
50 CONTINUE
X(1) = B0 + Y(1)
X(NO2+1) = Y(1) - B0
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: IFTSYM
C COMPUTE IDFT FOR REAL, SYMMETRIC, N-POINT SEQUENCE X(M) USING
C N/2-POINT FFT
C SYMMETRIC SEQUENCE MEANS X(M)=X(N-M), M=1,...,N/2-1
C NOTE: INDEX M IS SEQUENCE INDEX--NOT FORTRAN INDEX
C-----------------------------------------------------------------------
C
SUBROUTINE IFTSYM(X, N, Y)
DIMENSION X(1), Y(1)
C
C X = REAL ARRAY WHICH ON INPUT CONTAINS THE N/2+1 REAL POINTS OF THE
C TRANSFORM OF THE INPUT--I.E. THE ZERO VALUED IMAGINARY PARTS
C ARE NOT GIVEN AS INPUT
C ON OUTPUT X CONTAINS THE N/2+1 POINTS OF THE TIME SEQUENCE
C (SYMMETRICAL)
C N = TRUE SIZE OF INPUT
C Y = SCRATCH ARRAY OF SIZE N/2+2
C
C
C FOR N = 2, COMPUTE IDFT DIRECTLY
C
IF (N.GT.2) GO TO 10
T = (X(1)+X(2))/2.
X(2) = (X(1)-X(2))/2.
X(1) = T
RETURN
10 TWOPI = 8.*ATAN(1.0)
C
C FIRST COMPUTE X1=X(1) TERM DIRECTLY
C USE RECURSION ON THE SINE COSINE TERMS
C
NO2 = N/2
NO4 = N/4
TPN = TWOPI/FLOAT(N)
COSD = COS(TPN)
SIND = SIN(TPN)
COSI = 2.
SINI = 0.
X1 = X(1) - X(NO2+1)
DO 20 I=2,NO2
TEMP = COSI*COSD - SINI*SIND
SINI = COSI*SIND + SINI*COSD
COSI = TEMP
X1 = X1 + X(I)*COSI
20 CONTINUE
X1 = X1/FLOAT(N)
C
C SCRAMBLE ORIGINAL DFT (X(K)) TO GIVE Y(K)
C USE RECURSION RELATION TO GENERATE SIN(TPN*I) MULTIPLIER
C
COSI = COS(TPN)
SINI = SIN(TPN)
COSD = COSI
SIND = SINI
Y(1) = (X(1)+X(NO2+1))/2.
Y(2) = 0.
NIND = NO4 + 1
DO 30 I=2,NIND
IND = 2*I
NIND1 = NO2 + 2 - I
AK = (X(I)+X(NIND1))/2.
BK = (X(I)-X(NIND1))
Y(IND-1) = AK
Y(IND) = BK*SINI
TEMP = COSI*COSD - SINI*SIND
SINI = COSI*SIND + SINI*COSD
COSI = TEMP
30 CONTINUE
C
C TAKE N/2 POINT IDFT OF Y
C
CALL FSST(Y, NO2)
C
C FORM X SEQUENCE FROM Y SEQUENCE
C
X(1) = Y(1)
X(2) = X1
IF (N.EQ.4) GO TO 50
DO 40 I=2,NO4
IND = 2*I
IND1 = NO2 + 2 - I
X(IND-1) = (Y(I)+Y(IND1))/2.
T1 = (Y(I)-Y(IND1))/2.
X(IND) = T1 + X(IND-2)
40 CONTINUE
50 X(NO2+1) = Y(NO4+1)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FFTASM
C COMPUTE DFT FOR REAL, ANTISYMMETRIC, N-POINT SEQUENCE X(M) USING
C N/2-POINT FFT
C ANTISYMMETRIC SEQUENCE MEANS X(M)=-X(N-M), M=1,...,N/2-1
C NOTE: INDEX M IS SEQUENCE INDEX--NOT FORTRAN INDEX
C-----------------------------------------------------------------------
C
SUBROUTINE FFTASM(X, N, Y)
DIMENSION X(1), Y(1)
C
C X = REAL ARRAY WHICH ON INPUT CONTAINS THE N/2 POINTS OF THE
C INPUT SEQUENCE (ASYMMETRICAL)
C ON OUTPUT X CONTAINS THE N/2+1 IMAGINARY POINTS OF THE TRANSFORM
C OF THE INPUT--I.E. THE ZERO VALUED REAL PARTS ARE NOT RETURNED
C N = TRUE SIZE OF INPUT
C Y = SCRATCH ARRAY OF SIZE N/2+2
C
C
C FOR N = 2, ASSUME X(1)=0, X(2)=0, COMPUTE DFT DIRECTLY
C
IF (N.EQ.2) GO TO 30
TWOPI = 8.*ATAN(1.0)
C
C FORM NEW SEQUENCE, Y(M)=X(2*M)+(X(2*M+1)-X(2*M-1))
C
NO2 = N/2
NO4 = N/4
DO 10 I=2,NO4
IND = 2*I
T1 = X(IND) - X(IND-2)
Y(I) = X(IND-1) + T1
IND1 = NO2 + 2 - I
Y(IND1) = -X(IND-1) + T1
10 CONTINUE
Y(1) = 2.*X(2)
Y(NO4+1) = -2.*X(NO2)
C
C TAKE N/2 POINT (REAL) FFT OF Y
C
CALL FAST(Y, NO2)
C
C FORM ORIGINAL DFT BY UNSCRAMBLING Y(K)
C USE RECURSION RELATION TO GENERATE SIN(TPN*I) MULTIPLIER
C
TPN = TWOPI/FLOAT(N)
COSI = 2.*COS(TPN)
SINI = 2.*SIN(TPN)
COSD = COSI/2.
SIND = SINI/2.
NIND = NO4 + 1
DO 20 I=2,NIND
IND = 2*I
BK = Y(IND-1)/SINI
AK = Y(IND)
X(I) = AK - BK
IND1 = NO2 + 2 - I
X(IND1) = -AK - BK
TEMP = COSI*COSD - SINI*SIND
SINI = COSI*SIND + SINI*COSD
COSI = TEMP
20 CONTINUE
30 X(1) = 0.
X(NO2+1) = 0.
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: IFTASM
C COMPUTE IDFT FOR REAL, ANTISYMMETRIC, N-POINT SEQUENCE X(M) USING
C N/2-POINT FFT
C ANTISYMMETRIC SEQUENCE MEANS X(M)=-X(N-M), M=1,...,N/2-1
C NOTE: INDEX M IS SEQUENCE INDEX--NOT FORTRAN INDEX
C-----------------------------------------------------------------------
C
SUBROUTINE IFTASM(X, N, Y)
DIMENSION X(1), Y(1)
C
C X = IMAGINARY ARRAY WHICH ON INPUT CONTAINS THE N/2+1 REAL POINTS OF
C THE TRANSFORM OF THE INPUT--I.E. THE ZERO VALUED REAL PARTS
C ARE NOT GIVEN AS INPUT
C ON OUTPUT X CONTAINS THE N/2 POINTS OF THE TIME SEQUENCE
C (ANTISYMMETRICAL)
C N = TRUE SIZE OF INPUT
C Y = SCRATCH ARRAY OF SIZE N/2+2
C
C
C FOR N = 2, ASSUME X(1)=0, X(2)=0
C
IF (N.GT.2) GO TO 10
X(1) = 0
X(2) = 0
RETURN
10 TWOPI = 8.*ATAN(1.0)
C
C FIRST COMPUTE X1=X(1) TERM DIRECTLY
C USE RECURSION ON THE SINE COSINE TERMS
C
NO2 = N/2
NO4 = N/4
TPN = TWOPI/FLOAT(N)
C
C SCRAMBLE ORIGINAL DFT (X(K)) TO GIVE Y(K)
C USE RECURSION RELATION TO GIVE SIN(TPN*I) MULTIPLIER
C
COSI = COS(TPN)
SINI = SIN(TPN)
COSD = COSI
SIND = SINI
NIND = NO4 + 1
DO 20 I=2,NIND
IND = 2*I
IND1 = NO2 + 2 - I
AK = (X(I)-X(IND1))/2.
BK = -(X(I)+X(IND1))
Y(IND) = AK
Y(IND-1) = BK*SINI
TEMP = COSI*COSD - SINI*SIND
SINI = COSI*SIND + SINI*COSD
COSI = TEMP
20 CONTINUE
Y(1) = 0.
Y(2) = 0.
C
C TAKE N/2 POINT IDFT OF Y
C
CALL FSST(Y, NO2)
C
C FORM X SEQUENCE FROM Y SEQUENCE
C
X(2) = Y(1)/2.
X(1) = 0.
IF (N.EQ.4) GO TO 40
DO 30 I=2,NO4
IND = 2*I
IND1 = NO2 + 2 - I
X(IND-1) = (Y(I)-Y(IND1))/2.
T1 = (Y(I)+Y(IND1))/2.
X(IND) = T1 + X(IND-2)
30 CONTINUE
40 X(NO2) = -Y(NO4+1)/2.
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FFTOHM
C COMPUTE DFT FOR REAL, N-POINT, ODD HARMONIC SEQUENCES USING AN
C N/2 POINT FFT
C ODD HARMONIC MEANS X(2*K)=0, ALL K WHERE X(K) IS THE DFT OF X(M)
C NOTE: INDEX M IS SEQUENCE INDEX--NOT FORTRAN INDEX
C-----------------------------------------------------------------------
C
SUBROUTINE FFTOHM(X, N)
DIMENSION X(1)
C
C X = REAL ARRAY WHICH ON INPUT CONTAINS THE FIRST N/2 POINTS OF THE
C INPUT
C ON OUTPUT X CONTAINS THE N/4 COMPLEX VALUES OF THE ODD
C HARMONICS OF THE INPUT--STORED IN THE SEQUENCE RE(X(1)),IM(X(1)),
C RE(X(2)),IM(X(2)),...
C ****NOTE: X MUST BE DIMENSIONED TO SIZE N/2+2 FOR FFT ROUTINE
C N = TRUE SIZE OF X SEQUENCE
C
C FIRST COMPUTE REAL(X(1)) AND REAL(X(N/2-1)) SEPARATELY
C ALSO SIMULTANEOUSLY MULTIPLY ORIGINAL SEQUENCE BY SIN(TWOPI*(M-1)/N)
C SIN AND COS ARE COMPUTED RECURSIVELY
C
C
C FOR N = 2, ASSUME X(1)=X0, X(2)=-X0, COMPUTE DFT DIRECTLY
C
IF (N.GT.2) GO TO 10
X(1) = 2.*X(1)
X(2) = 0.
RETURN
10 TWOPI = 8.*ATAN(1.0)
TPN = TWOPI/FLOAT(N)
C
C COMPUTE X1=REAL(X(1)) AND X2=IMAGINARY(X(N/2-1))
C X(N) = X(N)*4.*SIN(TWOPI*(I-1)/N)
C
T1 = 0.
C
C COSD AND SIND ARE MULTIPLIERS FOR RECURSION FOR SIN AND COS
C COSI AND SINI ARE INITIAL CONDITIONS FOR RECURSION FOR SIN AND COS
C
COSD = COS(TPN*2.)
SIND = SIN(TPN*2.)
COSI = 1.
SINI = 0.
NO2 = N/2
DO 20 I=1,NO2,2
T = X(I)*COSI
X(I) = X(I)*4.*SINI
TEMP = COSI*COSD - SINI*SIND
SINI = COSI*SIND + SINI*COSD
COSI = TEMP
T1 = T1 + T
20 CONTINUE
C
C RESET INITIAL CONDITIONS (COSI,SINI) FOR NEW RECURSION
C
COSI = COS(TPN)
SINI = SIN(TPN)
T2 = 0.
DO 30 I=2,NO2,2
T = X(I)*COSI
X(I) = X(I)*4.*SINI
TEMP = COSI*COSD - SINI*SIND
SINI = COSI*SIND + SINI*COSD
COSI = TEMP
T2 = T2 + T
30 CONTINUE
X1 = 2.*(T1+T2)
X2 = 2.*(T1-T2)
C
C TAKE N/2 POINT (REAL) FFT OF PREPROCESSED SEQUENCE X
C
CALL FAST(X, NO2)
C
C FOR N = 4--SKIP RECURSION AND INITIAL CONDITIONS
C
IF (N.EQ.4) GO TO 50
C
C INITIAL CONDITIONS FOR RECURSION
C
X(2) = -X(1)/2.
X(1) = X1
C
C FOR N = 8, SKIP RECURSION
C
IF (N.EQ.8) GO TO 50
C
C UNSCRAMBLE Y(K) USING RECURSION FORMULA
C
NIND = NO2 - 2
DO 40 I=3,NIND,2
T = X(I)
X(I) = X(I-2) + X(I+1)
X(I+1) = X(I-1) - T
40 CONTINUE
50 X(NO2) = X(NO2+1)/2.
X(NO2-1) = X2
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: IFTOHM
C COMPUTE IDFT FOR REAL, N-POINT, ODD HARMONIC SEQUENCES USING AN
C N/2 POINT FFT
C ODD HARMONIC MEANS X(2*K)=0, ALL K WHERE X(K) IS THE DFT OF X(M)
C NOTE: INDEX M IS SEQUENCE INDEX--NOT FORTRAN INDEX
C-----------------------------------------------------------------------
C
SUBROUTINE IFTOHM(X, N)
DIMENSION X(1)
C
C X = REAL ARRAY WHICH ON INPUT CONTAINS THE N/4 COMPLEX VALUES OF THE
C ODD HARMONICS OF THE INPUT--STORED IN THE SEQUENCE RE(X(1)),
C IM(X(1)),RE(X(2)),IM(X(2)),...
C ON OUTPUT X CONTAINS THE FIRST N/2 POINTS OF THE INPUT
C ****NOTE: X MUST BE DIMENSIONED TO SIZE N/2+2 FOR FFT ROUTINE
C N = TRUE SIZE OF X SEQUENCE
C
C FIRST COMPUTE REAL(X(1)) AND REAL(X(N/2-1)) SEPARATELY
C ALSO SIMULTANEOUSLY MULTIPLY ORIGINAL SEQUENCE BY SIN(TWOPI*(M-1)/N)
C SIN AND COS ARE COMPUTED RECURSIVELY
C
C
C FOR N = 2, ASSUME X(1)=X0, X(2)=-X0, COMPUTE IDFT DIRECTLY
C
IF (N.GT.2) GO TO 10
X(1) = 0.5*X(1)
X(2) = -X(1)
RETURN
10 TWOPI = 8.*ATAN(1.0)
TPN = TWOPI/FLOAT(N)
NO2 = N/2
NO4 = N/4
NIND = NO2
C
C SOLVE FOR X(0)=X0 DIRECTLY
C
X0 = 0.
DO 20 I=1,NO2,2
X0 = X0 + 2.*X(I)
20 CONTINUE
X0 = X0/FLOAT(N)
C
C FORM Y(K)=J*(X(2K+1)-X(2K-1))
C OVERWRITE X ARRAY WITH Y SEQUENCE
C
XPR = X(1)
XPI = X(2)
X(1) = -2.*X(2)
X(2) = 0.
IF (NO4.EQ.1) GO TO 40
DO 30 I=3,NIND,2
TI = X(I) - XPR
TR = -X(I+1) + XPI
XPR = X(I)
XPI = X(I+1)
X(I) = TR
X(I+1) = TI
30 CONTINUE
40 X(NO2+1) = 2.*XPI
X(NO2+2) = 0.
C
C TAKE N/2 POINT (REAL) IFFT OF PREPROCESSED SEQUENCE X
C
CALL FSST(X, NO2)
C
C SOLVE FOR X(M) BY DIVIDING BY 4*SIN(TWOPI*M/N) FOR M=1,2,...,N/2-1
C FOR M=0 SUBSTITUTE PRECOMPUTED VALUE X0
C
COSI = 4.
SINI = 0.
COSD = COS(TPN)
SIND = SIN(TPN)
DO 50 I=2,NO2
TEMP = COSI*COSD - SINI*SIND
SINI = COSI*SIND + SINI*COSD
COSI = TEMP
X(I) = X(I)/SINI
50 CONTINUE
X(1) = X0
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FFTSOH
C COMPUTE DFT FOR REAL, SYMMETRIC, ODD HARMONIC, N-POINT SEQUENCE
C USING N/4-POINT FFT
C SYMMETRIC SEQUENCE MEANS X(M)=X(N-M), M=1,...,N/2-1
C ODD HARMONIC MEANS X(2*K)=0, ALL K, WHERE X(K) IS THE DFT OF X(M)
C X(M) HAS THE PROPERTY X(M)=-X(N/2-M), M=0,1,...,N/4-1, X(N/4)=0
C NOTE: INDEX M IS SEQUENCE INDEX--NOT FORTRAN INDEX
C-----------------------------------------------------------------------
C
SUBROUTINE FFTSOH(X, N, Y)
DIMENSION X(1), Y(1)
C
C X = REAL ARRAY WHICH ON INPUT CONTAINS THE N/4 POINTS OF THE
C INPUT SEQUENCE (SYMMETRICAL)
C ON OUTPUT X CONTAINS THE N/4 REAL POINTS OF THE ODD HARMONICS
C OF THE TRANSFORM OF THE INPUT--I.E. THE ZERO VALUED IMAGINARY
C PARTS ARE NOT GIVEN NOR ARE THE ZERO-VALUED EVEN HARMONICS
C N = TRUE SIZE OF INPUT
C Y = SCRATCH ARRAY OF SIZE N/4+2
C
C
C HANDLE N = 2 AND N = 4 CASES SEPARATELY
C
IF (N.GT.4) GO TO 20
IF (N.EQ.4) GO TO 10
C
C FOR N=2, ASSUME X(1)=X0, X(2)=-X0, COMPUTE DFT DIRECTLY
C
X(1) = 2.*X(1)
RETURN
C
C N = 4 CASE, COMPUTE DFT DIRECTLY
C
10 X(1) = 2.*X(1)
RETURN
20 TWOPI = 8.*ATAN(1.0)
C
C FORM NEW SEQUENCE, Y(M)=X(2*M)+(X(2*M+1)-X(2*M-1))
C
NO2 = N/2
NO4 = N/4
NO8 = N/8
IF (NO8.EQ.1) GO TO 40
DO 30 I=2,NO8
IND = 2*I
T1 = X(IND) - X(IND-2)
Y(I) = X(IND-1) + T1
IND1 = N/4 + 2 - I
Y(IND1) = -X(IND-1) + T1
30 CONTINUE
40 Y(1) = X(1)
Y(NO8+1) = -2.*X(NO4)
C
C THE SEQUENCE Y (N/4 POINTS) HAS ONLY ODD HARMONICS
C CALL SUBROUTINE FFTOHM TO EXPLOIT ODD HARMONICS
C
CALL FFTOHM(Y, NO2)
C
C FORM ORIGINAL DFT FROM COMPLEX ODD HARMONICS OF Y(K)
C BY UNSCRAMBLING Y(K)
C
TPN = TWOPI/FLOAT(N)
COSI = 2.*COS(TPN)
SINI = 2.*SIN(TPN)
COSD = COS(TPN*2.)
SIND = SIN(TPN*2.)
DO 50 I=1,NO8
IND = 2*I
BK = Y(IND)/SINI
TEMP = COSI*COSD - SINI*SIND
SINI = COSI*SIND + SINI*COSD
COSI = TEMP
AK = Y(IND-1)
X(I) = AK + BK
IND1 = N/4 + 1 - I
X(IND1) = AK - BK
50 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: IFTSOH
C COMPUTE IDFT FOR REAL, SYMMETRIC, ODD HARMONIC, N-POINT SEQUENCE
C USING N/4-POINT FFT
C SYMMETRIC SEQUENCE MEANS X(M)=X(N-M), M=1,...,N/2-1
C ODD HARMONIC MEANS X(2*K)=0, ALL K, WHERE X(K) IS THE DFT OF X(M)
C X(M) HAS THE PROPERTY X(M)=-X(N/2-M), M=0,1,...,N/4-1, X(N/4)=0
C NOTE: INDEX M IS SEQUENCE INDEX--NOT FORTRAN INDEX
C-----------------------------------------------------------------------
C
SUBROUTINE IFTSOH(X, N, Y)
DIMENSION X(1), Y(1)
C
C X = REAL ARRAY WHICH ON INPUT CONTAINS THE N/4 REAL POINTS OF
C THE ODD HARMONICS OF THE TRANSFORM OF THE ORIGINAL TIME SEQUENCE
C I.E. THE ZERO VALUED IMAGINARY PARTS ARE NOT GIVEN NOR ARE THE
C ZERO VALUED EVEN HARMONICS
C ON OUTPUT X CONTAINS THE FIRST N/4 POINTS OF THE ORIGINAL INPUT
C SEQUENCE (SYMMETRICAL)
C N = TRUE SIZE OF INPUT
C Y = SCRATCH ARRAY OF SIZE N/4+2
C
C
C HANDLE N = 2 AND N = 4 CASES SEPARATELY
C
IF (N.GT.4) GO TO 10
C
C FOR N=2, 4 ASSUME X(1)=X0, X(2)=-X0, COMPUTE IDFT DIRECTLY
C
X(1) = X(1)/2.
RETURN
C
C CODE FOR VALUES OF N WHICH ARE MULTIPLES OF 8
C
10 TWOPI = 8.*ATAN(1.0)
NO2 = N/2
NO4 = N/4
NO8 = N/8
TPN = TWOPI/FLOAT(N)
C
C FIRST COMPUTE X1=X(1) TERM DIRECTLY
C USE RECURSION ON THE SINE COSINE TERMS
C
COSD = COS(TPN*2.)
SIND = SIN(TPN*2.)
COSI = 2.*COS(TPN)
SINI = 2.*SIN(TPN)
X1 = 0.
DO 20 I=1,NO4
X1 = X1 + X(I)*COSI
TEMP = COSI*COSD - SINI*SIND
SINI = COSI*SIND + SINI*COSD
COSI = TEMP
20 CONTINUE
X1 = X1/FLOAT(N)
C
C SCRAMBLE ORIGINAL DFT (X(K)) TO GIVE Y(K)
C USE RECURSION RELATION TO GIVE SIN MULTIPLIERS
C
COSI = COS(TPN)
SINI = SIN(TPN)
DO 30 I=1,NO8
IND = 2*I
IND1 = NO4 + 1 - I
AK = (X(I)+X(IND1))/2.
BK = (X(I)-X(IND1))
Y(IND-1) = AK
Y(IND) = BK*SINI
TEMP = COSI*COSD - SINI*SIND
SINI = COSI*SIND + SINI*COSD
COSI = TEMP
30 CONTINUE
C
C THE SEQUENCE Y(K) IS THE ODD HARMONICS DFT OUTPUT
C USE SUBROUTINE IFTOHM TO OBTAIN Y(M), THE INVERSE TRANSFORM
C
CALL IFTOHM(Y, NO2)
C
C FORM X(M) SEQUENCE FROM Y(M) SEQUENCE
C USE X1 INITIAL CONDITION ON THE RECURSION
C
X(1) = Y(1)
X(2) = X1
IF (NO8.EQ.1) RETURN
DO 40 I=2,NO8
IND = 2*I
IND1 = NO4 + 2 - I
T1 = (Y(I)+Y(IND1))/2.
X(IND-1) = (Y(I)-Y(IND1))/2.
X(IND) = T1 + X(IND-2)
40 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FFTAOH
C COMPUTE DFT FOR REAL, ANTISYMMETRIC, ODD HARMONIC, N-POINT SEQUENCE
C USING N/4-POINT FFT
C ANTISYMMETRIC SEQUENCE MEANS X(M)=-X(N-M), M=1,...,N/2-1
C ODD HARMONIC MEANS X(2*K)=0, ALL K, WHERE X(K) IS THE DFT OF X(M)
C X(M) HAS THE PROPERTY X(M)=X(N/2-M), M=0,1,...,N/4-1, X(0)=0
C NOTE: INDEX M IS SEQUENCE INDEX--NOT FORTRAN INDEX
C-----------------------------------------------------------------------
C
SUBROUTINE FFTAOH(X, N, Y)
DIMENSION X(1), Y(1)
C
C X = REAL ARRAY WHICH ON INPUT CONTAINS THE (N/4+1) POINTS OF THE
C INPUT SEQUENCE (ANTISYMMETRICAL)
C ON OUTPUT X CONTAINS THE N/4 IMAGINARY POINTS OF THE ODD
C HARMONICS OF THE TRANSFORM OF THE INPUT--I.E. THE ZERO
C VALUED REAL PARTS ARE NOT GIVEN NOR ARE THE ZERO-VALUED
C EVEN HARMONICS
C N = TRUE SIZE OF INPUT
C Y = SCRATCH ARRAY OF SIZE N/4+2
C
C
C HANDLE N = 2 AND N = 4 CASES SEPARATELY
C
IF (N.GT.4) GO TO 20
IF (N.EQ.4) GO TO 10
C
C FOR N=2, ASSUME X(1)=0, X(2)=0, COMPUTE DFT DIRECTLY
C
X(1) = 0.
RETURN
C
C N = 4 CASE, ASSUME X(1)=X(3)=0, X(2)=-X(4)=X0, COMPUTE DFT DIRECTLY
C
10 X(1) = -2.*X(2)
RETURN
20 TWOPI = 8.*ATAN(1.0)
C
C FORM NEW SEQUENCE, Y(M)=X(2*M)+(X(2*M+1)-X(2*M-1))
C
NO2 = N/2
NO4 = N/4
NO8 = N/8
IF (NO8.EQ.1) GO TO 40
DO 30 I=2,NO8
IND = 2*I
T1 = X(IND) - X(IND-2)
Y(I) = X(IND-1) + T1
IND1 = N/4 + 2 - I
Y(IND1) = X(IND-1) - T1
30 CONTINUE
40 Y(1) = 2.*X(2)
Y(NO8+1) = X(NO4+1)
C
C THE SEQUENCE Y (N/4 POINTS) HAS ONLY ODD HARMONICS
C CALL SUBROUTINE FFTOHM TO EXPLOIT ODD HARMONICS
C
CALL FFTOHM(Y, NO2)
C
C FORM ORIGINAL DFT FROM COMPLEX ODD HARMONICS OF Y(K)
C BY UNSCRAMBLING Y(K)
C
TPN = TWOPI/FLOAT(N)
COSI = 2.*COS(TPN)
SINI = 2.*SIN(TPN)
COSD = COS(TPN*2.)
SIND = SIN(TPN*2.)
DO 50 I=1,NO8
IND = 2*I
BK = Y(IND-1)/SINI
TEMP = COSI*COSD - SINI*SIND
SINI = COSI*SIND + SINI*COSD
COSI = TEMP
AK = Y(IND)
X(I) = AK - BK
IND1 = N/4 + 1 - I
X(IND1) = -AK - BK
50 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: IFTAOH
C COMPUTE IDFT FOR REAL, ANTISYMMETRIC, ODD HARMONIC, N-POINT SEQUENCE
C USING N/4-POINT FFT
C ANTISYMMETRIC SEQUENCE MEANS X(M)=-X(N-M), M=1,...,N/2-1
C ODD HARMONIC MEANS X(2*K)=0, ALL K, WHERE X(K) IS THE DFT OF X(M)
C X(M)HAS THE PROPERTY X(M)=X(N/2-M), M=0,1,...,N/4-1, X(0)=0
C NOTE: INDEX M IS SEQUENCE INDEX--NOT FORTRAN INDEX
C-----------------------------------------------------------------------
C
SUBROUTINE IFTAOH(X, N, Y)
DIMENSION X(1), Y(1)
C
C X = REAL ARRAY WHICH ON INPUT CONTAINS THE N/4 IMAGINARY POINTS
C OF THE ODD HARMONICS OF THE TRANSFORM OF THE ORIGINAL TIME
C SEQUENCE--I.E. THE ZERO VALUED REAL PARTS ARE NOT INPUT NOR
C ARE THE ZERO-VALUED EVEN HARMONICS
C ON OUTPUT X CONTAINS THE FIRST (N/4+1) POINTS OF THE ORIGINAL
C TIME SEQUENCE (ANTISYMMETRICAL)
C N = TRUE SIZE OF INPUT
C Y = SCRATCH ARRAY OF SIZE N/4+2
C
C
C HANDLE N = 2 AND N = 4 CASES SEPARATELY
C
IF (N.GT.4) GO TO 20
IF (N.EQ.4) GO TO 10
C
C FOR N=2 ASSUME X(1)=0, X(2)=0, COMPUTE IDFT DIRECTLY
C
X(1) = 0.
RETURN
C
C FOR N=4, ASSUME X(1)=X(3)=0, X(2)=-X(4)=X0, COMPUTE IDFT DIRECTLY
C
10 X(2) = -X(1)/2.
X(1) = 0.
RETURN
C
C CODE FOR VALUES OF N WHICH ARE MULTIPLES OF 8
C
20 TWOPI = 8.*ATAN(1.0)
NO2 = N/2
NO4 = N/4
NO8 = N/8
TPN = TWOPI/FLOAT(N)
C
C SCRAMBLE ORIGINAL DFT (X(K)) TO GIVE Y(K)
C USE RECURSION TO GIVE SIN MULTIPLIERS
C
COSI = COS(TPN)
SINI = SIN(TPN)
COSD = COS(TPN*2.)
SIND = SIN(TPN*2.)
DO 30 I=1,NO8
IND = 2*I
IND1 = NO4 + 1 - I
AK = (X(I)-X(IND1))/2.
BK = -(X(I)+X(IND1))
Y(IND) = AK
Y(IND-1) = BK*SINI
TEMP = COSI*COSD - SINI*SIND
SINI = COSI*SIND + SINI*COSD
COSI = TEMP
30 CONTINUE
C
C THE SEQUENCE Y(K) IS AN ODD HARMONIC SEQUENCE
C USE SUBROUTINE IFTOHM TO GIVE Y(M)
C
CALL IFTOHM(Y, NO2)
C
C FORM X SEQUENCE FROM Y SEQUENCE
C
X(2) = Y(1)/2.
X(1) = 0.
IF (N.EQ.8) RETURN
DO 40 I=2,NO8
IND = 2*I
IND1 = NO4 + 2 - I
X(IND-1) = (Y(I)+Y(IND1))/2.
T1 = (Y(I)-Y(IND1))/2.
X(IND) = T1 + X(IND-2)
40 CONTINUE
X(NO4+1) = Y(NO8+1)
RETURN
END
Return to Main Software Page