计算卷积的最快方法答案

【问题标题】：Fastest method for calculating convolution计算卷积的最快方法
【发布时间】：2013-12-31 12:33:24
【问题描述】：

有人知道计算卷积的最快方法吗？不幸的是，我处理的矩阵非常大（500x500x200），如果我在 MATLAB 中使用 convn 需要很长时间（我必须在嵌套循环中迭代这个计算）。因此，我将卷积与 FFT 结合使用，现在速度更快了。但是，我仍在寻找一种更快的方法。有什么想法吗？

【问题讨论】：

CUFFT 非常好，但可能无法做一个不是 2 的幂对齐的矩阵。此外，您还需要硬件并知道自己在做什么。

标签： c++ matlab signal-processing convolution template-matching

【解决方案1】：

如果您的内核是可分离的，那么通过执行多个连续的一维卷积将实现最大的速度增益。

MathWorks 的 Steve Eddins 在his blog 上描述了当内核在 MATLAB 上下文中可分离时如何利用卷积的关联性来加速卷积。对于P-by-Q 内核，与 2D 卷积相比，执行两个单独和顺序卷积的计算优势是 PQ/(P+Q)，对应于 9x9 内核的 4.5x 和 15x15 内核的 ~11x。编辑：在this Q&A 中给出了一个有趣的不知情的证明。

要确定内核是否可分离（即两个向量的外积），请参阅博客 goes on to describe 如何检查您的内核是否可与 SVD 分离以及如何获得一维内核。他们的示例是针对 2D 内核的。 N维可分离卷积的解决方案，请查看this FEX submission。

另一个值得指出的资源是this SIMD (SSE3/SSE4) implementation of 3D convolution by Intel，它包括source 和presentation。该代码适用于 16 位整数。除非您转向 GPU（例如 cuFFT），否则可能很难比 Intel 的实现更快，其中还包括 Intel MKL。 this page of the MKL documentation 底部有一个 3D 卷积（单精度浮点）示例（链接已修复，现在在 https://stackoverflow.com/a/27074295/2778484 中镜像）。

【讨论】：

有趣的是，imfilter 函数实际上隐含地执行了此操作。它需要一个二维数组作为内核，但在应用过滤器之前会检查它是否可分离。此外，如前所述，如果您进行循环卷积，FFT 也会很快。
@jucestain 这是一个很好的观点。 I have noted this before 是 imfilter 在循环中调用时更快的原因，如果您尝试过滤每个具有相同 2D 内核的 2D 图像堆栈，而不是给它提供图像堆栈，即使它支持这样做。如果它检测到 3D 数据，即使 2D 内核是可分离的（功能或错误？），它也会声明内核不可分离。
不幸的是我的矩阵似乎是不可分离的！并且无法使用该功能。
@Nicole 内核是什么？数据矩阵不必是任何特别的东西，只要是内核。

【解决方案2】：

您可以尝试重叠添加和重叠保存方法。它们涉及将您的输入信号分解成更小的块，然后使用上述任何一种方法。

FFT 最有可能（我可能错了）是最快的方法，尤其是当您使用 MATLAB 中的内置例程或 C++ 中的库时。除此之外，将输入信号分成更小的块应该是一个不错的选择。

【讨论】：

由于我想用卷积来进行模式匹配，所以我认为分解矩阵是有问题的！
@Nicole 如果您可以使用信号处理工具箱，fftfilt 应该可以为您完成繁重的工作。 mathworks.de/de/help/signal/ref/fftfilt.html
@blackbird：但是我如何在 matlab 中使用该命令而不是 convn 呢？假设我有 a = rand(500,500,100) 和 b = rand(20,20,20)
@Nicole 对不起，我没有注意到您要计算三维卷积。我必须传递那个。

【解决方案3】：

我有两种方法来计算 fastconv

2 比 1 好

1- 犰狳您可以使用 armadillo 库通过此代码计算 conv

cx_vec signal(1024,fill::randn);
cx_vec code(300,fill::randn);
cx_vec ans = conv(signal,code);

2-使用 fftw ans sigpack 和 armadillo 库以这种方式计算快速转换，您必须在构造函数中初始化代码的 fft

FastConvolution::FastConvolution(cx_vec inpCode)
{
    filterCode = inpCode;
    fft_w = NULL;
}


cx_vec FastConvolution::filter(cx_vec inpData)
{
int length = inpData.size()+filterCode.size();
    if((length & (length - 1)) == 0)
    {

    }
    else
    {
        length = pow(2 , (int)log2(length) + 1);
    }
    if(length != fftCode.size())
        initCode(length);

    static cx_vec zeroPadedData;
    if(length!= zeroPadedData.size())
    {
        zeroPadedData.resize(length);
    }
    zeroPadedData.fill(0);
    zeroPadedData.subvec(0,inpData.size()-1) = inpData;


    cx_vec fftSignal = fft_w->fft_cx(zeroPadedData);
    cx_vec mullAns = fftSignal % fftCode;
    cx_vec ans = fft_w->ifft_cx(mullAns);
    return ans.subvec(filterCode.size(),inpData.size()+filterCode.size()-1);
}

void FastConvolution::initCode(int length)
{
    if(fft_w != NULL)
    {
        delete fft_w;
    }
    fft_w = new sp::FFTW(length,FFTW_ESTIMATE);
    cx_vec conjCode(length,fill::zeros);
    fftCode.resize(length);
    for(int i = 0; i < filterCode.size();i++)
    {
        conjCode.at(i) = filterCode.at(filterCode.size() - i - 1);
    }
    conjCode = conj(conjCode);
    fftCode = fft_w->fft_cx(conjCode);
}

【讨论】：