program FourierTrans
  implicit none

  integer i
  real,parameter:: pi=3.14159265358979
  integer,parameter:: m=17
  integer,parameter:: N=2**m			! size of vectors, power of 2
  real a(2,0:N-1)			        !  input vector
  real b(2,0:N-1)				! output vector
  real tb,te                                    
  integer :: ic
  logical :: yes = .true.
  !----------------------------------------------------------------
  call random_seed()

  !*** random generation of vector 'a'***
  call random_number(a(1,0:N-1)); a(2,0:N-1)=0.0
  a = 0.
  do i=0,N-1
!     a(1,i) = 1.
!     a(1, i) = sin(2. *pi *i / n) 
!     a(1, i) = sin(2. *pi *i / n) + 0.001* sin(20. *pi *i / n) 
          !+ 20* sin(250. *pi *i / n) + 0.001*s
!     if( i > 2*n/5 .and. i  < 3*n/5) a(1,i) = 1.
     if( i > 1*n/7 .and. i  < 2*n/7) a(1,i) = 1.
     if( i > 3*n/7 .and. i  < 4*n/7) a(1,i) = 1.
     if( i > 5*n/7 .and. i  < 6*n/7) a(1,i) = 1.
!     if( i < 2*n/5 .or. i  > 3*n/5) a(1,i) = 1.
     write(1,*) a(:,i), i, sqrt(a(1,i)**2 + a(2,i)**2)
  enddo
  print*,'N = ', N

  !***  (no fast) Fourier transformation ***
  if(.not. yes) then
     !forward transformation 'a' to 'b'
     call cpu_time(tb)
     b=DFT(N,a,1)
     call cpu_time(te)
     write(*,'(a30,f14.6,a2)') 'CPU time for DFT forward:  ', te-tb, ' s'
     do i=0,N-1
        write(10,*) b(:,i),i, sqrt(b(1,i)**2 + b(2,i)**2)
     enddo

     !bacward transformation of'b'
     call cpu_time(tb)
     b=DFT(N,b,-1)
     call cpu_time(te)
     write(*,'(a30,f14.6,a2)') 'CPU time for DFT backward: ', te-tb, ' s'
     do i=0,N-1
        write(11,*) b(:,i),i, sqrt(b(1,i)**2 + b(2,i)**2)
     enddo
  endif

  !*** FFT -recursive ***
  if(yes) then
    !forward transformation 'a' to 'b'
      call cpu_time(tb)
     call FFTrec(N,a,b,1)
     call cpu_time(te)
     write(*,'(a30,f14.6,a20)') 'CPU time for FFT forward:  ', te-tb, ' s  (with recursion)'
     do i=0,N-1
        write(20,*) b(:,i),i, sqrt(b(1,i)**2 + b(2,i)**2)
     enddo

     ic = 0 
     do i=0,N-1
        if( norm2(b(1:2,i)) < 10.) then
           b(:,i) = 0.0
           ic = ic + 1
        endif
        
     enddo

     write(*,'(a30, 2i9, f6.3)') 'Size, reduction, ratio::', N, ic, 1.*ic/N
     
     !bacward transformation of'b'
     call cpu_time(tb)
     call FFTrec(N,b,b,-1)
     call cpu_time(te)
     write(*,'(a30,f14.6,a20)') 'CPU time for FFT backward: ', te-tb, ' s  (with recursion)'
     do i=0,N-1
        write(21,*) b(:,i),i, sqrt(b(1,i)**2 + b(2,i)**2)
     enddo
  endif

  !*** FFT - direct method***
  if(yes) then
      !forward transformation 'a' to 'b'
     call cpu_time(tb)
     b=FFT(N,a,1)
     call cpu_time(te)
     write(*,'(a30,f14.6,a23)') 'CPU time for FFT forward:  ', te-tb, ' s  (without recursion)'

     do i=0,N-1
        write(30,*) b(:,i),i, sqrt(b(1,i)**2 + b(2,i)**2)
     enddo

     !bacward transformation of'b'
     call cpu_time(tb)
     b=FFT(N,b,-1)
     call cpu_time(te)
     write(*,'(a30,f14.6,a23)') 'CPU time for FFT backward: ', te-tb, ' s  (without recursion)'
     do i=0,N-1
        write(31,*) b(:,i),i, sqrt(b(1,i)**2 + b(2,i)**2)
     enddo
  endif


