2011-08-30 95 views
5

我需要对2个大数组进行一维卷积。我在C#中使用这段代码,但它需要很长时间才能运行。1D无FFT快速卷积

我知道,我知道! 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?

+5

虽然你说过不要问,为什么你不能使用FFT?如果这是针对明确禁止的类项目,则应该将其标记为家庭作业。 – templatetypedef

+0

C#可以调用CUDA吗?如果可以的话,你可以使用并行指令,这相当大地加速了幼稚卷积。或者你可以使用Winograd变换或其他东西(不是Cooley-Tukey经典FFT,如果这足够满足你的“不FFT”规则)。或者,如果您知道输入或滤波器的某些信息(例如只有某些频率存在或某种情况),则可以使用该知识。你将不得不更加具体地了解你的限制以及你可能拥有的任何外部知识。 – mtrw

回答

5

卷积是数字相同的与额外的缠绕步骤的多项式乘法。因此,可以使用所有的多项式和大整数乘法算法来执行卷积。

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

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

0

这里有两种可能会导致较小的加速,但你需要测试以确保。

  1. 展开内部循环以删除一些测试。如果知道过滤器长度将总是一个倍数(如果有的话),这将更容易。
  2. 颠倒循环的顺序。 filter.length是否通过整个阵列。这在内部循环中的解引用较少,但可能具有较差的缓存行为。
4

你可以减少索引访问次数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 < filterLengthinputLengthresultLength ≤)

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; 
} 
0

您可以使用特殊的IIR滤波器。 然后处理,如:

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

我认为它更快。