I try to validate my understanding of Numpy's FFT with an example: the Fourier transform of exp(-pi*t^2)
should be exp(-pi*f^2)
when no scaling is applied on the direct transform.
However, I find that to obtain this result I need to multiply the result of FFT by a factor dt
, which is the time interval between two sample points on my function. I don't understand why. Can anybody help ?
Here is a sample code:
# create data
N = 4097
T = 100.0
t = linspace(-T/2,T/2,N)
f = exp(-pi*t**2)
# perform FT and multiply by dt
dt = t[1]-t[0]
ft = fft(f) * dt
freq = fftfreq( N, dt )
freq = freq[:N/2+1]
# plot results
plot(freq,abs(ft[:N/2+1]),'o')
plot(freq,exp(-pi*freq**2),'r')
legend(('numpy fft * dt', 'exact solution'),loc='upper right')
xlabel('f')
ylabel('amplitude')
xlim(0,1.4)
Be careful, you are not computing the continuous time Fourier transform, computers work with discrete data, so does Numpy, if you take a look to numpy.fft.fft
documentation it says:
numpy.fft.fft(a, n=None, axis=-1)[source]
Compute the one-dimensional discrete Fourier Transform.
This function computes the one-dimensional n-point discrete Fourier Transform (DFT) with the efficient Fast Fourier Transform (FFT) algorithm
That means that your are computing the DFT which is defined by equation:
the continuous time Fourier transform is defined by:
And if you do the maths to look for the relationship between them:
As you can see there is a constant factor 1/N
which is exactly your scale value dt
(x[n] - x[n-1]
where n is in [0,T] interval is equivalent to 1/N
).
Just a comment on your code, it is not a good practice to import everything from numpy import *
instead use:
import numpy as np
import matplotlib.pyplot as plt
# create data
N = 4097
T = 100.0
t = np.linspace(-T/2,T/2,N)
f = np.exp(-np.pi*t**2)
# perform FT and multiply by dt
dt = t[1]-t[0]
ft = np.fft.fft(f) * dt
freq = np.fft.fftfreq(N, dt)
freq = freq[:N/2+1]
# plot results
plt.plot(freq, np.abs(ft[:N/2+1]),'o')
plt.plot(freq, np.exp(-np.pi * freq**2),'r')
plt.legend(('numpy fft * dt', 'exact solution'), loc='upper right')
plt.xlabel('f')
plt.ylabel('amplitude')
plt.xlim([0, 1.4])
plt.show()