# How to get the inverse FFt in this Fortran code?

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!

• Thank you to getting the code question computational experts. I would still ask you to show us the inverse transform you tried and the non-matching results you got. One would normally just change the sign in exp and normalize. Commented Mar 23, 2023 at 7:19
• I would create a short test array with data of known frequency, i.e. fill an array with sin(i*omega) and then debug Commented Mar 23, 2023 at 7:21
• Apply the forward FFT twice and divide by N, and you get back to your original data, but in a reflected spacial arrangement. So you could implement the iFFT with minimal effort. Or introduce a direction parameter for the sign in the exponents so that fft(x,-1) is the analysis direction and fft(x,+1)/N is the synthesis direction. Commented Mar 23, 2023 at 16:00
• @LutzLehmann, would you mind expanding your comment into an answer? Commented Mar 23, 2023 at 16:32

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.