2
$\begingroup$

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 enter image description here

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

$\endgroup$
2
  • $\begingroup$ 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 $\endgroup$
    – Hilmar
    Commented Jun 27, 2023 at 18:49
  • $\begingroup$ 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. $\endgroup$
    – pyatsysh
    Commented Jun 28, 2023 at 11:45

1 Answer 1

3
$\begingroup$

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.

Re-running and zooming into the region of interest.

Instead of using

$$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).

Illustration of windowing of exp(-x^2).


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();
$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – pyatsysh
    Commented Jun 28, 2023 at 11:44
  • $\begingroup$ @pyatsysh How does that look now? $\endgroup$
    – Peter K.
    Commented Jun 29, 2023 at 18:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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