python - 在 3-D numpy 数组的小切片上高效使用 1-D pyfftw

Question

我有一个 3D 数据立方体，其大小值约为 10,000x512x512。我想沿着 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 数组；变量channelChunks。我将每个 channel 保存为二进制文件(未包含在上面的代码中，但保存部分工作正常)。

此过程适用于我们的解调项目，然后将 4D 数据立方体 ChannelChunks 解析为 eval(numChannel) 3D 数据立方体(电影)，并且根据我们的实验设置，我们可以根据颜色将电影分开。我希望我可以绕过编写一个通过 pyfftw 调用矩阵上的 fft 的 C++ 函数。

实际上，我在给定索引 1 和 2 轴处沿着 dataChunk 的 0 轴获取 windowSize=6 个元素，并执行 1D FFT。我需要在整个 3D 数据 block 中执行此操作以生成解调的电影。谢谢。

Answer 1

FFTW advanced plans可以由 pyfftw 自动构建。代码可以按以下方式修改:

• 实数到复数变换可以用来代替复数到复数变换。使用 pyfftw，它通常会这样写:

  ftwIn = ftw.empty_aligned(windowSize, dtype='float64')
  ftwOut = ftw.empty_aligned(windowSize//2+1, dtype='complex128')
  fftObject = ftw.FFTW(ftwIn,ftwOut)
  

• 向 FFTW 规划器添加一些标志。例如，FFTW_MEASURE 将对不同的算法进行计时并选择最佳的。 FFTW_DESTROY_INPUT 表示可以修改输入数组:可以使用一些实现技巧。

  fftObject = ftw.FFTW(ftwIn,ftwOut, flags=('FFTW_MEASURE','FFTW_DESTROY_INPUT',))
  

• 限制分区数量。除法比乘法花费更多。

  scale=1.0/windowSize
  for ...
      for ...
          2*np.abs(ftwOut[:,:,:])*scale  #instead of /windowSize
  

• 通过使用 FFTW advanced plan 避免多个 for 循环通过 pyfftw。

  nbwindow=numFrames//windowSize
  # create fftw arrays
  ftwIn = ftw.empty_aligned((nbwindow,windowSize,dataChunk.shape[2]), dtype='float64')
  ftwOut = ftw.empty_aligned((nbwindow,windowSize//2+1,dataChunk.shape[2]), dtype='complex128')
  fftObject = ftw.FFTW(ftwIn,ftwOut, axes=(1,), flags=('FFTW_MEASURE','FFTW_DESTROY_INPUT',))
  
  ...
  for yy in range(data.shape[1]):
      ftwIn[:] = np.reshape(data[0:nbwindow*windowSize,yy,:],(nbwindow,windowSize,data.shape[2]),order='C')
      fftObject()
      channelOut[:,:,yy,:]=np.transpose(2*np.abs(ftwOut[:,:,:])*scale, (1,0,2))
  

这是修改后的代码。我还将帧数减少到 100，设置随机生成器的种子以检查结果是否未修改并评论 tkinter。 窗口的大小可以设置为2的幂，或者是2、3、5或7相乘的数字，从而可以有效地应用Cooley-Tuckey算法。避免使用大质数。

import numpy as np
import pyfftw as ftw
#from tkinter import simpledialog
from math import ceil
import multiprocessing
import time


ftw.config.NUM_THREADS = multiprocessing.cpu_count()
ftw.interfaces.cache.enable()
ftw.config.PLANNER_EFFORT = 'FFTW_MEASURE'

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
    np.random.seed(seed=42)
    dataChunk = np.random.random((100,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?'))
    windowSize=32
    numChannels = windowSize//2+1

    nbwindow=numFrames//windowSize
    # create fftw arrays
    ftwIn = ftw.empty_aligned((nbwindow,windowSize,dataChunk.shape[2]), dtype='float64')
    ftwOut = ftw.empty_aligned((nbwindow,windowSize//2+1,dataChunk.shape[2]), dtype='complex128')

    #ftwIn = ftw.empty_aligned(windowSize, dtype='complex128')
    #ftwOut = ftw.empty_aligned(windowSize, dtype='complex128')
    fftObject = ftw.FFTW(ftwIn,ftwOut, axes=(1,), flags=('FFTW_MEASURE','FFTW_DESTROY_INPUT',))
    # 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
    printed=0
    nbwindow=data.shape[0]//windowSize
    scale=1.0/windowSize
    for yy in range(data.shape[1]):
        #for xx in range(data.shape[2]):
            index = 0

            ftwIn[:] = np.reshape(data[0:nbwindow*windowSize,yy,:],(nbwindow,windowSize,data.shape[2]),order='C')
            fftObject()
            channelOut[:,:,yy,:]=np.transpose(2*np.abs(ftwOut[:,:,:])*scale, (1,0,2))
            #for i in range(nbwindow):
                #channelOut[:,i,yy,xx] = 2*np.abs(ftwOut[i,:])*scale

            if printed==0:
                      for j in range(channelOut.shape[0]):
                          print j,channelOut[j,0,yy,0]
                      printed=1

    return channelOut

if __name__ == '__main__':
    seconds=time.time()
    runme()
    print "time: ", time.time()-seconds

让我们知道它可以加快您的计算速度!我在电脑上的时间从 24 秒缩短到不到 2 秒...