Fortran 结果中的 FFTW 仅包含零答案

【问题标题】：FFTW in Fortran result contains only zerosFortran 结果中的 FFTW 仅包含零
【发布时间】：2018-09-07 13:52:04
【问题描述】：

我一直在尝试编写一个简单的程序来使用 fftw3 对一维输入数组执行 fft。在这里，我使用地震图作为输入。然而，输出数组只包含零。

我知道输入是正确的，因为我也尝试在 MATLAB 中对同一个输入文件进行 fft，这给出了正确的结果。没有编译错误。我正在使用 f95 来编译它，但是，gfortran 也给出了几乎相同的结果。这是我写的代码：-

program fft

    use functions
    implicit none
    include 'fftw3.f90'
    integer nl,row,col
    double precision, allocatable :: data(:,:),time(:),amplitude(:)
    double complex, allocatable :: out(:)
    integer*8 plan


    open(1,file='test-seismogram.xy')
    nl=nlines(1,'test-seismogram.xy')
    allocate(data(nl,2))
    allocate(time(nl))
    allocate(amplitude(nl))
    allocate(out(nl/2+1))
    do row = 1,nl
        read(1,*,end=101) data(row,1),data(row,2)
        amplitude(row)=data(row,2)
    end do
    101 close(1)


    call dfftw_plan_dft_r2c_1d(plan,nl,amplitude,out,FFTW_R2HC,FFTW_PATIENT)
    call dfftw_execute_dft_r2c(plan, amplitude, out)
    call dfftw_destroy_plan(plan)


    do row=1,(nl/2+1)
        print *,out(row)
    end do


    deallocate(data)
    deallocate(amplitude)
    deallocate(time)
    deallocate(out)
end program fft

nlines() 函数是一个用于计算文件行数的函数，它可以正常工作。它在称为函数的模块中定义。

这个程序几乎试图效仿http://www.fftw.org/fftw3_doc/Fortran-Examples.html的例子

我可能只是犯了一个非常简单的逻辑错误，但我严重无法弄清楚这里出了什么问题。任何指针都会非常有帮助。

这几乎是整个输出的样子：-

           .
           .
           .
           (0.0000000000000000,0.0000000000000000)
           (0.0000000000000000,0.0000000000000000)
           (0.0000000000000000,0.0000000000000000)
           (0.0000000000000000,0.0000000000000000)
           (0.0000000000000000,0.0000000000000000)
           .
           .
           .

我的疑问是直接关于 fftw，因为 SO 上有 fftw 的标签，所以我希望这个问题不是题外话

【问题讨论】：

在 2018 年，我强烈建议使用现代 (Fortran 2003) 界面。您是否也打印了输入数据？从文件中读取的数据是否正确？我更强烈建议使用大于 1 的单元号和action="read", status="old"。
是的，我确实打印了输入数据，并且读取正确。关于现代 Fortran，由于我实验室中的大部分代码都是用旧版 Fortran 编写的，所以我什至必须学习一些 f77 语法。不过，我正在慢慢适应现代界面。
我会记住您的建议。谢谢你。但是，我认为这不是解决这个问题的方法吗？
仅与链接页面相比，是否需要FFTW_R2HC？（链接页面中不存在）
调用 dfftw_plan_dft_r2c_1d() 时，输入和输出数组是否可能被覆盖以测量性能？（根据this page）。 This 和 this 页面也可能是相关的（稍后）

标签： fortran fft gfortran fftw

【解决方案1】：

正如 @roygvib 和 @Ross 在 cmets 中首先解释的那样，计划子例程会覆盖输入数组，因为它们会使用不同的参数多次尝试转换。我将添加一些实际使用注意事项。

您声称您确实关心性能。那么有两种可能：

您只进行一次转换，正如您在代码中显示的那样。那么使用FFTW_MEASURE就没有意义了。计划子例程比实际计划执行子例程慢许多倍。使用FFTW_ESTIMATE 会快很多。

FFTW_MEASURE 告诉 FFTW 找到一个优化的计划 实际上计算几个 FFT 并测量它们的执行时间。取决于在您的机器上，这可能需要一些时间（通常是几秒钟）。 FFTW_MEASURE 是默认的计划选项。

FFTW_ESTIMATE 指定，而不是实际测量不同的算法，一个简单的启发式用于选择一个（可能次优）计划快速。使用这个标志，输入/输出数组是在规划期间不会被覆盖。

http://www.fftw.org/fftw3_doc/Planner-Flags.html

您多次对不同的数据执行相同的转换。然后，您必须在第一次转换之前只进行一次计划，然后重新使用该计划。只需先制定计划，然后才用第一个输入数据填充数组。在每次运输之前制定计划会使计划变得非常缓慢。

【讨论】：