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)
The frequency axis is identical to that of the two-sided power spectrum. The amplitude of the FFT is related to the number of points in the time-domain signal.
Let X = fft(x) . Both x and X have length N . Suppose X has two peaks at n0 and N-n0 . Then the sinusoid frequency is f0 = fs*n0/N Hertz.
This function computes the one-dimensional n-point discrete Fourier Transform (DFT) with the efficient Fast Fourier Transform (FFT) algorithm [CT]. Parameters aarray_like. Input array, can be complex. nint, optional. Length of the transformed axis of the output.
The fft function which uses the functionality of the SciPy package works in a way that, it uses the basic data structures that are used in the numpy arrays, in order to create a module that is required for scientific calculations and programming.
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()
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With