FFT Convolution and Zero-Padding

Continuous and Discrete Convolution

Given two continuous functions $f$ and $g$, their convolution is usually defined as:

Quite naturally, for discrete functions $f,g$ defined on the set $\mathbb{Z}$ of integers, the discrete convolution is defined as the sum

Cyclic Convolution with FFTs

The cyclic convolution of two periodic sequences $f_n, g_n$ of period N, denoted by $(f_n *_N g_n)$, will also have period N, and for convenience, it is defined for $0 \le n \le N-1$. 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 $O(n \log n)$ operations, insted of $O(n^2)$.

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 $f*g$, if $f$ is supported on $[0,N-1]$ and $g$ is supported on $[0,M-1]$.

To do this, let’s note that the support of $g[n-\cdot]$ is $[n-(M-1),n]$. Then, for the product $f[m] g[n-m]$ to be nonzero for some value of $m$, it must be that $[n-(M-1),n] \cap [0,N-1] \ne \emptyset$. Equating opposite extremes of both intervals and solving for $n$, we obtain $n \in [0, N + M - 2]$.

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 makes g periodic.
• Extende f and g 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.