提问人：bproxauf 提问时间：10/16/2016 最后编辑：ROMANIA_engineerbproxauf 更新时间：5/4/2019 访问量：5053
相对较慢的 python numpy 3D 傅里叶变换
Comparatively slow python numpy 3D Fourier Transformation
问：
对于我的工作，我需要对大图像执行离散傅里叶变换 （DFT）。在当前示例中，我需要 1921 x 512 x 512 图像的 3D FT（以及 512 x 512 图像的 2D FFT）。现在，我正在使用 numpy 包和关联的函数 np.fft.fftn（）。下面的代码片段示例性地显示了大小相等/稍小的 2D/3D 随机数生成网格上的 2D 和 3D FFT 时间，具体方式如下：
import sys
import numpy as np
import time
tas = time.time()
a = np.random.rand(512, 512)
tab = time.time()
b = np.random.rand(100, 512, 512)
tbfa = time.time()
fa = np.fft.fft2(a)
tfafb = time.time()
fb = np.fft.fftn(b)
tfbe = time.time()
print "initializing 512 x 512 grid:", tab - tas
print "initializing 100 x 512 x 512 grid:", tbfa - tab
print "2D FFT on 512 x 512 grid:", tfafb - tbfa
print "3D FFT on 100 x 512 x 512 grid:", tfbe - tfafb
输出：
initializing 512 x 512 grid: 0.00305700302124
initializing 100 x 512 x 512 grid: 0.301637887955
2D FFT on 512 x 512 grid: 0.0122730731964
3D FFT on 100 x 512 x 512 grid: 3.88418793678
我遇到的问题是我经常需要这个过程，所以每张图像花费的时间应该很短。在我自己的计算机上进行测试时（中间段笔记本电脑，分配给虚拟机的 2GB RAM（-->因此测试网格较小）），如您所见，3D FFT 需要 ~ 5 秒（数量级）。现在，在工作中，机器变得更好了，集群/网格架构系统和 FFT 速度要快得多。在这两种情况下，2D 的都是准瞬间完成的。
但是，对于 1921x512x512，np.fft.fftn（） 需要 ~ 5 分钟。由于我猜 scipy 的实现速度并不快，并且考虑到在相同大小网格的 MATLAB FFT 上在 ~ 5 秒内完成，我的问题是是否有一种方法可以将该过程加速到或几乎达到 MATLAB 时间。我对 FFT 的了解有限，但显然 MATLAB 使用 FFTW 算法，而 python 没有。使用一些 pyFFTW 包，我得到类似的时间有什么合理的机会吗？此外，1921 年似乎是一个不吉利的选择，只有 2 个质因数（17、113），所以我认为这也起了作用。另一方面，512 是非常适合 2 的幂。如果可能的话，是否可以在不用零填充到 2048 的情况下实现类似 MATLAB 的时间？
我之所以这么问，是因为我必须大量使用 FFT（这种差异将产生巨大影响！），并且如果无法减少 python 中的计算时间，我将不得不切换到其他更快的实现。
答：
4赞
francis
10/16/2016
#1
是的，与 或 相比，通过接口使用 FFTW 可能会减少您的计算时间。这些 DFT 算法实现的性能可以在基准测试中进行比较，例如： 在 Python 中提高 FFT 性能中报告了一些有趣的结果pyfftwnumpy.fftscipy.fftpack
我建议使用以下代码进行测试：
import pyfftw
import numpy
import time
import scipy
f = pyfftw.n_byte_align_empty((127,512,512),16, dtype='complex128')
#f = pyfftw.empty_aligned((33,128,128), dtype='complex128', n=16)
f[:] = numpy.random.randn(*f.shape)
# first call requires more time for plan creation
# by default, pyfftw use FFTW_MEASURE for the plan creation, which means that many 3D dft are computed so as to choose the fastest algorithm.
fftf=pyfftw.interfaces.numpy_fft.fftn(f)
#help(pyfftw.interfaces)
tas = time.time()
fftf=pyfftw.interfaces.numpy_fft.fftn(f) # here the plan is applied, nothing else.
tas = time.time()-tas
print "3D FFT, pyfftw:", tas
f = pyfftw.n_byte_align_empty((127,512,512),16, dtype='complex128')
#f = pyfftw.empty_aligned((33,128,128), dtype='complex128', n=16)
f[:] = numpy.random.randn(*f.shape)
tas = time.time()
fftf=numpy.fft.fftn(f)
tas = time.time()-tas
print "3D FFT, numpy:", tas
tas = time.time()
fftf=scipy.fftpack.fftn(f)
tas = time.time()-tas
print "3D FFT, scipy/fftpack:", tas
# first call requires more time for plan creation
# by default, pyfftw use FFTW_MEASURE for the plan creation, which means that many 3D dft are computed so as to choose the fastest algorithm.
f = pyfftw.n_byte_align_empty((128,512,512),16, dtype='complex128')
fftf=pyfftw.interfaces.numpy_fft.fftn(f)
tas = time.time()
fftf=pyfftw.interfaces.numpy_fft.fftn(f) # here the plan is applied, nothing else.
tas = time.time()-tas
print "3D padded FFT, pyfftw:", tas
对于 127*512*512 的尺寸，在我普通的计算机上，我得到了：
3D FFT, pyfftw: 3.94130897522
3D FFT, numpy: 16.0487070084
3D FFT, scipy/fftpack: 19.001199007
3D padded FFT, pyfftw: 2.55221295357
所以比 和 快得多。使用填充甚至更快，但计算的内容不同。pyfftwnumpy.fftscipy.fftpack
最后，由于它根据文档使用标志，因此在第一次运行时可能看起来较慢。当且仅当连续计算许多相同大小的 DFT 时，这是一件好事。pyfftwFFTW_MEASURE
评论
0赞
bproxauf
10/18/2016
好的，首先感谢您的回答。作为我工作的一部分，我需要做一个方位角平均，为此，我将两个维度为 1921x512x512 的立方体相乘。起初大约花了 25 秒（太长了，因为我必须经常这样做）。我发现这与大步有关（直到今天我才知道）。numpy FFT 会自动将其从 C 样式更改为 Fortran 样式。有什么方法可以防止这种情况发生（副本除外）？使用相同的 （C） 样式的步幅，时间缩短到 ~ 4 秒。
0赞
bproxauf
10/18/2016
将 axes 参数指定为 （2,1,0） 而不是 （0,1,2） 可以保留步幅顺序，但应该有一种比该解决方法更简单的方法......
0赞
francis
10/19/2016
我不确定是否理解您所说的“numpy FFT 自动将其从 C 更改为 Fortran 样式”是什么意思。您可以使用它来检查尺寸是否反转：事实并非如此。实际上，如果输入的形状为 127x512x512，则输出的形状为 127x512x512。此外，我还定时执行逐元素乘法：它比大小为 127x512x512 的 pyfftw dft 快约 10 倍。因此，如果瓶颈被证明是元素乘法，我会感到惊讶！print fftf.shapenumpy.multiply(f,fftf)
0赞
bproxauf
10/19/2016
我指的是改变的步伐：另请参阅我的新问题：
0赞
bproxauf
10/19/2016
stackoverflow.com/questions/40109915/......我收到了几个答案，说我可以使用 scipy.fftpack 而不是 numpy.fft（我还不知道 pyFFTW，因为我必须等到星期一才能集中安装软件包（我没有 sudo 权限））。显然，那里保留了步幅结构。但我仍然看不出 numpy.fft.fftn 首先改变结构的原因。
0赞
marcin
5/4/2019
#2
您可以尝试英特尔 MKL（数学内核库）的 FFT，它比 FFTW 更快。 英特尔为 Python 提供了 mkl-fft，它取代了 numpy.fft。您需要做的就是键入：
pip install mkl-fft
然后再次运行您的程序，无需任何更改。
此外，numpy 1.17（即将发布）将有新的FFT实现：
用 pocketfft 库替换基于 fftpack 的 FFT 模块
两种实现具有相同的祖先（Paul 的 Fortran77 FFTPACK N. Swarztrauber），但 pocketfft 包含其他修改 在某些情况下，这提高了准确性和性能。为 包含大质因数的 FFT 长度，pocketfft 使用 Bluestein 的 算法，它保持 O（N log N） 运行时复杂度，而不是 对于素长度，向 O（N*N） 恶化。此外，精度 具有接近素数长度的实值 FFT 有所改善，并且与此相当 具有复值 FFT。
评论