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

Question

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

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

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

的积分

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

然而，结果并不是动量的高斯函数space。

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 频率位于频谱的中间。

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

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 )

执行上页中显示的 "type-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 更改结果如下，因此修改似乎有效（但请参阅下面的网格定义）。

补充说明

我在我的代码（第 (2-a) 和 (2-b) 行）中稍微修改了 r 和 k 的网格定义，发现这提高了准确性。但是我还是不确定上面的定义是否符合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) for r = 0 -> infinite 的积分是pi^(3/2) * exp(-k^2 / 4)（从Wolfram Alpha 或者注意到这实际上是 exp(-(x^2 + y^2 + z^2)) 的 3-D 傅里叶变换 exp(-i*(k1 x + k2 y + k3 z)) 和 k =(k1,k2,k3))。所以，虽然有点违反直觉，但结果变成了正高斯分布。
我想 r-grid 可以选择得更粗糙（例如 delta 高达 0.4），只要它覆盖变换函数的频域（这里 exp(-r^2)).

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

Evaluating the fast Fourier transform of Gaussian function in FORTRAN using FFTW3 library

fortran

physics

fft

fftw