C    -------------------------------------------------------------------------
C    |  File name: fft.for               Sean O'Connor               9-23-87 |
C    ------------------------------------------------------------------------- 

C    -------------------------------------------------------------------------
C    |  VAX 11/785                 VMS 4.5         VAX EXTENDED FORTRAN 77   |
C    -------------------------------------------------------------------------

C 
C  ****************************************************************************
C  *                                                                          *
C  *                            FFT                                           *
C  *                                                                          *
C  ****************************************************************************
C
C
C  DESCRIPTION
C
C      This program does an in-place one dimensional complex discrete Fourier
C      transform.  It uses an FFT algorithm for a number of input points which
C      is a power of two.
C
C
C  INPUT 
C
C      N,  The number of complex points to be transformed.
C          N is a power of 2.  N must be at least 2;  it
C          must not exceed the dimension of array X, or the 
C          maximum size of the sine and cosine table.  (See 
C          below).
C
C      X,  The array holding the complex data X( k )
C          for k = 0, 1, . . ., N - 1 which is to be
C          transformed.  It is dimensioned in the calling
C          program.  Only the array address is passed to 
C          this subroutine.
C
C      D,  The direction of the FFT. D is 1 for a 
C          forward transform and -1 for an inverse
C          transform.
C
C
C  OUTPUT
C
C      X,  On output, X holds the complex Fourier coefficients
C          ^
C          X( k ) for k = 0, 1, . . ., N - 1.
C
C      Warning messages to the terminal under these conditions:
C
C          (1)  D is not 1 or -1
C          (2)  N is less than 2.
C          (3)  N is greater than the sine table can handle.
C          (4)  N is not a power of 2.
C
C
C  CALLING SEQUENCE
C
C         CALL FFT ( N, X, D )
C
C
C  EXAMPLE 
C
C      DIMENSION X( 0 : 3 )        ! 4 point FFT
C
C      N = 4                       ! 4 point FFT
C      D = 1                       ! Forward transform
C
C      CALL FFT( N, X, D )         ! Do the FFT
C
C      D = -1                      ! Do an inverse FFT
C
C      CALL FFT( N, X, D )
C
C   Suppose X contains the data
C
C   X ( 0 )  =  1  +  i
C   X ( 1 )  =  2  +  2 i
C   X ( 2 )  =  3  +  3 i
C   X ( 3 )  =  4  +  4 i
C 
C   then the forward FFT returns
C
C   X ( 0 ) =   5  +  5 i
C   X ( 1 ) =  -2
C   X ( 2 ) =  -1  -    i
C   X ( 3 ) =      -  2 i
C
C   and the inverse FFT returns the original data.
C
C
C  METHOD
C                                ^
C       The Fourier coefficients X( k ) are defined for k = 0, . . ., N - 1
C       by the formula
C
C                             N - 1               2 Pi i
C       ^             1       ____             -  ------ j k
C       X( k ) =  ---------   \                     N
C                   -----     /      X( j )  e
C                \ /  N       ----
C                 v           j = 0
C
C       where e is 2.1718 . . ., Pi = 3.14159. . . and i is the square root
C       of -1.
C
C       This is the formula for the forward transform.  The inverse transform
C       is the same, except that the minus sign in the exponent above is a
C       plus sign instead.
C                                    ---
C       The normalizing factor 1 / \/ N  in front was picked so that the
C       forward and inverse transforms look the same.
C
C
C  The particular type of FFT algorithm we use is a
C
C                             complex,
C                           two to the N,
C                     decimation in  frequency,
C               input in order, output bit reversed
C
C  type of FFT.  It is described in the book
C
C               THE FAST FOURIER TRANSFORM,
C               F. Oran Brigham, Prentice Hall, 1974
C
C
C  BUGS
C
C      The bit-reversal routine could be made faster by using the
C      VAX FORTRAN 77 library routines for testing, setting, clearing and
C      shifting bits:  BTEST, IBSET, IBCLR, ISHFT.
C
C      You might want to use double precision complex numbers to get more
C      accuracy.
C
C  ****************************************************************************
C  *                                                                          *
C  ****************************************************************************


        SUBROUTINE FFT( N, X, D )

        IMPLICIT NONE   ! Force declaration of all variables.


C    -------------------------------------------------------------------------
C    |                            MACRO DEFINITIONS                          |
C    ------------------------------------------------------------------------- 

        INTEGER MAX_TABLE_SIZE
        REAL    PI

        PARAMETER ( PI = 3.1415926535897932385 )

