使用 FFTW3 库评估 FORTRAN 中高斯函数的快速傅里叶变换

Question

我正在尝试编写一个 FORTRAN 代码来评估使用 FFTW3 库的高斯函数 f(r)=exp(-(r^2)) 的快速傅里叶变换。众所周知，高斯函数的傅里叶变换是另一个高斯函数。

我考虑评估球坐标中高斯函数的傅里叶变换积分。

因此得到的积分可以简化为[r*exp(-(r^2))*sin(kr)]dr的积分。

我编写了以下 FORTRAN 代码来评估离散正弦变换 DST，它是使用纯实输入数组的离散傅里叶变换 DFT。 DST由FFTW3中存在的C_FFTW_RODFT00执行，考虑到位置空间中的离散值是r=i*delta(i=1,2,...,1024)，DST的输入数组是函数r*exp(-(r^2))不是高斯函数。 [r*exp(-(r^2))*sin(kr)]dr 积分中的正弦函数由 SPHERICAL 坐标上的 INTEGRATION 得到，通常在进行解析傅里叶变换时出现的不是 exp(ik.r) 的虚部。

但是，结果不是动量空间中的高斯函数。

Module FFTW3
 use, intrinsic :: iso_c_binding
include 'fftw3.f03'
end module  

program sine_FFT_transform
use FFTW3
implicit none
integer, parameter :: dp=selected_real_kind(8)

real(kind=dp), parameter :: pi=acos(-1.0_dp)
integer, parameter :: n=1024 
real(kind=dp) :: delta, k
real(kind=dp) :: numerical_F_transform
integer :: i
type(C_PTR) ::  my_plan
real(C_DOUBLE), dimension(1024) :: y
real(C_DOUBLE), dimension(1024) :: yy, yk
integer(C_FFTW_R2R_KIND) :: C_FFTW_RODFT00

my_plan= fftw_plan_r2r_1d(1024,y,yy,FFTW_FORWARD, FFTW_ESTIMATE)

delta=0.0125_dp
do i=1, n        !inserting the input one-dimension position function
y(i)= 2*(delta)*(i-1)*exp(-((i-1)*delta)**2) 
! I multiplied by 2 due to the definition of C_FFTW_RODFT00 in FFTW3
end do

call fftw_execute_r2r(my_plan, y,yy)   
do i=2, n
k = (i-1)*pi/n/delta 
yk(i) = 4*pi*delta*yy(i)/2  !I divide by 2 due to the definition of 
                            !C_FFTW_RODFT00
numerical_F_transform=yk(i)/k
write(11,*) i,k,numerical_F_transform
end do
call fftw_destroy_plan(my_plan)

end program 

执行前面的代码会得到以下不适用于高斯函数的图。 谁能帮我理解问题是什么？我猜这个问题主要是由于FFTW3。也许我没有正确使用它，尤其是在边界条件方面。

Answer 1

当然，有限高斯谱的 FFT 的实部存在负分量。您只是在使用变换的实部。所以你的情节是绝对正确的。

您似乎将实数部分误认为幅度，这当然不会是负数。为此，您需要 fftw_plan_dft_r2c_1d 然后计算复系数的绝对值。或者您可能将傅里叶变换误认为是有限 DFT。

您可能需要检查此处以说服自己上述计算的正确性：

http://docs.mantidproject.org/nightly/algorithms/FFT-v1.html

请记住，上页的图表发生了偏移，因此 0 频率位于频谱的中间。

引用您自己的话，如果将最高频率归一化为 0，[r*exp(-(r^2))*sin(kr)]dr 的数值积分将对所有 k>1 产生负分量。

TLDR：您的情节绝对是最先进的，并且符合离散和有限的功能分析。

Answer 2

查看 FFTW 站点中的相关页面（Real-to-Real Transforms、transform kinds、Real-odd DFT (DST)）和 Fortran 的头文件，似乎 FFTW 期望 FFTW_RODFT00 等而不是 FFTW_FORWARD 指定种类的 真实到真实的转换。例如，

