1D Fast Convolution without FFT

Question

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?

Solution 1

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.

Solution 2

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

Comments

[templatetypedef]

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.