Problem implementing convolutions exactly with the FFT

Question

I'm trying to perform convolutions as defined mathematically $f \star g (\tau)= \int_{\mathcal{R}}f(t-\tau)g(t) dt$ in a numerical simulation. Hence, my signal is a sampling of points $f(x_i)$.

I need the computation to be precise, unbiased. Indeed, the simulation is there to show very refined effects, so I cannot discard small errors of principle. Note also that this convolution will be performed very many times over large arrays, hence it needs be efficient.

My question is if and how can I correctly compute this convolution? I have tried with the FFT algorithms already implemented in Julia (or Python, same story), but I seem to have systematic errors, as shown below. Are these the product of my lack of understanding? How else should I do?

This code is a 'test' to check that convoluting two boxes numerically gives the analytical result:

using Plots
using FFTW
using DSP


function box(x)
    if abs(x) < 0.5
        return 1
    else
        return 0
    end
end

function analytical_convolved_box(x)
    if x <= 0 && x > -1
        return x + 1
    elseif x > 0  && x < 1
        return 1 - x
    else
        return 0
    end
end

x = range(-2, stop = 2, length = 20)
dx = x[2] - x[1]

input = box.(x)

l_half = convert(Int,length(input)/2)

output = dx *conv(input, input)[l_half+1:3l_half]

plot(x, box.(x), color = :blue)
plot!(x, box.(x), seriestype = :scatter, color = :blue)
plot!(x, analytical_convolved_box.(x), color = :green)
plot!(x, analytical_convolved_box.(x), seriestype = :scatter , color = :green)
plot!(x, output, color = :black)
plot!(x, output, seriestype = :scatter, color = :black)

and for 11 points:

x = range(-2, stop = 2, length = 11)
dx = x[2] - x[1]

input = box.(x)

l_half = convert(Int,length(input)/2 - 0.5)

output = dx *conv(input, input)[l_half+1:3l_half+1]

I'd like to understand things clearly. Note that in both examples it wasn't really clear which indexes to use, so I had to essentially try and guess.

Answer 1

It looks to me that the FFT convolution algorithm is doing what is expected here. Remember that you work with discretized signals, and therefore, discretized convolution.

There are a few ways to implement the discrete convolution numerically. The slight differences basically come from the way to deal with the finiteness of each signal.

The discrete convolution for two infinite arrays $f[m]$ and $g[m]$ (with $m \in \mathbb{Z}$) can be written as:

$$\left( f *g\right)[n] = \sum_{m\in \mathbb{Z}} f[m] g[n-m] $$

For discrete arrays of size $N$, the simplest method to use (which takes full advantage of the FFT), is the circular convolution:

$$\left( f *g\right)[n] = \sum_{m = 0}^{N-1} f[m] g[(n-m) \, \mathrm{mod}(N)], $$

where the indexes taken modulo $N$ (so $g(-1) \equiv g(N-1)$ and so on).

This is just the way discrete convolution work. You can treat the boundaries of each array differently (for instance, extending each array with zeroes, or repeating the last value of the array indefinitely), but for the most part, it won't change the result (especially in the center).

From this, what is the expected result of your calculation? If you take an array $f$ of length $11$, with $f[5] = f[6] = f[7] = 1$ and $f[m] = 0$ otherwise (the array from your second example), then you can convince yourself that $(f * f)[0] = 3$ (because the three ones in the middle "overlap" nicely), $(f*f)[1] = 2$ (only two out of the three ones "overlap"), $(f*f)[2] = 1$, $(f*f)[3] = 0$, ...

Up to shifting the array you obtain to center it, and multiplying by your step size $dx = 0.4$, this is exactly what the FFT algorithm gives you.

If you want to work with the continuous version of the convolution, you should probably just take more points. Right now you are working with a gate function so it is easy to find the exact result theoretically, but imagine that you work with anything else (trigonometric functions, power functions, gaussian, etc...). You are still limited by your sampling rate (the number of points per unit distance). How is the program supposed to interpolate for the missing values between two points?