Basic use of FFTPACK

[zoziha]

Yes, fortran-lang’s fftpack is based on netlib’s dfftpack 4.0, they all only support FFT of rank-1 arrays, and it may be possible to expand FFT of higher rank in the future, here is a simple example to support fft2, you can learn from this code, extend FFT for higher rank.
Or simply do this

use fftpack
do i = 1, n
    y(:,i) = fft(x(:,i))  !! rank2-fft step1: do fft for dim=1
end do

do i = 1, m
    y(i, :) = fft(y(i,:)) !! rank2-fft step2: do fft for dim=2
end do

In addition, you can also use ncar’s fftpack5.1 or fftw’s fortran interface, which supports fft, fft2, fft3, fftn.

[Shahid]

Ok. Now I want to do fft of the random number. My code is

program rand2D_test
    use fftpack   
    use fftpack_kind, only: rk
    implicit none

    integer(kind=4), parameter    :: Nx = 4, Ny = 4
    real(kind=8),dimension(Nx,Ny) :: r
    integer(kind=4)               :: i,j

    call random_number(r)

    do i=1,Nx
        r(:,i) = fft(r(:,i))
    end do
    do j=1,Ny
       r(j,:) = fft(r(i,:))
    end do

    do i=1,Nx
       print*,(r(i,j), j=1,Ny)
    end do

end program rand2D_test

I get this error:

[mecej4]

The interface name fft represents one or more specific functions that are defined in the module source file(s).

Look through fftpack.f90, and see if you can find the interface name fft. If you do, see if the argument(s) to the function match the actual argument type(s) (array of type real).

If there is no match, look for another function with the same interface name, and see if it matches. If you can find no such function, that explains the error message.

[zoziha]

Currently, the fft input argument here only supports complex(8) types:

program rand2D_test
    use fftpack   
    use fftpack_kind, only: rk
    implicit none

    integer(kind=4), parameter    :: Nx = 4, Ny = 4
    complex(kind=rk),dimension(Nx,Ny) :: r
    integer(kind=4)               :: i,j

    call random_number(r%re)
    r%im = 0.0_rk

    do i=1,Nx
        r(:,i) = fft(r(:,i))
    end do
    do j=1,Ny
       r(j,:) = fft(r(i,:))
    end do

    do i=1,Nx
       print*,(r(i,j), j=1,Ny)
    end do

end program rand2D_test

From a practical point of view, adding fft’s support for real numbers (real(8)) may be more user-friendly; but if only complex numbers are supported here, the consistency of fft and ifft interfaces can be maintained.

[Shahid]

Ok. That is helpful. Thanks. How about the inverse fft in this program?

program rand2D_test
   use fftpack   
   use fftpack_kind, only: rk
   implicit none
 
   integer(kind=4), parameter        :: Nx = 2, Ny = 2
   complex(kind=rk),dimension(Nx,Ny) :: r
   integer(kind=4)                   :: i,j

  call random_number(r%re)
  r%im = 0.0_rk

  print*, ''
  print*, 'Random Numbers are'
  print*, ''

  do i=1,Nx
     print*,(r(i,j), j=1,Ny)
  end do

  do i=1,Nx
     r(:,i) = fft(r(:,i))
  end do
  do j=1,Ny
     r(j,:) = fft(r(i,:))
  end do

  print*, ''
  print*, 'FFT is'
  print*, ''

  do i=1,Nx
     print*,(r(i,j), j=1,Ny)
  end do

  do i=1,Nx
     r(:,i) = ifft(r(:,i))
  end do
  do j=1,Ny
     r(j,:) = ifft(r(i,:))
  end do

  print*, ''
  print*, 'Inverse FFT is'
  print*, ''

  do i=1,Nx
     print*,(r(i,j), j=1,Ny)
  end do

end program rand2D_test

By running it I am not getting the same original values. Is ifft here not the desired function for inverse fft?

[zoziha]

I’m sorry, I copied your first code directly and didn’t notice that the index was not written correctly, I corrected it now.
Perform fft on rank-1 arrays of the same length, fft internal working array wsave has the same initialization, if you want to avoid it and improve the running efficiency, please use zffti, zfftf, not fft.

program rand2D_test
    use fftpack
    use fftpack_kind, only: rk
    implicit none

    integer(kind=4), parameter :: Nx = 2, Ny = 2
    complex(kind=rk), dimension(Nx, Ny) :: r
    real(kind=rk) :: wsave(4*Nx + 15)
    integer(kind=4) :: i, j

    call random_number(r%re)
    r%im = 0.0_rk

    print *, ''
    print *, 'Random Numbers are'
    print *, ''

    do i = 1, Nx
        print *, (r(i, j), j=1, Ny)
    end do

    do i = 1, Nx
        r(:, i) = fft(r(:, i))
    end do
    do j = 1, Ny
        r(j, :) = fft(r(j, :)) ! i -> j
    end do

    print *, ''
    print *, 'FFT is'
    print *, ''

    do i = 1, Nx
        print *, (r(i, j), j=1, Ny)
    end do

    do i = 1, Nx
        r(:, i) = ifft(r(:, i))
    end do
    do j = 1, Ny
        r(j, :) = ifft(r(j, :)) ! i -> j
    end do

    print *, ''
    print *, 'Inverse FFT is'
    print *, ''

    r = r/Nx/Ny ! normalize
    do i = 1, Nx
        print *, (r(i, j), j=1, Ny)
    end do

    call zffti(Nx, wsave)
    do i = 1, Nx
        call zfftf(Ny, r(i, :), wsave)
    end do

    do i = 1, Ny
        call zfftf(Nx, r(:, i), wsave)
    end do

    print *, ''
    print *, 'zfftf is'
    print *, ''
    
    do i = 1, Nx
        print *, (r(i, j), j=1, Ny)
    end do

end program rand2D_test

 Random Numbers are

              (0.96918105010483302,0.0000000000000000)              (0.61750876263731336,0.0000000000000000)
              (0.45185016840516656,0.0000000000000000)              (0.84288263491537585,0.0000000000000000)

 FFT is

               (2.8814226160626886,0.0000000000000000)        (-3.93601790426896248E-002,0.0000000000000000)
              (0.29195700942160396,0.0000000000000000)              (0.74270475397772895,0.0000000000000000)

 Inverse FFT is

              (0.96918105010483302,0.0000000000000000)              (0.61750876263731325,0.0000000000000000)
              (0.45185016840516651,0.0000000000000000)              (0.84288263491537574,0.0000000000000000)

 zfftf is

               (2.8814226160626886,0.0000000000000000)        (-3.93601790426894582E-002,0.0000000000000000)
              (0.29195700942160396,0.0000000000000000)              (0.74270475397772895,0.0000000000000000)

[zoziha]

When I started to learn to use fftpack, I was also confused. I only need to use fft and ifft in octave/MatLab, without fordward and backword. Why does fftpack need to initialize coefficients, forward fft, and backword fft? This seems to be an efficiency issue.

[Shahid]

I have the same idea. I am trying to use fft and ifft for my code originally written in Octave/Matlab. It is so convenient there to use it. Just type fft and the answer is there! Here, in Fortran, so far it took more than two days and I am still learning. But thanks to your help. I think I can read it more and learn from it. In case, I would let you know!