1

我正在尝试使用 FFT 和 pyfftw 实现 3d 卷积。我使用在 SO 的另一篇文章中发布的代码作为基础：

class CustomFFTConvolution(object):

def __init__(self, A, B, threads=1):

    shape = (np.array(A.shape) + np.array(B.shape))-1
    #shape=np.array(A.shape) - np.array(B.shape)+1
    if np.iscomplexobj(A) and np.iscomplexobj(B):
        self.fft_A_obj = pyfftw.builders.fftn(
                A, s=shape, threads=threads)
        self.fft_B_obj = pyfftw.builders.fftn(
                B, s=shape, threads=threads)
        self.ifft_obj = pyfftw.builders.ifftn(
                self.fft_A_obj.get_output_array(), s=shape,
                threads=threads)

    else:
        self.fft_A_obj = pyfftw.builders.rfftn(
                A, s=shape, threads=threads)
        self.fft_B_obj = pyfftw.builders.rfftn(
                B, s=shape, threads=threads)
        self.ifft_obj = pyfftw.builders.irfftn(
                self.fft_A_obj.get_output_array(), s=shape,
                threads=threads)

def __call__(self, A, B):
    s1=np.array(A.shape)
    s2=np.array(B.shape)

    fft_padded_A = self.fft_A_obj(A)
    fft_padded_B = self.fft_B_obj(B)

    ret= self.ifft_obj(fft_padded_A * fft_padded_B)

    return self._centered(ret, s1 - s2 + 1)

def _centered(self,arr, newshape):
    # Return the center newshape portion of the array.
    newshape = np.asarray(newshape)
    currshape = np.array(arr.shape)
    startind = (currshape - newshape) // 2
    endind = startind + newshape
    myslice = [slice(startind[k], endind[k]) for k in range(len(endind))]
    return arr[tuple(myslice)]

我的数据 A 的形状为 (931, 411, 806)，我的过滤器 B 的形状为 (32, 32, 32)。如果我在 24 核机器上使用 24 个线程运行此代码，则该操作需要 263 秒。现在，如果我在同一台机器上运行相同的实验，但这次 A 的形状为 (806, 411, 931)只是轴的交换，代码只需要 16 秒。这是什么原因？是否有获得最佳性能的经验法则？也许填充尺寸之一？谢谢！

4

1 回答 1

1

既然考虑了填充，那么填充的大小可以增加到偶数，还是小素数的倍数？选择均匀尺寸可以将挂钟时间除以 3。

根据维度，某些 DFT 算法可能不可用或效率不高。例如，执行 DFT 的最有效算法之一是Cooley-Tukey 算法。它在于将复合大小 N=N1*N2 的信号的 DFT 划分为大小为 N2 的 N1 个 DTF。因此，它对于通过乘以小素因数（2、3、5、7）获得的复合尺寸效果更好， FFTW 中为此提供了专门的高效算法。从FFTW 的文档中：

例如，标准 FFTW 分布对大小可以分解为小素数（2、3、5 和 7）的数组最有效，否则它使用较慢的通用例程。如果您需要其他大小的有效转换，您可以使用 FFTW 的代码生成器，它为您可能关心的任何特定数组大小生成快速 C 程序（“codelet”）。例如，如果您需要大小为 513 = 19*33 的变换，则可以自定义 FFTW 以有效地支持因子 19。

您的填充尺寸具有较高的主要因素：

931=>962=2*13*37
411=>442=2*13*17
806=>837=3*3*3*31

可以扩展填充以更接近具有小质数的数字，例如 980、448 和 864。然而，填充 3D 图像会导致内存占用显着增加，以至于它并不总是可能的。

为什么改变维度的顺序会改变计算时间？差异可能是由于输入数组是真实的。因此，在一个维度上执行 R2C DFT，然后在第二个和第三个维度上执行 C2C 以计算 3D DFT。如果要变换的第一个维度的大小是偶数，则可以将 R2C 变换变成一半大小的复 DFT，如下所示。这个技巧不适用于奇数大小。因此，一些快速算法可能会随着 962 和 837 的翻转而变得可用。

这是一个测试它的代码：

import pyfftw
import matplotlib.pyplot as plt
import multiprocessing
import numpy as np
from timeit import default_timer as timer

def listofgoodsizes():
    listt=[]
    p2=2
    for i2 in range(11):
        p3=1
        for i3 in range(7):
            p5=1
            for i5 in range(2):

                listt.append(p2*p3*p5)
                p5*=5
            p7=1
            for i7 in range(2):
                listt.append(p2*p3*p7)
                p7*=7

            p3*=3
        p2*=2
    listt.sort()
    return listt

def getgoodfftwsize(n,listt):
    for i in range(len(listt)):
        if listt[i]>=n:
            return listt[i]
    return n

def timea3DR2CDFT(n,m,p):
    bb = pyfftw.empty_aligned((n,m, p), dtype='float64')
    bf= pyfftw.empty_aligned((n,m, (p/2+1)), dtype='complex128')
    pyfftw.config.NUM_THREADS = 1 #multiprocessing.cpu_count()
    fft_object_b = pyfftw.FFTW(bb, bf,axes=(0,1,2))

    print n,m,p
    start = timer()
    fft_object_b(bb)
    end = timer()
    print end - start

#three prime numbers !      
n=3*37
m=241
p=5*19

timea3DR2CDFT(n,m,p)



# to even size :
neven=2*((n+1)/2)
meven=2*((m+1)/2)
peven=2*((p+1)/2)

timea3DR2CDFT(neven,meven,peven)


#to nearest multiple of prime
listt=listofgoodsizes()

ngood=getgoodfftwsize(n,listt)
mgood=getgoodfftwsize(m,listt)
pgood=getgoodfftwsize(p,listt)

timea3DR2CDFT(ngood,mgood,pgood)

在我的电脑上，它打印：

111 241 95
0.180601119995
112 242 96
0.0560319423676
112 252 96
0.0564918518066
于 2019-04-18T21:41:17.640 回答