C       We must set the size of the sine and cosine table used in the FFT.
C       The table is called starting_multiplier.
C     
C       Suppose MAX is the largest FFT we would ever want
C       to perform.  Then 
C
C       MAX_TABLE_SIZE  =  LOG ( MAX ) - 1
C                             2
C
C       is the size of the table.  This table can be used by the
C       FFT for any number of points from 2 up to MAX.
C
C       For example, if MAX_TABLE_SIZE = 14, then we can transform
C       anywhere from 2 to 2 ** 15 = 32,768 points, using the same 
C       sine and cosine table.
C

        PARAMETER ( MAX_TABLE_SIZE  = 20 )


C -------------------------------------------------------------------------
C |                          CALLING ARGUMENTS                            |
C -------------------------------------------------------------------------

        INTEGER N            !  Number of points in transform

        COMPLEX X( 0 : * )   !  Complex data to be transformed

C    Array X is dimensioned in the calling program.  Only its starting 
C    address is passed into this subroutine.

        INTEGER D            !  Transform direction:  1 for a forward 
                             !  transform, -1 for an inverse transform.


C -------------------------------------------------------------------------
C |                          LOCAL VARIABLES                              |
C -------------------------------------------------------------------------

C  A column refers to a column of the FFT signal flow graph.  A group is
C  a collection of butterflies which have the same multiplier W^p.  A
C  butterfly is the basic computational unit of the FFT.

        INTEGER log_of_n,          ! Log base 2 of the number of points, N.
     $          column,            ! Current column.
     $          number_of_groups,  ! Number of groups in a column.
     $          group,             ! Current group number.
     $          start_of_group,    ! Array index of first butterfly in group.
     $          end_of_group,      ! Array index of last butterfly in group.
     $          butterfly_height,  ! Height of a butterfly.
     $          offset,            ! Lower branch of a butterfly.
     $          i,                 ! Loop counter.  Upper branch of butterfly.
     $          rev_i,             ! Bit-reversed array index i.
     $          REVERSE_BITS       ! Bit reversal function.


        COMPLEX product,     ! Temporary storage in a butterfly calculation.
     $          temp,        ! Same as above.
     $          increment,   ! Difference between current multiplier and next.
     $          multiplier   ! The multiplier W^p for a group of butterflies.

        COMPLEX starting_multiplier( 0 : MAX_TABLE_SIZE )   ! Table of sines
                                                            ! and cosines.

        REAL normalizer,     ! 1/\/N normalizing factor.
     $       angle           ! Current angle in sine and cosine table.


        LOGICAL called_already   ! .TRUE. if this is the second call or
                                 ! later call to this subroutine.


C -------------------------------------------------------------------------
C |                          SUBROUTINE CODE                              |
C -------------------------------------------------------------------------

C
C  Initialize the variable called_already to .FALSE. at compile time.
C  Change it to .TRUE. after the first call to this subroutine.
C  Save the value so it does not get destroyed after returning from
C  this subroutine.
C
        SAVE called_already
        DATA called_already /.FALSE./


C
C       Compute log base 2 of N. Perturb slightly to avoid roundoff error.
C       For example, INT( 4.99999999999 ) = 4 (wrong!), but
C       INT( 4.999999999 + .01 ) = 5 (correct!)
C
        log_of_n = INT( ALOG( FLOAT( N ) ) / ( ALOG( 2.0 ) )  + .01 )



C
C  Error checking.  Print a message to the terminal if any one of these
C  conditions is true:
C
C      (1) D is not 1 or -1.
C      (2) N < 2.
C      (3) log2( N ) > MAX_TABLE_SIZE + 1.
C      (4) N is not a power of 2.
C

        IF ((D .NE. 1) .AND. (D .NE. -1)) THEN

            PRINT *, ' ERROR ---  D is not 1 or -1'
            RETURN

        ELSE IF (N .LT. 2) THEN

            PRINT *, ' ERROR ---  N is less than 2. '
            RETURN

        ELSE IF ( (log_of_n - 1) .GT. MAX_TABLE_SIZE) THEN

            PRINT *, ' ERROR ---  N is larger than the sine table can'
            PRINT *, ' handle.  Increase the parameter MAX_TABLE_SIZE'
            PRINT *, ' in this subroutine and recompile.'

            RETURN

        ELSE IF ( (2 ** log_of_n) .NE. N ) THEN

            PRINT *, ' ERROR --- N is not a power of 2.'
            RETURN

        END IF


C
C  If called_already is false, generate the sine and cosine table for
C  the first time, then set called_already to .TRUE.  Now later calls to 
C  the subroutine don't recreate the sine and cosine table.
C

        IF (.NOT. called_already) THEN

            DO i = 0, MAX_TABLE_SIZE - 1

                angle = - PI / ( 2 ** i )
                starting_multiplier( i ) = CEXP( CMPLX(0.0 , angle) )

            END DO

            called_already = .TRUE.

        END IF



