FFTW：真实到复杂和复杂到真实 2D 变换的问题答案

【问题标题】：FFTW: Trouble with real to complex and complex to real 2D tranfsormsFFTW：真实到复杂和复杂到真实 2D 变换的问题
【发布时间】：2011-09-22 18:14:25
【问题描述】：

正如标题所述，我正在使用带有 Fortran 90/95 的 FFTW（版本 3.2.2）来执行真实数据的 2D FFT（实际上是一个随机数字段）。我认为向前迈出了一步（至少我得到了一些输出）。但是我想通过 IFFT 检查所有内容，看看我是否可以重新构建原始输入。不幸的是，当我将复杂程序调用为真正的例程时，什么也没有发生，也没有得到错误输出，所以我有点困惑。下面是一些sn-ps的代码：

implicit none

include "fftw3.f"

! - im=501, jm=401, and lm=60

real*8    :: u(im,jm,lm),recov(im,jm,lm)
complex*8 :: cu(1+im/2,jm)
integer*8 :: planf,planb    
real*8    :: dv

! - Generate array of random numbers
dv=4.0
call random_number(u)
u=u*dv
recov=0.0

k=30

! - Forward step (FFT)

call dfftw_plan_dft_r2c_2d(planf,im,jm,u(:,:,k),cu,FFTW_ESTIMATE)
call dfftw_execute_dft_r2c(planf,u(:,:,k),cu)
call dfftw_destroy_plan(planf)

! - Backward step (IFFT)

call dfftw_plan_dft_c2r_2d(planb,im,jm,cu,recov(:,:,k),FFTW_ESTIMATE)
call dfftw_execute_dft_c2r(planb,cu,recov(:,:,k))
call dfftw_destroy_plan(planb)

上面的前进步骤似乎有效（r2c），但后退步骤似乎无效。我通过区分 u 和 recov 数组来检查这一点 - 它最终不为零。此外，recov 数组的最大值和最小值都为零，这似乎表明没有任何变化。

我查看了 FFTW 文档，并将我的实现基于以下页面 http://www.fftw.org/fftw3_doc/Fortran-Examples.html#Fortran-Examples 。我想知道这个问题是否与索引有关，至少这是我倾向于的方向。无论如何，如果有人可以提供一些帮助，那就太好了！

谢谢！

【问题讨论】：

需要注意的一点是，当您执行 FFT 后跟 IFFT 时，结果将按相对于原始输入的常数因子进行缩放。比例因子的值取决于 FFT 和 IFFT 的定义/实现方式，但通常为 N，其中 N 是 FFT 中的点数。当然你这里也可能有其他不相关的问题。
我发现了这个问题 - 它与 FFT 调用完全无关。在我的代码（不是上面的示例）中，我将 k 设为一维数组，但一直将其作为标量索引值进行访问。上面的代码应该没有什么问题。哎呀！就这样这对其他人仍然有用，我最终在网上找到了一些代码，其中包含许多使用 Fortran 调用 FFTW 例程的非常有用的示例：people.sc.fsu.edu/~jburkardt/f_src/fftw3/fftw3_prb.f90
@Jen 在上面的代码中，您通过混合 real*8 和 complex*8 变量而失去精度：complex*8 对应于 real*4，而不是 real*8。请改用 complex*16。
仁，实际上我什至无法让该示例代码工作。您是否在开始时启用了任何东西以使其与 gfortran 一起使用？他还调用了当前版本中不存在的文件，因为只有 .f 包装器而不是 .f90。
@Thomas James，我没有尝试使用 gfortran 的示例 - 我使用的是 xlf。你确定你的编译步骤有正确的命令/语法吗？ .f 包装器应该没问题（替换 .f90）。如果您愿意，我可以将用于编译代码的 Makefile 发送给您——请告诉我。我尝试将其粘贴在这里，但它不适合作为“评论”。

标签： fortran fft fortran90 fftw ifft

【解决方案1】：

不确定这是否是所有问题的根源，但你声明变量的方式可能是罪魁祸首。

对于大多数编译器（这显然甚至不是标准），Complex*8 是单精度的旧语法：复数变量共占用 8 个字节，在实部和虚部之间共享（4+4 字节)。

[在 Vladimir F 对我的回答发表评论后编辑 1，有关详细信息，请参阅他的链接：] 根据我的经验（即我曾经使用过的系统/编译器），Complex(Kind=8) 对应于双精度复数的声明（a实部和虚部，均占 8 个字节）。

在任何系统/编译器上，Complex(Kind=Kind(0.d0)) 应该声明一个双精度复数。

简而言之，您的复杂数组的大小不合适。将出现的Real*8 和Complex*8 分别替换为Real(kind=8) 和Complex(Kind=8)（或Complex(Kind=kind(0.d0)) 以获得更好的可移植性）。

【讨论】：

你是对的，这是问题所在。但是您不正确地建议 kind=8 是一件好事并且它总是意味着 8 字节组件。见stackoverflow.com/questions/838310/fortran-90-kind-parameter 一些编译器会拒绝编译 kind=8。
谢谢。我的回答确实充满了不准确。希望这会更好。
好多了。