在 Fortran 混乱输出中使用 FFTW 对正弦 (x) 进行 DFT

Question

以下是我编写的用于查找一段时间内正弦 (x) 的 DFT 的代码。

program fftw_test

implicit none

INTEGER FFTW_MEASURE
    PARAMETER (FFTW_MEASURE=0)
INTEGER FFTW_ESTIMATE
    PARAMETER (FFTW_ESTIMATE=64)
INTEGER FFTW_FORWARD
    PARAMETER (FFTW_FORWARD=-1)

integer, parameter :: n = 8
integer :: i
double complex, dimension(0:n-1) :: input, output
double precision, parameter :: pi = 3.141592653, h = 2.0d0*pi/(n)
integer*8 :: plan

call dfftw_plan_dft_1d(plan, n, input, output, fftw_forward, fftw_measure)

do i = 0, n-1
    input(i) = cmplx(sin(h*i), 0)
end do

call dfftw_execute_dft(plan, input, output)

output = output/n 

output(0) = cmplx(0,0) ! setting oddball wavenumber to be 0

call dfftw_destroy_plan(plan)

do i = -n/2, n/2-1, 1
    write(*, *) i, output(i+(n/2))
end do

end program

我知道 FFTW 库中的 r2c（实数到复数）函数。但有人建议我使用普通的 c2c 功能。所以我将输入函数定义为实部 = sine(x) 和复数部分 0 的复数。

sine(x) 的 DFT 应该是 fk(-1) = cmplx(0, -0.5) 和 fk(1) = cmplx(0, 0.5) 其中 fk(k) 表示傅里叶系数第 k 个波数

我收到的输出如下。

      -4 (  0.0000000000000000     ,  0.0000000000000000     )
      -3 (  3.2001271327131496E-008,-0.49999998518472011     )
      -2 ( -1.0927847071684482E-008,  1.4901161193847656E-008)
      -1 ( -1.0145577183762535E-008,  1.4815279864022202E-008)
       0 ( -1.0927847071684482E-008,  0.0000000000000000     )
       1 ( -1.0145577183762535E-008, -1.4815279864022202E-008)
       2 ( -1.0927847071684482E-008, -1.4901161193847656E-008)
       3 (  3.2001271327131496E-008, 0.49999998518472011     )

我得到 fk(-3) = cmplx(~0, -0.5) 和 fk(3) = cmplx(~0, 0.5)。如果我将网格大小增加到 16、32 左右，我会得到具有所需值的 -n/2 -1 和 n/2 -1 波数，而不是 -1 和 1 波数。

这与 FFTW 在输出数组中存储输出的方式有关吗？还是我哪里出错了？

此外，我似乎没有 'proper 0' 我应该去的地方。相反，它是 10^(-8) 数量级的数字，我认为这是我的数据类型 double 可以容纳的最小数字。这是我应该担心的事情吗？

Answer 1

就像@VladimirF 已经说过的，值的顺序与您预期的有点不同。数组的前半部分包含正频率，后半部分以相反的顺序包含负频率（参见 this link）。您可能需要检查 FFTW 使用的符号约定。

准确性问题源于 pi 的单精度值和使用 cmplx 生成单精度复数（使用关键字参数 kind）。在这种情况下，您可以简单地将您的实际值分配给复杂变量。应用这两个更改会产生 ~1e-10 的精度。这可以通过为 pi 提供更好的近似值（即超过 10 位数字）来改进。

例如值 pi = 3.141592653589793d0 产生的结果精度为 1e-16.

在 Fortran 混乱输出中使用 FFTW 对正弦 (x) 进行 DFT

DFT of sine(x) using FFTW in Fortran muddled output

fortran

fft

fftw

dft