12

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)
jabaldonedo
  • 25,822
  • 8
  • 77
  • 77
user1163676
  • 155
  • 1
  • 1
  • 7

1 Answers1

19

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

jabaldonedo
  • 25,822
  • 8
  • 77
  • 77
  • 1
    Given \omega is typically used to denote angular frequency, you might want to drop the 2pi from your exponential, or change \omega to f. Also, you're missing a `j` in the initial definition. Otherwise, great answer! – Henry Gomersall Nov 14 '13 at 11:37
  • You are right Henry, I don't even know how I passed my signal processing course at uni with those mistakes! I will correct them. Thanks – jabaldonedo Nov 14 '13 at 11:39
  • 2
    Thank you very much. It should have been quite obvious indeed. By the way, do you know why the factor is not taken into account in standard fft routines ? My guess would be that most people use fft, then do something, then use ifft and so the factor is dropped to spare computational time. Right ? – user1163676 Nov 14 '13 at 13:14
  • 1
    Various "standard" FFT/IFFT libraries vary in the distribution or placement of this factor between the FFT and IFFT.. In your example, the factor placement preserves Parseval's energy equality. – hotpaw2 Nov 14 '13 at 17:48