contains
  !================================================================
  function DFT(N,a,i)		! compute discrete Fourier trasformation of 'a' 
    !i == 1 forward, i= -1 backward
    integer,intent(in):: N,i
    real,intent(in):: a(2,0:N-1)	! index (1,:) -- real part,  index(2,:) -- imaginary part
    real DFT(2,0:N-1)			! index (1,:) -- real part,  index(2,:) -- imaginary part
    real factor,arg,re,im
    integer j,k
    if(i/=1 .and. i/=-1) stop 'Factor i must be from {-1;1}'
    factor=6.28318530717959/N	!2*PI/N
    do j=0,N-1
       DFT(:,j)=0.0
       do k=0,N-1
          arg=factor*j*k
          re=cos(arg); im=i*sin(arg);
          DFT(1,j)=DFT(1,j)+re*a(1,k)-im*a(2,k)
          DFT(2,j)=DFT(2,j)+im*a(1,k)+re*a(2,k)
       enddo
    enddo
    if(i==-1) DFT=DFT/N
  endfunction DFT
  !================================================================
  recursive subroutine FFTrec(N,a,b,i)	! compute fast Fourier trasformation of 'a' by recursion
    !i == 1 forward, i= -1 backward
    ! FFT with recursion --> N must be power of 2
    integer,intent(in):: N,i
    ! the following arrays:  index (1,:) -- real part,  index(2,:) -- imaginary part
    real,intent(in):: a(2,0:N-1)       	
    real,intent(out):: b(2,0:N-1)       
    real s(2,0:N/2-1),l(2,0:N/2-1)	! even (s) and odd (l) components of input arrays
    real factor,arg,re,im
    integer j
    if(iand(N,N-1)/=0) stop 'N is not the power of 2'
    if(i/=1 .and. i/=-1) stop 'Factor i must be from  {-1;1}'
    if(N==1) then
       b=a
       return
    endif
    factor=6.28318530717959/N	!2*PI/N
    call FFTrec(N/2,a(:,0:N-2:2),s,i)	! recursive calling of for even components
    call FFTrec(N/2,a(:,1:N-1:2),l,i)	! recursive calling of for odd components
    do j=0,N/2-1
       arg=factor*j
       re=cos(arg); im=i*sin(arg);
       !*** first part of input  array
       b(1,j)=s(1,j)+re*l(1,j)-im*l(2,j)
       b(2,j)=s(2,j)+im*l(1,j)+re*l(2,j)
       !*** second part of input  array
       b(1,j+N/2)=s(1,j)-(re*l(1,j)-im*l(2,j))
       b(2,j+N/2)=s(2,j)-(im*l(1,j)+re*l(2,j))
    enddo
    if(i==-1) b=b/2
  endsubroutine FFTrec
  !================================================================
  function FFT(N,a,i)	! compute fast Fourier trasformation of 'a' without recursion
    !i == 1 forward, i= -1 backward
    ! FFT with recursion --> N must be power of 2
    integer,intent(in):: N,i
    ! the following arrays:  index (1,:) -- real part,  index(2,:) -- imaginary part
    real,intent(in):: a(2,0:N-1)
    real FFT(2,0:N-1)
    real FFT1(2),FFT2(2),FFTtemp(2)
    integer j,k,l,m
    integer n_i,n_d
    integer n_s,n_e
    real factor,arg,re,im
    if(iand(N,N-1)/=0 .or. N<1) stop 'N is not the power of 2'
    if(i/=1 .and. i/=-1) stop 'Factor i must be from  {-1;1}'
    !Finding the power of N
    do j=0,bit_size(N)-2
       if(btest(N,j)) then
          m=j
          exit
       endif
    enddo
    !***FFT bit reversion ***
    do j=0,N-1
       k=bit_reversed(m,j)
       if(k==j) then
          FFT(:,j)=a(:,j)
       elseif(k>j) then
          FFT(:,j)=a(:,k)
          FFT(:,k)=a(:,j)
       endif
    enddo
    !*** transformation***
    n_d=N
    n_i=1
    do j=1,m	!transformation of size block  2, 4, 8, ..., N
       n_d=n_d/2
       n_s=0
       n_e=n_i
       factor=3.14159265358979/n_i	!2*PI/(size of the step)
       do k=1,n_d	!transformation for each block, their number is  N/2, N/4, N/8, ..., 1	
          do l=0,n_i-1	!transformation over pairs in block
             arg=factor*l
             re=cos(arg); im=i*sin(arg);
             FFTtemp(1)=re*FFT(1,n_e)-im*FFT(2,n_e)
             FFTtemp(2)=im*FFT(1,n_e)+re*FFT(2,n_e)
             FFT1=FFT(:,n_s)+FFTtemp
             FFT2=FFT(:,n_s)-FFTtemp
             FFT(:,n_s)=FFT1
             FFT(:,n_e)=FFT2
             n_s=n_s+1
             n_e=n_e+1
          enddo
          n_s=n_s+n_i
          n_e=n_e+n_i
       enddo
       n_i=n_i*2
    enddo
    if(i==-1) FFT=FFT/N
  endfunction FFT
  !================================================================
  integer function bit_reversed(m,number)
    !*** bit reversion 
    integer,intent(in):: m,number
    integer i,j
    i=number
    bit_reversed=0
    do j=1,m-1
       bit_reversed=(bit_reversed+mod(i,2))*2;
       i=i/2
    enddo
    bit_reversed=bit_reversed+mod(i,2)
  endfunction bit_reversed
  !================================================================
  integer function bit_reversed2(m,number)
    integer,intent(in):: m,number
    integer j
    bit_reversed2=0
    do j=0,m-1
       if(btest(number,j)) bit_reversed2=ibset(bit_reversed2,m-1-j)
    enddo
  endfunction bit_reversed2
  !================================================================

end program FourierTrans