[通用算法] 二维傅里叶变换

Module FFT2D_Mod

  Implicit None

  Integer, parameter :: DP = Selected_Real_Kind(15)

  Integer, parameter :: FFT_Forward = 1

  Integer, parameter :: FFT_Inverse = -1

contains

  Subroutine fc2DFFT(x, forback)

    ! 2D Fast Fourier Transform using Cooley-Tukey algorithm

    Real(Kind=DP), parameter :: PI = 3.141592654_DP

    Complex(Kind=DP), Intent(INOUT) :: x(:,:)

    Integer, Intent(IN) :: forback

    Integer :: nx, ny

    Integer :: i, j, k, ncur, ntmp, itmp

    Real(Kind=DP) :: ex, ey

    Complex(Kind=DP) :: ctmp

    nx = Size(x, 1)

    ny = Size(x, 2)

   

    ! Perform FFT along rows

    Do j = 1, ny

      ncur = nx

      Do

        ntmp = ncur

        ex = 2.0_DP * PI / ncur

        ncur = ncur / 2

        if (ncur < 1) exit

        Do i = 1, nx, ntmp

          Do k = 1, ntmp / 2

            itmp = i + k

            ctmp = x(itmp, j) * Exp(forback * CMPLX(0.0_DP, ex * (k - 1), Kind=DP))

            x(itmp, j) = x(i, j) - ctmp

            x(i, j) = x(i, j) + ctmp

          End Do

        End Do

      End Do

    End Do

   

    ! Perform FFT along columns

    Do i = 1, nx

      ncur = ny

      Do

        ntmp = ncur

        ey = 2.0_DP * PI / ncur

        ncur = ncur / 2

        if (ncur < 1) exit

        Do j = 1, ny, ntmp

          Do k = 1, ntmp / 2

            itmp = j + k

            ctmp = x(i, itmp) * Exp(forback * CMPLX(0.0_DP, ey * (k - 1), Kind=DP))

            x(i, itmp) = x(i, j) - ctmp

            x(i, j) = x(i, j) + ctmp

          End Do

        End Do

      End Do

    End Do

   

    ! If performing inverse FFT, divide by the total number of points

    if (forback == FFT_Inverse) then

      x = x / (nx * ny)

    end if

   

  End Subroutine fc2DFFT

End Module FFT2D_Mod

Program TwoDFFT

  use FFT2D_Mod

  Implicit None

  Integer :: i, j

  Integer, parameter :: Nx = 4

  Integer, parameter :: Ny = 4

  Complex(Kind=DP) :: x(Nx, Ny) = reshape([36.d0, 21.d0, 33.d0, 44.d0, 55.d0, 63.d0, 73.d0, 38.d0, &

                                            36.d0, 21.d0, 33.d0, 44.d0, 55.d0, 63.d0, 73.d0, 38.d0], [Nx, Ny])

  

  ! Perform 2D FFT

  Call fc2DFFT(x, FFT_Forward)

  

  ! Print the result

  Do j = 1, Ny

    Do i = 1, Nx

      Write(*, *) x(i, j)

    End Do

  End Do

  

  ! Perform inverse 2D FFT

  Call fc2DFFT(x, FFT_Inverse)

  

  ! Normalize and print the result

  Do j = 1, Ny

    Do i = 1, Nx

      Write(*, *) x(i, j)

    End Do

  End Do

  

End Program TwoDFFT