1D Fast Convolution without FFT

11
votes

I need an 1D Convolution against 2 big arrays. I'm using this code in C# but it takes a loooong time to run.

I know, i know! FFT convolutions is very fast. But in this project i CAN'T use it. It is a constraint of the project to not use FFT (please don't ask why :/).

This is my code in C# (ported from matlab, by the way): 

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

So, anyone knows any fast convolution algorithm widthout FFT?

4
c#optimizationconvolution
Although you've said not to ask, why can't you use the FFT? If this is for a class project where it's explicitly prohibited, you should probably tag this as homework. - templatetypedef
Can C# call CUDA? You could use parallel instructions if so, which speeds up naive convolutions considerably. Or you could use the Winograd transform or something (not the Cooley-Tukey classic FFT, if that's far enough away to satisfy your "no FFT" rule). Or if you know something about the input or filter (like only certain frequencies are present or something) you can use that knowledge. You're going to have to be more specific about your constraints and any outside knowledge you might have. - mtrw

4 Answers

6
votes

Convolution is numerically the same as a polynomial multiplication with an extra wrap-around step. Therefore, all the polynomial and large integer multiplication algorithms can be used to perform convolution.

FFT is the only way to get the fast O(n log(n)) run-time. But you can still get sub-quadratic run-time using the divide-and-conquer approaches like Karatsuba's algorithm.

Karatsuba's algorithm is fairly easy to implement once you understand how it works. It runs in O(n^1.585), and will probably be faster than trying to super-optimize the classic O(n^2) approach.

6
votes

You could reduce the number of indexed accesses to result, as well as the Length properties:

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

If you further split the outer loop, you can get rid of some repeating conditionals. (This assumes 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;
}
0
votes

You can use a special IIR filter. Then process that as like:

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

I think it's faster.

0
votes

Here are two possibilities that may give minor speedups, but you'd need to test to be sure.

  1. Unroll the inner loop to remove some tests. This will be easier if you know the filter length will always be a multiple if some N.
  2. Reverse the order of the loops. Do filter.length passes over the whole array. This does less dereferencing in the inner loop but may have worse caching behavior.