Make your own free website on Tripod.com
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