【发布时间】：2019-08-29 19:45:10
【问题描述】：
我有一个大小约为 10,000x512x512 的 3D 数据立方体。我想重复解析沿 dim[0] 的向量窗口（例如 6）并有效地生成傅立叶变换。我想我正在将数组复制到 pyfftw 包中，这给了我巨大的开销。我现在正在查看文档，因为我认为我需要设置一个选项，但我可以在语法上使用一些额外的帮助。
这段代码最初是由另一个人用 numpy.fft.rfft 编写的，并用 numba 加速。但是该实现在我的工作站上不起作用，所以我重新编写了所有内容并选择使用 pyfftw。
import numpy as np
import pyfftw as ftw
from tkinter import simpledialog
from math import ceil
import multiprocessing
ftw.config.NUM_THREADS = multiprocessing.cpu_count()
ftw.interfaces.cache.enable()
def runme():
# normally I would load a file, but for Stack Overflow, I'm just going to generate a 3D data cube so I'll delete references to the binary saving/loading functions:
# load the file
dataChunk = np.random.random((1000,512,512))
numFrames = dataChunk.shape[0]
# select the window size
windowSize = int(simpledialog.askstring('Window Size',
'How many frames to demodulate a single time point?'))
numChannels = windowSize//2+1
# create fftw arrays
ftwIn = ftw.empty_aligned(windowSize, dtype='complex128')
ftwOut = ftw.empty_aligned(windowSize, dtype='complex128')
fftObject = ftw.FFTW(ftwIn,ftwOut)
# perform DFT on the data chunk
demodFrames = dataChunk.shape[0]//windowSize
channelChunks = np.zeros([numChannels,demodFrames,
dataChunk.shape[1],dataChunk.shape[2]])
channelChunks = getDFT(dataChunk,channelChunks,
ftwIn,ftwOut,fftObject,windowSize,numChannels)
return channelChunks
def getDFT(data,channelOut,ftwIn,ftwOut,fftObject,
windowSize,numChannels):
frameLen = data.shape[0]
demodFrames = frameLen//windowSize
for yy in range(data.shape[1]):
for xx in range(data.shape[2]):
index = 0
for i in range(0,frameLen-windowSize+1,windowSize):
ftwIn[:] = data[i:i+windowSize,yy,xx]
fftObject()
channelOut[:,index,yy,xx] = 2*np.abs(ftwOut[:numChannels])/windowSize
index+=1
return channelOut
if __name__ == '__main__':
runme()
我得到了一个 4D 数组；变量通道块。我将每个通道保存到二进制文件中（上面的代码中不包含，但保存部分工作正常）。
这个过程适用于我们的解调项目，然后将 4D 数据立方体 channelChunks 解析为 eval(numChannel) 3D 数据立方体（电影），并且根据我们的实验设置，我们能够通过颜色分离电影。我希望我可以绕过编写一个通过 pyfftw 调用矩阵上的 fft 的 C++ 函数。
实际上，我在给定的 1 轴和 2 轴索引处沿 dataChunk 的 0 轴获取 windowSize=6 个元素并执行 1D FFT。我需要在整个 3D 数据块中执行此操作以生成解调电影。谢谢。
【问题讨论】：
-
您对此进行了分析吗？这应该让您很好地了解实际花费的时间。在大多数情况下，编写 C++ 函数不会加快速度。
-
是的，我运行了 cProfiler。大约 98% 的时间都花在了 getDFT 调用上。问题似乎出在 for 循环中：for i in range(0,frameLen-windowSize+1,windowSize): 但是，分析器没有详细说明，所以我需要获得更好的统计数据，也许可以调用另一个选项。无论如何，我认为问题出在 ftwIn[:] = data[i:i+windowSize,yy,xx] 行的 for 循环中。在axis0中10,000帧需要40分钟，我很肯定如果内存管理得当，我可以得到1-2个数量级的改进。
-
cProfiler.run应该提供每个函数调用所花费的时间——这只是了解那里有什么的问题。函数getDFT是最明显的，因为它循环，与其他所有函数相比，它需要很长时间。您可以在 fft 函数中用axis参数替换xx和yy上的循环。
标签： python numpy numpy-slicing pyfftw