! my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_FORWARD, FFTW_ESTIMATE )
my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_RODFT00, FFTW_ESTIMATE )

执行上页所示的“I 型”离散正弦变换 (DST-I)。这种修改似乎解决了这个问题（即，使傅里叶变换成为具有正值的高斯变换）。

---

以下是 OP 的代码稍作修改的版本，用于试验上述修改：

! ... only the modified part is shown...
real(dp) :: delta, k, r, fftw, num, ana
integer :: i, j, n
type(C_PTR) ::  my_plan
real(C_DOUBLE), allocatable :: y(:), yy(:)

delta = 0.0125_dp ; n = 1024   ! rmax = 12.8
! delta = 0.1_dp    ; n = 128    ! rmax = 12.8
! delta = 0.2_dp    ; n = 64    ! rmax = 12.8
! delta = 0.4_dp    ; n = 32    ! rmax = 12.8

allocate( y( n ), yy( n ) )

! my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_FORWARD, FFTW_ESTIMATE )
my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_RODFT00, FFTW_ESTIMATE )

! Loop over r-grid
do i = 1, n
    r = i * delta              ! (2-a)
    y( i )= r * exp( -r**2 )
end do

call fftw_execute_r2r( my_plan, y, yy )

! Loop over k-grid
do i = 1, n

    ! Result of FFTW
    k = i * pi / ((n + 1) * delta)    ! (2-b)
    fftw = 4 * pi * delta * yy( i ) / k / 2   ! the last 2 due to RODFT00

    ! Numerical result via quadrature
    num = 0
    do j = 1, n
        r = j * delta
        num = num + r * exp( -r**2 ) * sin( k * r )
    enddo
    num = num * 4 * pi * delta / k

    ! Analytical result
    ana = sqrt( pi )**3 * exp( -k**2 / 4 )

    ! Output
    write(10,*) k, fftw
    write(20,*) k, num
    write(30,*) k, ana
end do

编译（使用 gfortran-8.2 + FFTW3.3.8 + OSX10.11）：

$ gfortran -fcheck=all -Wall sine.f90 -I/usr/local/Cellar/fftw/3.3.8/include -L/usr/local/Cellar/fftw/3.3.8/lib -lfftw3

如果我们在原始代码中使用FFTW_FORWARD，我们得到

具有负瓣（其中 fort.10、fort.20 和 fort.30 对应于 FFTW、正交和分析结果）。修改代码以使用FFTW_RODFT00 会更改如下结果，因此修改似乎正在工作（但请参阅下面的网格定义）。

---

补充说明

• 我稍微修改了代码中 r 和 k 的网格定义（第 (2-a) 和 (2-b) 行），发现这可以提高准确性。但是我还是不确定上面的定义是否与FFTW使用的定义相符，所以请阅读手册了解详情...

• fftw3.f03 头文件给出了fftw_plan_r2r_1d 的接口

  type(C_PTR) function fftw_plan_r2r_1d(n,in,out,kind,flags) bind(C, name='fftw_plan_r2r_1d')
    import
    integer(C_INT), value :: n
    real(C_DOUBLE), dimension(*), intent(out) :: in
    real(C_DOUBLE), dimension(*), intent(out) :: out
    integer(C_FFTW_R2R_KIND), value :: kind
    integer(C_INT), value :: flags
  end function fftw_plan_r2r_1d
  

• （由于没有 Tex 支持，这部分非常难看...）4 pi r^2 * exp(-r^2) * sin(kr)/(kr) 对于 r = 0 -> 无限的积分是 pi^(3/2) * exp(-k^2 / 4)（从 Wolfram Alpha 获得或通过注意实际上是 exp(-(x^2 + y^2 + z^2)) 通过 exp(-i*(k1 x + k2 y + k3 z)) 与 k =(k1,k2) 的 3-D 傅里叶变换,k3))。所以，虽然有点反直觉，但结果变成了正高斯。

• 我猜 r-grid 可以选择得更粗一些（例如 delta 高达 0.4），只要它覆盖转换函数的频域（这里是 exp(-r^2)），它就可以提供几乎相同的精度。