intel fortran调用FFTW，傅里叶正变换值N/2+1之后为0

放耗子的绅 · 发表于 2019-11-20 11:32:58

各位大佬们，请教一个问题，我用intel fortran去调用FFTW，但傅里叶正变换值N/2+1之后为0，通过对比matlab算的结果，发现偏差很大；
麻烦问一下，如果我要得到N/2+1以后的数据，应该怎么做呢？代码如下：
MODULE globalvars
IMPLICIT NONE
! Simulation cell parameters:
integer,parameter:: Nx = 4, Ny = 4
integer(kind = 8) :: p_up, p_dn
END MODULE globalvars

!!!-------------------傅里叶正变换
subroutine backward(output,input)
USE globalvars
!calculates backward fourier transform (complex to real) using FFTW
implicit none
COMPLEX, INTENT(IN) :: output(nx,ny)
REAL, INTENT(OUT) :: input(nx,ny)
COMPLEX :: tempk(nx,ny)
INTEGER :: i,j
REAL :: sizescale

sizescale = 1.0/FLOAT(nx * ny)
tempk = 0.0
DO j=1,ny
DO i=1,nx
tempk(i,j) = output(i,j)
END DO
END DO

CALL sfftw_execute_dft_c2r(p_dn,tempk,input)
DO j=1,ny
DO i=1,nx
input(i,j) = sizescale*input(i,j)
END DO
END DO
! RETURN
END subroutine backward

!!!-------------------傅里叶逆变换
subroutine forward(input,output)
USE globalvars
!calculates forward fourier transform (real to complex) using FFTW
implicit none

real, intent(IN) :: input(nx,ny)
complex, intent(OUT) :: output(nx,ny)

call sfftw_execute_dft_r2c(p_up,input,output)
! RETURN
end subroutine forward

!---------------------主程序

program main

USE globalvars

implicit none
include 'fftw3.f'
integer :: i,j
real,allocatable  :: eta(:,:), input(:,:)
complex,allocatable :: etak(:,:), output(:,:)
allocate (eta(nx,ny), input(nx,ny))
allocate (etak(nx,ny), output(nx,ny))
eta = reshape([3,0,-1,0,&
   -1,-2,1,1,&
   -2,2,2,2,&
   2,-1,4,3],shape(eta))
do i=1,nx
do j=1,ny
write(*,*) eta(j,i)
enddo
enddo
CALL sfftw_plan_dft_r2c_2d(p_up, nx, ny, input, output, FFTW_EXHAUSTIVE)
CALL sfftw_plan_dft_c2r_2d(p_dn, nx, ny, output, input, FFTW_EXHAUSTIVE)

call forward(eta,etak)
do i=1,nx
do j=1,ny
write(*,*) etak(j,i)
enddo
enddo

call backward(etak,eta)
do i=1,nx
do j=1,ny
write(*,*) eta(j,i)
enddo
enddo
deallocate (eta, input)
deallocate (etak, output)
end program

运行结果见图。

li913 · 发表于 2019-11-23 23:33:40

本帖最后由 li913 于 2019-11-23 23:35 编辑

结果没错，只是这个函数省略了共轭对称的元素，你把它恢复就是。

放耗子的绅 · 发表于 2019-12-4 11:11:27

li913 发表于 2019-11-23 23:33
结果没错，只是这个函数省略了共轭对称的元素，你把它恢复就是。

非常感谢您的解答，我测试了下，结果没问题！

青衣巷 · 发表于 2019-12-4 22:40:41

时间域频率域信息量不变哟~

[数学库] intel fortran调用FFTW，傅里叶正变换值N/2+1之后为0

规矩勋章

新人勋章

元老勋章