Wrong amplitude of convolution using numpy fft

Question

I try to convolve a rectangle function in [-1/2, 1/2] with itself using fft. The convolution should be a tent shaped function, see figure below.

The code is below. In the 3rd to last line I add /50 so it appears correct - I have no idea what the normalization factor should be. And my question is how to normalize the power of a fft result, so the subsequent application of ifft gives the correct result?

import numpy as np
from numpy.fft import fft, ifft, fftshift, ifftshift
N = 1000
t=np.linspace(-10, 10, N)
dt = t[1]-t[0]

data = 1.0*(np.abs(t)<0.5)
data_fft = np.fft.fft(data)

plt.hold(1)
plt.plot(t, (ifft(data_fft**1)), 'r')
plt.plot(t, ifftshift(ifft(data_fft**2)) / 50.0, 'g')
plt.xlim(-4,4)
plt.show()

Answer 1

You actually do recover the convolution, but as it is discussed in the comments, there is a normalization issue due to discretization.

According to the documentation, fft is implemented like this:

$$ A_k = \sum_{m=0}^{n-1} a_m \exp \{ - 2\pi i \frac{mk}{n} \} $$

with $A_k$ being the Fourier-coefficients, $a_m$ the $m$-th element of your signal vector and $n$ the length of the signal.

Squaring this gives you

$$ A_k^2 = \sum_{m=0}^{n-1} \sum_{m'=0}^{n-1} a_m a_{m'} \exp\{ - 2\pi i \frac{(m+m')k}{n} \} \} $$

Now, applying ifft to the squared Fourier-transform gives you, using the ifft-definition from the documentation:

$$ \text{ifft}(A_k^2)_{m''} = \frac{1}{n} \sum_{k=0}^{n-1} \sum_{m=0}^{n-1} \sum_{m'=0}^{n-1}a_m \exp\{ - 2\pi i \frac{(m+m'-m'')k}{n} \} \}$$

With the observation, that

$$ \frac{1}{n} \sum_{k=0}^{n-1} \exp\{ - 2\pi i \frac{(m+m'-m'')k}{n} \} = \delta_{m+m', m''} $$

you end up with

$$ \text{ifft}(A_k^2)_{m''} = \sum_{m=0}^{n-1} \sum_{m'=0}^{n-1} a_m a_{m'} \delta_{m+m', m''} = \sum_{m=0}^{n-1} a_m a_{m'' - m} $$

This is actually how np.convolve is defined (except for some padding). If you use np.convolve on your data, you end up with the same result (except for some padding), so within the numpy-world, you did exactly what you set out to do, i.e. verify, the convolution property of the Fourier transform. As noted in the comments however, neither fft nor convolve "know" anything about your descretization, so you have to take care of that manually by multiplying the results with dt.