Question
I'm trying to use Python's scipy library to compute the DFT of the Ricker wavelet function and compare it with the analytical Fourier-transformed version of the same function. When I compare the results, they do not match. I want to understand why.
Perhaps related, but also unanswered: DFT is not a sampling of FT?
Details
I'm using Python for the experiment. The Ricker pulse function in time domain is defined as follows:
def ricker_pulse_time(w0, t):
return (1. - .5 * w0**2 * t**2) * np.e ** (-.25 * w0**2 * t**2)
And, in the frequency domain:
def ricker_pulse_freq(w0, w):
return (2. * w**2) / (np.sqrt(np.pi) * w0**3) * \
np.e ** (-w**2 / w0**2)
As an example, I'm computing the function with 1Hz frequency:
hz_to_rad = lambda x: 2. * np.pi * x
w0 = hz_to_rad(1.0) # [Rad/s]
time_step = 1e-6 # [s]
t = np.arange(-2.0, 2.0, time_step)
y = ricker_pulse_time(w0, t)
N = len(y)
omega = np.arange(0.0, 40.0, 1e-4)
yy = ricker_pulse_freq(w0, omega)
plt.plot(t, y);
plt.show();
plt.plot(omega, yy);
Finally, to compute the DFT, I'm using scipy's implementation as follows:
y_f = fftpack.fft(y)
sample_freq = fftpack.fftfreq(N, d=time_step)
And comparing with the analytical solution previously commented:
plt.plot(omega, yy, 'r--')
plt.plot(sample_freq[0:N//2], np.abs(y_f[0:N//2]), 'b-')
plt.xlabel('Frequency [Rad/s]')
plt.ylabel('y_f')
plt.xlim(0, 40);
plt.show();
Question: Why are the results different?
Note 1: I'm only getting half the points, ignoring the negative part of the solution.
Note 2: I also tryied scalonating the result using 2/N * abs(y_f), but results are also not correct.



