1
$\begingroup$

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

enter image description here

but working for it in Fortran is so complicated!

$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$ Mar 23, 2023 at 7:19
  • $\begingroup$ I would create a short test array with data of known frequency, i.e. fill an array with sin(i*omega) and then debug $\endgroup$
    – MPIchael
    Mar 23, 2023 at 7:21
  • 3
    $\begingroup$ 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. $\endgroup$ Mar 23, 2023 at 16:00
  • 2
    $\begingroup$ @LutzLehmann, would you mind expanding your comment into an answer? $\endgroup$
    – nicoguaro
    Mar 23, 2023 at 16:32

1 Answer 1

1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.