Question

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)
Was it helpful?

Solution

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:

enter image description here

the continuous time Fourier transform is defined by:

enter image description here

And if you do the maths to look for the relationship between them:

enter image description here

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()

enter image description here

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top