How to get the inverse FFt in this Fortran code?

Question

I find this fft algorithm on the link

The code looks simple and easy to implement. But it does not have inverse fast Fourier transformation. A brief search on the internet shows that to get the inverse of it, the data should be divided by the length of the array. Dividing this, however, does not produce the original input data. So which part is missing in getting the inverse of FFT?

module fft_mod
  implicit none
  integer,       parameter :: dp=selected_real_kind(15,300)
  real(kind=dp), parameter :: pi=3.141592653589793238460_dp
contains

  ! In place Cooley-Tukey FFT
  recursive subroutine fft(x)
    complex(kind=dp), dimension(:), intent(inout)  :: x
    complex(kind=dp)                               :: t
    integer                                        :: N
    integer                                        :: i
    complex(kind=dp), dimension(:), allocatable    :: even, odd

    N=size(x)

    if(N .le. 1) return

    allocate(odd((N+1)/2))
    allocate(even(N/2))

    ! divide
    odd =x(1:N:2)
    even=x(2:N:2)

    ! conquer
    call fft(odd)
    call fft(even)

    ! combine
    do i=1,N/2
       t=exp(cmplx(0.0_dp,-2.0_dp*pi*real(i-1,dp)/real(N,dp),kind=dp))*even(i)
       x(i)     = odd(i) + t
       x(i+N/2) = odd(i) - t
    end do

    deallocate(odd)
    deallocate(even)

  end subroutine fft

end module fft_mod

program test
  use fft_mod
  implicit none
  complex(kind=dp), dimension(8) :: data = (/1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0/)
  integer :: i

  call fft(data)

  do i=1,8
     write(*,'("(", F20.15, ",", F20.15, "i )")') data(i)
  end do

end program test

I did write

data = data/8       ! 8 is the dimension

In matlab, i could just use it like

but working for it in Fortran is so complicated!

Answer 1

The theory of the discrete Fourier transform is based on the geometric sum formula $$ 1+q+...+q^{n-1}=\begin{cases}0,&q^n=1 \\ n,&q=1 \\ \frac{1-q^n}{1-q},& else.\end{cases} $$ The last case does not occur explicitly in the DFT, as the $q$ are $n$th unit roots and their powers.

The DFT of a sequence $x_0,...,x_{n-1}$ gives the numbers $$ \hat x_k=\sum_{j=0}^{n-1}x_jq^{jk} $$ The interesting magic happens if the transform is applied again, as then \begin{align} \hat{\hat x}_m&=\sum_{k=0}^{n-1}\sum_{j=0}^{n-1}x_jq^{jk}q^{km} \\ &=\sum_{j=0}^{n-1}x_j\sum_{k=0}^{n-1}q^{(j+m)k} =\begin{cases}nx_0&m=0,\\nx_{n-m},&m>0.\end{cases} \end{align} So essentially the inverse transform is obtained as $x_0=\frac1n\hat{\hat x}_0$ and $x_j=\frac1n\hat{\hat x}_{n-j}$.

Now observing that $q^{n-j}=q^{-j}=\bar q^j$ one can also formulate the inverse transform as extra procedure $$ x_j=\frac1n\sum_{k=0}^{n-1}\hat x_k q^{-jk}. $$

So how you implement the inverse transform is a matter of taste, you can use the same FFT procedure and apply the indicated post-processing, or copy the FFT routine and just modify the angle factors and divide by $n$ at some point, or modify the FFT to have a switch that applies these modification in one procedure.