Why computing the DFT from the Ricker Pulse isn't the same as the continuous FT function?

94 Views Asked by At

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

time_domain

frequency domain

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

comparison

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.

comparison 2