C
C       Compute the FFT according to its signal flow graph with 
C       computations arranged into butterflies, groups of butterflies
C       with the same multiplier, and columns.
C
        butterfly_height = N / 2
        number_of_groups = 1


        DO column = 1, log_of_n

            multiplier = starting_multiplier( column - 1 )

C
C       Modify starting multiplier W^p when computing the inverse transform.
C
            IF ( D .EQ. -1 ) THEN

                multiplier = CONJG( multiplier )

            END IF

            increment = multiplier


C
C       Compute the first group.
C
            DO i = 0, butterfly_height - 1

                offset = i + butterfly_height

                temp        =  X( i ) + X( offset )
                X( offset ) =  X( i ) - X( offset )
                X( i )      =  temp

            END DO  ! i


C
C       Now do the remaining groups.
C
            DO group = 2, number_of_groups

                start_of_group = REVERSE_BITS( group - 1, log_of_n )
                end_of_group = start_of_group +
     $                         butterfly_height - 1


C       Compute all butterflies in a group.

                DO i = start_of_group, end_of_group

                    offset   =  i + butterfly_height
                    product  =  x( offset ) * multiplier

                    temp         = X( i ) + product
                    X( offset )  = X( i ) - product
                    X( i )       = temp

                END DO  ! i

                multiplier = multiplier * increment

            END DO ! group

            butterfly_height = butterfly_height / 2
            number_of_groups = number_of_groups * 2

        END DO ! column


C
C       Normalize the results and permute them into order, two at a time.
C
        normalizer = 1 / SQRT( FLOAT( N ) )


        DO i = 0, N - 1

            rev_i = REVERSE_BITS( i, log_of_n )

            IF ( rev_i .GT. i ) THEN   ! exchange not yet done.

                temp       =  X( i )
                X( i )     =  X( rev_i ) * normalizer
                X( rev_i ) =  temp       * normalizer

           ELSE IF ( rev_i .EQ. i ) THEN  ! exchange with itself.

                X( i ) = X( i ) * normalizer

           END IF

        END DO ! i


        RETURN
        END   !  subroutine FFT



C  ****************************************************************************
C  *                                                                          *
C  *                            REVERSE_BITS                                  *
C  *                                                                          *
C  ****************************************************************************
C
C  DESCRIPTION
C
C      This function reverses the order of the bits in a number.
C
C  
C  INPUT
C      
C       n              The number whose bits are to be reversed.
C
C       num_bits       The number of bits in N, counting leading zeros.
C
C
C  OUTPUT
C  
C       REVERSE_BITS   The number N with all its num_bits bits reversed.
C
C
C  EXAMPLE
C
C       REVERSE_BITS( (10)  , 2 ) = (01) = 1, but
C                         2             2
C
C       REVERSE_BITS( (10)  , 3 ) = REVERSE_BITS( (010)  , 3 ) = (010)  = (10)
C                         2                            2              2       2
C
C       = 10
C
C
C  METHOD
C
C       The algorithm works by a method similar to base conversion.
C       As the low order bits are broken off, one by one, they are
C       multiplied by 2 and accumulated to generate the bit-reversed
C       number.
C
C
C  BUGS
C
C      The bit-reversal routine could be made faster by using the
C      VAX FORTRAN 77 library routines for testing, setting, clearing and
C      shifting bits:  BTEST, IBSET, IBCLR, ISHFT.
C
C
C  ****************************************************************************
C  *                                                                          *
C  ****************************************************************************

        INTEGER FUNCTION REVERSE_BITS( n, num_bits )

        IMPLICIT NONE


C -------------------------------------------------------------------------
C |                          CALLING ARGUMENTS                            |
C -------------------------------------------------------------------------

        INTEGER n, num_bits


C -------------------------------------------------------------------------
C |                          LOCAL VARIABLES                              |
C -------------------------------------------------------------------------

        INTEGER bit_num,      ! Current bit number.
     $          remainder,    ! Bits remaining to be reversed.
     $          reversed_n,   ! Bit-reversed value of n.
     $          bit,          ! The current least significant bit.
     $          MOD, INT      ! Function return types.


C -------------------------------------------------------------------------
C |                          SUBROUTINE CODE                              |
C -------------------------------------------------------------------------

        remainder = n                      ! Load n.
        reversed_n = 0                     ! Clear buffer for bit-reversed n.

        DO bit_num = 1, num_bits

            bit        = MOD( remainder, 2 )   ! Next least significant bit.
            reversed_n = bit + 2 * reversed_n  ! Shift left and add to buffer.
            remainder  = INT( remainder / 2 )  ! Remaining bits.

        END DO

        REVERSE_BITS = reversed_n

        RETURN
        END     !  subroutine REVERSE_BITS