# Convolution error when using DFT for non-periodic functions

I need help understanding why FFT-based computation of convolution on a finite domain, between two non-periodic functions often gets it wrong. Specific questions:

(1) why does FFT-based computation have error when interval [a, b] is short in the specific example below but no error when interval [a, b] is wider ;

(2) how can I anticipate this error when applying FFT-based approach to convolve functions?;

(3) Can this error be mitigated/reduced, e.g. different zero-padding options? Increasing discretisation size does not help at all.

Let's convolve two Gaussians on $$x\in [a, b]$$ for two choices of $$[a, b]$$: (1) $$[a,b]=[0,5]$$ and (2) $$[a,b]=[-1, 5]$$. We will compute convolution via DFT and compare against an exact result and a real-space evaluation.

Specific functions:

$$f_1(x)=\exp{(-x^2)}$$

$$f_2(x)=\exp{(-x^2)}$$

I want to compute the convolution on $$[a, b]$$ $$f(x)=\int\limits_a^b f_1(x) f_2(x-y) dx$$

Implementation:

import numpy as np

import matplotlib.pyplot as plt

import scipy.fft as fft
from scipy.special import erf

N = 64
a = 0
b = 5
x = np.linspace(a, b, N)
dx = np.diff(x)[0]

f1 = np.exp(-x**2)
f2 = np.exp(-x**2)

# exact convolution
fe = .5*np.exp(-x**2/2)*\
np.sqrt(np.pi/2)*\
(-erf((2*a-x)/np.sqrt(2))+erf((2*b-x)/np.sqrt(2)))

# computation on real domain
c = np.zeros(x.size)
for i in range(x.size):
c[i]=np.trapz(np.exp(-x**2)*np.exp(-(x-x[i])**2), x)

f3 = fft.ifft(fft.fft(f1)*fft.fft(f2)).real*dx

f31 = fft.ifft(fft.fft(np.concatenate((f1, np.zeros(N-1))))*\
fft.fft(np.concatenate((f2, np.zeros(N-1))))
).real*dx

x1 = np.arange(x[0], x[0] + f31.size*dx, dx)

plt.figure()

plt.subplot(2,1,1)
plt.plot(x, f1, label = 'f1(x)');
plt.plot(x, f2, '.', label = 'f2(x)');
plt.plot(x, fe, label = 'Exact f1*f2 (x)');
plt.plot(x, c, '.', label = 'f1*f2 (x) computed directly');
plt.legend()
plt.title(f"Exact convolution on [a,b]=[{a}, {b}]");

plt.subplot(2,1,2)
plt.plot(x, fe,'.', label = 'f1*f2: exact result');
plt.plot(x1+a, f31, label = 'f1*f2: computed with FFT')
plt.axvline(a)
plt.axvline(b)
plt.legend()
plt.title('FFT computation of convolution');
plt.tight_layout();



Results for $$[a, b]=[0, 5]$$: NOTICE THAT FFT APPROACH HAS HUGE ERROR

Results for $$[a, b]=[-2, 5]$$: NOTICE THAT FFT APPROACH HAS ALMOST NO ERROR

• Your math uses the continuous Fourier Transform but your code uses the Discrete Fourier Transform (DFT). These are fairly different animals: specifically multiplying DFTs is equivalent to CIRCULAR (not linear) convolution in the time domain. You also need to pick a sampling frequency that doesn't result in unacceptable amount of aliasing. Finite time domain signals have unlimited bandwidth, so you will always end up with some amount of aliasing Jun 27 at 18:49
• First of all, as you can see I applied zero-padding. So convolution I am taking is linear indeed. Secondly, changing sampling frequency will not affect the outcome. I can compute with N=64 or N=1024 discretisation points, and the problem persists. Jun 28 at 11:45

This seems to be because the significant non-zero component of the Gaussian before the region of interest $$[0,5]$$ is neglected from the fft calculation but not in the direct calculations.

If I do the calculation over $$[-10,10]$$ but just look at the range of interest, it works fine.

$$f(x) = e^{-x^2}$$

in the DFT calculation, this windowed version

$$f(x) = e^{-x^2}(u[x-a] - u[x-b])$$

is being used (where $$u[n]$$ is the unit step).

## Updated code below

import numpy as np

import matplotlib.pyplot as plt

import scipy.fft as fft
from scipy.special import erf

N = 64
a = 0
b = 5
aa = -10
bb = 10
x = np.linspace(aa, bb, 2*N)
dx = np.diff(x)[0]

f1 = np.exp(-x**2)
f2 = np.exp(-x**2)

# exact convolution
fe = .5*np.exp(-x**2/2)*\
np.sqrt(np.pi/2)*\
(-erf((2*aa-x)/np.sqrt(2))+erf((2*bb-x)/np.sqrt(2)))

# computation on real domain
c = np.zeros(x.size)
for i in range(x.size):
c[i]=np.trapz(np.exp(-x**2)*np.exp(-(x-x[i])**2), x)

f3 = fft.ifft(fft.fft(f1)*fft.fft(f2)).real*dx

#f31 = fft.ifft(fft.fft(np.concatenate((f1, np.zeros(N-1))))*\
#               fft.fft(np.concatenate((f2, np.zeros(N-1))))
#               ).real*dx

f31 = fft.ifft(np.multiply(fft.fft(np.concatenate((f1, np.zeros(2*N)))),
fft.fft(np.concatenate((f2, np.zeros(2*N)))))
).real*dx

x1 = np.arange(x[0], x[0] + f31.size*dx, dx)

plt.figure()

plt.subplot(3,1,1)
plt.plot(x, f1, label = 'f1(x)');
plt.plot(x, f2, '.', label = 'f2(x)');
plt.plot(x, fe, label = 'Exact f1*f2 (x)');
plt.plot(x, c, '.', label = 'f1*f2 (x) computed directly');
plt.legend()
plt.title(f"Exact convolution on [a,b]=[{a}, {b}]");

plt.subplot(3,1,2)
plt.plot(x, fe,'.', label = 'f1*f2: exact result');
plt.plot(x1+aa, f31, label = 'f1*f2: computed with FFT')
plt.axvline(a)
plt.axvline(b)
plt.legend()
plt.title('FFT computation of convolution');
plt.tight_layout();

plt.subplot(3,1,3)
plt.plot(x, fe,'.', label = 'f1*f2: exact result');
plt.plot(x1+aa, f31, label = 'f1*f2: computed with FFT')
plt.axvline(a)
plt.axvline(b)
plt.legend()
plt.title('FFT computation of convolution');
plt.xlim([a,b])
plt.tight_layout();

• Sorry, I don't know what you mean by "neglected". Let's forget that the data comes from a Gaussian. We have a well-defined function on [a, b]. We can convolve it with itself. Using FFT to compute such convolution leads to error. Why? Jun 28 at 11:44
• @pyatsysh How does that look now?
– Peter K.
Jun 29 at 18:22