无 FFT 的一维快速卷积答案

【问题标题】：1D Fast Convolution without FFT无 FFT 的一维快速卷积
【发布时间】：2011-08-30 01:47:40
【问题描述】：

我需要针对 2 个大阵列的一维卷积。我在 C# 中使用此代码，但运行需要很长时间。

我知道，我知道！ FFT 卷积非常快。但在这个项目中我不能使用它。不使用 FFT 是项目的一个限制（请不要问为什么：/）。

这是我的 C# 代码（顺便说一下，从 matlab 移植）：

var result = new double[input.Length + filter.Length - 1];
for (var i = 0; i < input.Length; i++)
{
    for (var j = 0; j < filter.Length; j++)
    {
        result[i + j] += input[i] * filter[j];
    }
}

那么，有人知道任何快速卷积算法widthout FFT吗？

【问题讨论】：

虽然你说不要问，为什么不能用FFT？如果这是针对明确禁止的课堂项目，您可能应该将其标记为家庭作业。
C#可以调用CUDA吗？如果是这样，您可以使用并行指令，这会大大加快简单的卷积。或者您可以使用 Winograd 变换或其他东西（不是 Cooley-Tukey 经典 FFT，如果它足够远以满足您的“无 FFT”规则）。或者，如果您对输入或滤波器有所了解（例如仅存在某些频率或其他内容），您可以使用该知识。您必须更具体地了解您的限制条件以及您可能拥有的任何外部知识。

标签： c# optimization convolution

【解决方案1】：

卷积在数值上与带有额外回绕步骤的多项式乘法相同。因此，所有的多项式和大整数乘法算法都可以用来进行卷积。

FFT 是获得快速 O(n log(n)) 运行时间的唯一方法。但是您仍然可以使用 Karatsuba's algorithm 等分而治之的方法获得亚二次运行时间。

一旦您了解了 Karatsuba 的算法是如何工作的，它就很容易实现。它在 O(n^1.585) 中运行，并且可能比尝试超级优化经典的 O(n^2) 方法更快。

【讨论】：

【解决方案2】：

您可以减少对result 以及Length 属性的索引访问次数：

int inputLength = filter.Length;
int filterLength = filter.Length;
var result = new double[inputLength + filterLength - 1];
for (int i = resultLength; i >= 0; i--)
{
    double sum = 0;
    // max(i - input.Length + 1,0)
    int n1 = i < inputLength ? 0 : i - inputLength + 1;
    // min(i, filter.Length - 1)
    int n2 = i < filterLength ? i : filterLength - 1;
    for (int j = n1; j <= n2; j++)
    {
        sum += input[i - j] * filter[j];
    }
    result[i] = sum;
}

如果你进一步拆分外循环，你可以去掉一些重复的条件。（假设 0 filterLength &leq; inputLength &leq; resultLength）

int inputLength = filter.Length;
int filterLength = filter.Length;
int resultLength = inputLength + filterLength - 1;

var result = new double[resultLength];

for (int i = 0; i < filterLength; i++)
{
    double sum = 0;
    for (int j = i; j >= 0; j--)
    {
        sum += input[i - j] * filter[j];
    }
    result[i] = sum;
}
for (int i = filterLength; i < inputLength; i++)
{
    double sum = 0;
    for (int j = filterLength - 1; j >= 0; j--)
    {
        sum += input[i - j] * filter[j];
    }
    result[i] = sum;
}
for (int i = inputLength; i < resultLength; i++)
{
    double sum = 0;
    for (int j = i - inputLength + 1; j < filterLength; j++)
    {
        sum += input[i - j] * filter[j];
    }
    result[i] = sum;
}

【讨论】：

【解决方案3】：

您可以使用特殊的 IIR 过滤器。然后像这样处理：

y(n)= a1*y(n-1)+b1*y(n-2)...+a2*x(n-1)+b2*x(n-2)......

我认为它更快。

【讨论】：

【解决方案4】：

这里有两种可能会稍微加快速度，但您需要进行测试才能确定。

展开内部循环以删除一些测试。如果您知道过滤器长度将始终是 N 的倍数，这将更容易。
颠倒循环的顺序。做filter.length 传递整个数组。这会减少内部循环中的取消引用，但可能会有更糟糕的缓存行为。

【讨论】：