Fortran Coder

标题: 二维傅里叶变换 [打印本页]
作者: zhaomin    时间: 2023-10-8 17:08
标题: 二维傅里叶变换
可以帮忙看一下这段代码哪里出问题了嘛，算出来的数据不太对
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





欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2