2
$\begingroup$

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.

rect wave and its convolution

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()
$\endgroup$
  • 1
    $\begingroup$ The np.fft computations are correct; what is incorrect is that you expect these computations to give different results. It appears that you trying to verify Fourier transform properties of continuous-time signals by discretizing the latter and applying discrete Fourier transform (FFT). I would not recommend this approach due to subtle but critical differences between the continuous and discrete time domains. $\endgroup$ – Stelios Jan 17 '17 at 12:25
  • $\begingroup$ @Stelios You are exact right, I am trying to verify the convolution property of the continuous Fourier transform. Could you please elaborate on the "subtle but critical differences"? I was expecting a straightforward relation between the two, no more than some normalization. $\endgroup$ – Taozi Jan 17 '17 at 14:14
  • $\begingroup$ This is a rather deep topic. Just to give you an idea, consider for example the rectangular pulse signal. Whereas the continuous-time version has a well defined notion of "area", which is obtained by integration of the singal and is equal to 1, no such notion exists for the discrete-time version. The only "similar" quantity in the discrete-time version is the sum of the sample values, which is equal to 50. However this value will be different if you choose a different sampling period (dt). $\endgroup$ – Stelios Jan 17 '17 at 15:25
  • $\begingroup$ @Stelios Thank you for your reply. But how about the sampling "interval" of the rectangle function. If you discretize it into 50 points, then the gap between the samples is 1/50. With numerical integration like the rectangle rule, you still have area 1 -- and this is independent of how many sample you use if you consider it from a numerical integration perspective. $\endgroup$ – Taozi Jan 17 '17 at 15:38
  • $\begingroup$ Sure, however, this observation is relevant in the case where you want to compute the Fourier integrals numerically. In that case you should do the integration manually (or using a numerical integration function) and not use (I)FFT. The (I)FFT algorithm does not incorporate any sampling period information, i.e., it always implicitly assumes a normalized sampling period of 1. $\endgroup$ – Stelios Jan 17 '17 at 15:47
2
$\begingroup$

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.

$\endgroup$

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.