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
        return 0

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

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)

enter image description here

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]

enter image description here

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.

  • 4
    $\begingroup$ Have you tried taking signals with more points? You are essentially trying to estimate the convolution of a discontinuous function using an FFT with only 10 or 11 modes. Using the FFT to compute convolutions on samples to estimate the convolution of the true signal only makes sense if the signal can be accurately reconsctructed using that number of Fourier modes. In this case, the true signal is discontinuous so it takes more Fourier modes to reconstruct and so Fourier-based methods may take more modes to provide accurate results $\endgroup$ – whpowell96 Apr 26 '20 at 23:32
  • $\begingroup$ @whpowell96 yes, it works better with more points... However, I feel I should I obtain the exact analytical equivalent up to floating point error... Something doesn't feel right $\endgroup$ – Comrad dau Apr 27 '20 at 14:00
  • $\begingroup$ You should not obtain the analytic solution up to floating point error if you are sampling a signal at only 10/11 points because your signal has significant reconstruction error when represented in only the first 10/11 Fourier mdoes. That is a property of your continuous signal and there is nothing you can do about it $\endgroup$ – whpowell96 Apr 27 '20 at 18:18

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?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.