# Strategy / Method for Implementation of the Fastest 1D Linear Convolution / Correlation

I am after the method and its implementation of the fastest convolution of 1d signals:

$$y \left[ n \right] = \left( x \ast h \right) \left[ n \right]$$

Where $$\ast$$ is the linear convolution operator.

Let's assume that we have infinite storage or memory but we have finite time. What strategy or algorithm can we use for a faster discrete time convolution operation?

Could someone compare the known methods for implementing convolution: Direct Convolution, Frequency Domain Convolution / FFT Convolution and Overlap and Save / Overlap and Add methods? Is there something even faster?

• If I had an algorithm for faster convolution, I'd publish it in a journal (before posting it here). This is a problem that has been worked on for decades; I'd wager that you won't get many answers. Let us know if you make progress, though :) – MBaz Oct 21 '18 at 23:14
• Bless you mate :D. I sure will. – Shkodrani Oct 21 '18 at 23:30
• Please consider also in your comparison of optimized approaches the well understood relationship using fft's: $conv = ifft(fft(a)*fft(b))$. Also note that cross correlation is also solved with a similar fft relationship: $xcorr = ifft(fft(a)*conj(fft(b))$, where "conj" is the complex conjugate. – Dan Boschen Oct 22 '18 at 0:23
• @DanBoschen Indeed, AFAIK for non-short sequences that approach is the fastest known. Since we don't know that current FFT algorithms are the fastest possible, the implication is that we don't know if there are faster algorithms for convolution either :) – MBaz Oct 22 '18 at 0:31
• @MBaz Yes indeed. That's why I am rooting for Shkodrani and look forward to his FCATW implementation! – Dan Boschen Oct 22 '18 at 0:42

I compared 3 implementations for Linear Convolution of 1D signals:

1. Direct - Using MATLAB's conv() funciton.
2. Overlap and Save - Implemented in MATLAB with tuned loop to prevent allocation and optimal choice of the DFT window.
3. Frequency Domain - Using MATLAB' fft() and proper padding to implement Linear Convolution using Circular Convolution.

For various lengths of the input signal and the kernel I got this result: I only compared cases the signal is not shorter than the kernel hence the upper triangle (Dark Blue) is invalid. For the lower triangle I coded the fastest method as following:

1. Bright Blue - Direct.
2. Green - Overlap and Save.
3. Yellow - Frequency Method.

So, there is nor practical reason to use Overlap and Save.
For kernels with length up to ~400 samples it is better to use Direct Method.
For longer kernels it is better to use Frequency Domain.

Now, the actual number of samples as the border between Frequency Domain and Direct will be different for different CPU's and Memory Bandwidth (One can run the script on his computer to see).
Yet as a guideline I'd say:

1. If you can hold all the data in Memory - Don't use Overlap and Save.
2. Unless your kernel is longer few hundreds of samples use direct.
3. For kernels with more than few hundreds samples use Frequency Domain.

The full MATLAB code is available on my StackExchange Signal Processing Q52760 GitHub Repository.

• The MATLAB code allocates a lot of memory so it might be faster even as MEX with direct call to FFTW. – David Sep 11 at 6:02
• You're right. I'd do it in case I need it in my code. – Royi Sep 11 at 6:06

It's not very clear whether you ask for :

• a more efficient definition of the standard convolution operator between two sequences $$x[n]$$ and $$h[n]$$ given as:

$$y[n] = x[n] \star h[n] = \sum_{k=-\infty}^{\infty} h[k] x[n-k]$$ which would yield the same $$y[n]$$ per given $$n$$, but with a different formulation which would require less number of arithmetic operations and hence will be more efficient to compute.

or:

• a more efficient implementation of the standard convolution operator, (instead of using the direct sum approach) which would produce the same $$y[n]$$ per given $$n$$, based on the same formula above but using an architecture that lends the same result with less number of operations.

Note that whether these two approaches are indeed distinct or just the same thing stated differently is also a matter of debate.

Neverthless, for the second approach, as Dan Boschen has already stated, the Fast Fourier Transform (FFT) algorithm, which is indeed a fast compuation of DFT (discrete Fourier transform) can be used to implement the standard convolution operator at a much faster rate (due to less number of MACs) when sequences are long enough.

Furthermore, again for the second approach, based on your statement that you have infinite memory (and I would add that also an infinite (very large) number of FPU ALUs) then you can implement the direct sum approach using a massively parallel architecture. You would still be performing the same number of MACs but since they are in parallel, the execution steps will be much shorter, and therefore faster.

For example, assume those two sequences above are defined in $$0 \leq n \leq N_x$$ and $$0\leq n\leq N_h$$ respectively. The standard convolution will produce $$y[n]$$ which is defined for $$0\leq n \leq N_y = N_x + N_h -1$$.

Then a direct sum approach performed on a single core (single ALU) serial execution architecture would compute about $$\min\{N_h,N_x\}$$ MACs per each output sample of $$y[n]$$, totalling a $$\min\{N_h,N_x\}\times N_y$$ MACs. (The exact number should take the partially overlapping edge cases in detail and is less than this upper bound.) Due to serial processing, each output time is added onto.

Now if all those $$N_y$$ output samples are computed in parallel, then $$N_y \times$$ times faster execution is possible (still doing the same number of MACs though).

Furthermore, computation of each output sample would normally take $$\min\{N_h,N_x\}$$ MACs and in a serial computation that would also require $$\min\{N_h,N_x\}$$ steps (apprx CPU cycles). But a parallel computation of that dot product would also yield substantial execution (step) time reduction.

As a result, the parallel execution time of a direct sum approach of the standard convolution operator reduces into that of a single dot product time between two vectors of length $$\min\{N_h,N_x\}$$, which can be done (again in parallely) in $$\log_2 \min\{N_h,N_x\}$$ MAC steps. Quite fast compared to serial $$N_y \times \min\{N_h,N_x\}$$ steps.

Note that we have ignored any kind of overhead; data movement, memory management, synchronization, etc. These can very easily overwhelm any gains you would expect from a parallel execution.

• I did a full comparison of Run Time in my answer. Sometimes asymptotic properties are not a good guideline for real life. – Royi Apr 27 at 11:18
• @Royi Run time for the best possible implementations? Or Matlab? Counting adds and mult combined with profiling the best code you can get on the hardware that you are going to use is often a strong combination. – Knut Inge May 27 at 13:06