fortran - Fortran 结果中的 FFTW 仅包含零
问题描述
我一直在尝试编写一个简单的程序来使用 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 的标签,所以我希望这个问题不是题外话
解决方案
正如@roygvib 和@Ross 在评论中首先解释的那样,计划子例程会覆盖输入数组,因为它们会使用不同的参数多次尝试转换。我将添加一些实际使用注意事项。
你声称你确实关心性能。那么有两种可能:
- 正如您在代码中显示的那样,您只进行一次转换。那么就没有使用的意义了
FFTW_MEASURE
。计划子程序比实际计划执行子程序慢很多倍。使用FFTW_ESTIMATE
它会快得多。
FFTW_MEASURE 告诉 FFTW通过实际计算几个 FFT并测量它们的执行时间来找到优化的计划。根据您的机器,这可能需要一些时间(通常是几秒钟)。FFTW_MEASURE 是默认计划选项。
FFTW_ESTIMATE 指定,而不是不同算法的实际测量,一个简单的启发式用于快速选择一个(可能是次优的)计划。使用此标志,在规划期间不会覆盖输入/输出数组。
http://www.fftw.org/fftw3_doc/Planner-Flags.html
- 您对不同的数据多次执行相同的转换。然后,您必须在第一次转换之前只进行一次计划,然后重新使用该计划。只需先制定计划,然后才用第一个输入数据填充数组。在每次运输之前制定计划会使程序非常缓慢。
推荐阅读
- python - Python:如何用另一个列表索引一个列表,然后用另一个列表再次索引它?
- python - 如何在 Visual Studio 中使用 python 和 C# 构建应用程序
- ruby-on-rails - 如何在给定的 ruby 场景中使用 gsub?
- nuxt.js - Nuxt.js - ssr,错误在将标头发送到客户端后无法设置标头
- docker - 我通过 docker 在我的本地计算机上尝试了 start auditbeat。但是我从弹性搜索中拒绝连接
- javascript - 使用 app-module-path 库时,vscode 是否可以以非相对方式自动导入路径?
- android - 无论如何要限制 Android 设备策略中的设备注册?
- java - 为什么 Wifi 直接发现在 Android 10 上失败?
- web-scraping - Scrapy splash 不加载内容
- python - 如何忽略第一个屈服值?