Continuous and Discrete Convolution
Given two continuous functions and , their convolution is usually defined as:
Quite naturally, for discrete functions defined on the set of integers, the discrete convolution is defined as the sum
Cyclic Convolution with FFTs
The cyclic convolution of two periodic sequences of period N, denoted by , will also have period N, and for convenience, it is defined for . It is given by:
A direct implementation of this formula in Julia (taking into account that Julia counts from 1, and the periodicity of g
) is the following:
function DirectCyclicConv1D(f,g)
N = length(f)
Conv = zeros(N)
for n ∈ 1:N, m ∈ 1:N
if n-m+1 > 0
Conv[n] = Conv[n] + f[m] * g[n-m+1]
else
# make g periodic
Conv[n] = Conv[n] + f[m] * g[N+(n-m+1)]
end
end
return Conv
end
In view of the Convolution Theorem (see Wiki), we can employ FFTs to compute this transformation in operations, insted of .
FastCyclicConv1D(f,g) = ifft(fft(f).*fft(g))
DirectCyclicConv1D
should be equal (up to floating point precision) to FastCyclicConv1D
, which is something we can test as follows:
using Test
f = [1;3;2;5;2;3;2]
g = [3;6;4;5;3;4;2]
@test DirectCyclicConv1D(f,g) ≈ FastConv1D(f,g)
From Cyclic to Linear Convolution
The cyclic convolution can be useful in certain contexts, but more often than not we are interested in the linear (or ordinary) convolution of two sequences of finite length.
Importantly, if we want to employ FFTs to compute this convolution, we need to employ zero-padding.
How can we do this?
Well, let’s make sure that we know what we want to compute in the first place, by writing a direct convolution which will serve us as a test function for our FFT code.
As a first step, let’s consider which is the support of , if is supported on and is supported on .
To do this, let’s note that the support of is . Then, for the product to be nonzero for some value of , it must be that . Equating opposite extremes of both intervals and solving for , we obtain .
We therefore want to compute:
The most straightforward way to implement this is by modifying the cyclic convolution code written above as little as possible.
Actually, a convenient way to think about this problem is that we will implement a cyclic convolution (to use FFTs) in addition to zero padding, so the result coincides with the linear convolution we want to compute.
Linear Convolution Length and Zero-Padding
While the cyclic convolution of two length N
signals has length N
, the linear convolution of two signals of lenght N
and M
has length N+M-1
. So at minimum, it seems that we need to pad both signals with enough zeros so they both have length N+M-1
. That way, their cyclic convolution will have (at least) the correct size of N+M-1
.
As we are performing summation, extending the data with zeros clearly won’t affect the result. So we could do all the zero-padding we want to do.
As a second step, we might want to make our convolution function to resemble a cyclic convolution. This will get us closer to our goal of FFT acceleration of this function.
Let’s first write an easy implementation against which to test the FFT implementation we’ll do later. In order to do this, we’ll just do the following:
- Change all lengths to
N+M-1
- Remove the
else
statement that makesg
periodic. - Extende
f
andg
with zeros, which doesn’t alter the result, and allows for succinct code.
function DirectLinearConvolution(f,g)
N = length(f)
M = length(g)
g_pad = [ g; zeros(N-1) ] # length N+M-1
f_pad = [ f; zeros(M-1) ] # length N+M-1
Conv = zeros(N+M-1)
for n=1:N+M-1
for m=1:N+M-1
if n-m+1 > 0
Conv[n] = Conv[n] + f_pad[m] * g_pad[n-m+1]
end
# n+1 <= m
end
end
return Conv
end
The purpose of zero-padding g
is simply to be able to evaluate g_pad[n-m+1]
for the required values of n-m+1
while avoiding the programming error of trying to use an undefined vector. Zero-padding of f
allows us to use a larger inner loop for m=1:N+M-1
, which we need to traverse the full support of g_pad[n-m+1]
.
As I always like to test my code (even this one, which is intended for testing our FFT code!), in this case I’ve employed the conv
function in the DSP
package:
using Test, DSP
f = [1;3;2;5;2;3;2]
g = [3;6;4;5;3;4;2]
@test DSP.conv(f,g) ≈ DirectLinearConvolution(f,g)
Two take-aways:
- Our code for the linear convolution looks very similar to the cyclic convolution. Once signals are zero-padded, the only difference lies in the
else
statement of the cyclic convolution. - Zero-padding was added to compute the linear convolution with succinct code, and there were no FFTs or frequencies involved (yet).
New let’s revisit the else
statement of the cyclic convolution. It is triggered when n-m+1 < 0
, and in that case the computed quantity is f[m] * g[N+(n-m+1)] = 0
. It is easy to check that terms in this form are actually zero…
So the loop we just wrote actually coincides with a cyclic convolution!
As we already have a fast implementation for the cyclic convolution, we can now build our fast linear convolution as follows:
function FastLinearConvolution(f,g)
N = length(f)
M = length(g)
f_pad = [ f; zeros(M-1) ]
g_pad = [ g; zeros(N-1) ]
return FastCyclicConv1D( f_pad, g_pad )
end
And of course, test it with:
@test FastLinearConvolution(f,g) ≈ DirectLinearConvolution(f,g)
Conclusion
That was it for the basics of FFT convolution in one dimension. Many more things are possible from here, like applying filters, image processing, computing cross- and auto-correlation, etc.