program fft0

     ! Radix-2 DIF FFT.  Compiler: GNU Fortran 4.6.3-2.   
     ! John Bryan, November 2016.

     implicit none
     complex, parameter :: j=(0,1)
     real, parameter  :: pi=3.14159265
     integer, parameter :: N = 64      ! fft length 
     complex :: xa(0:((N/2)-1))
     complex :: xs(0:(N-1))
     complex :: e
     complex :: o
     real    :: Ts = 1.0/N   ! sampling interval 1/Fs
     integer :: m = 0
     integer :: Np=N
     integer :: Nph
     integer :: BaseE
     integer :: BaseO
     integer :: tss=1           ! twiddle step size
     integer :: n2              ! Number of butterflies/sub-block
     integer :: Bp=1            ! Number of blocks 
     integer :: b
     integer :: p=int(log(real(N))/log(2.0))   ! Perform p passes  

     do m=0,(N-1)
       xs(m) = sin(2*pi*5*m*Ts)  ! sampled input
     end do
     do m=0,(N/2)-1
       xa(m)=exp(-2*j*pi*m/N)   ! N/2 twiddles 
     end do
     ! start DIF FFT algorithm
     do m=0,(p-1)               ! pass loop
         Nph=Np/2               ! # butterflies/sub-block halves 
         BaseE=0;
         do b=0,(Bp-1)          ! block loop
            BaseO=BaseE+Nph
            do n2=0,(Nph-1)     ! butterfly loop
               e=xs(BaseE+n2)+xs(BaseO+n2)
               o=(xs(BaseE+n2)-xs(BaseO+n2))*xa(n2*tss)
               xs(BaseE+n2)=e
               xs(BaseO+n2)=o
            end do
            BaseE=BaseE+Np;
         end do
         Bp=Bp*2                  ! # sub-blocks doubles 
         Np=Np/2                  ! # butterflies/sub-block halves 
         tss=tss*2                ! twiddle-step-size doubles
     end do

     call bitreversal()
     call writedata()

     contains

        subroutine bitreversal()
        ! xs changes: out-of-order to in-order array. 
        complex :: temp
        integer :: r(0:(N-1))
        integer :: M1=1
        integer :: T1
        integer :: n1
        integer :: k
        ! initialize r() in while loop
        do while (M1 < N)
             k = 0
             do while (k < M1)
                   T1 = 2*r(k)
                   r(k) = T1
                   r(k+M1) = T1+1
                   k = k+1
             end do
             M1 = 2*M1
        end do
        ! end of initializing r()
        ! r(n1) is the bitreversal of n1.
        ! For n1=1,2,...,N-2 do: if r(n1) > n1, swap xs(n1),xs(r(n1)).
        n1=1
        do while (n1 < (N-1))
            if (r(n1)>n1) then
               temp=xs(n1)
               xs(n1)=xs(r(n1))
               xs(r(n1))=temp
            end if
            n1=n1+1
        end do
        end subroutine bitreversal

        subroutine writedata()
        ! Write normalized magnitude single-sided output to file.
        real    :: xsa(0:(N-1))
        integer :: a
        do a=0,N-1
           xsa(a)=abs(xs(a)/N)
        end do
        open (unit=41, file='abs.txt')
        do  a=0,N/2-1
           write(41,'(F6.5)') xsa(a)
        end do
        call execute_command_line ("octave -q  plotfortranfft.octave")
        end subroutine writedata

end program fft0