二维傅里叶变换
可以帮忙看一下这段代码哪里出问题了嘛，算出来的数据不太对{:4_116:}Module FFT2D_ModImplicit NoneInteger, parameter :: DP = Selected_Real_Kind(15)Integer, parameter :: FFT_Forward = 1Integer, parameter :: FFT_Inverse = -1contains 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 TwoDFFTuse FFT2D_ModImplicit NoneInteger :: i, jInteger, parameter :: Nx = 4Integer, parameter :: Ny = 4Complex(Kind=DP) :: x(Nx, Ny) = reshape( 36.d0, 21.d0, 33.d0, 44.d0, 55.d0, 63.d0, 73.d0, 38.d0], )! Perform 2D FFTCall fc2DFFT(x, FFT_Forward)! Print the resultDo j = 1, Ny Do i = 1, Nx Write(*, *) x(i, j) End DoEnd Do! Perform inverse 2D FFTCall fc2DFFT(x, FFT_Inverse)! Normalize and print the resultDo j = 1, Ny Do i = 1, Nx Write(*, *) x(i, j) End DoEnd DoEnd Program TwoDFFT
页:
[1]