Reconstruction of a Ricker Wavelet using inverse discrete fourier transform - signal cut in a half?

Question

I am new here and new to DSP, so maybe my question is really basic.

I have the formula for the Ricker wavelet (Mexican Hat) in frequency-domain and I wish to do an inverse Fourier transform to recover my original signal in time-domain. I am using python numpy.fft module for this.

For some reason, instead of a Ricker wavelet (https://wiki.seg.org/wiki/Dictionary:Ricker_wavelet), I am obtaining a divided version of the signal, like it is aliased or cut in half or lagged (yes, I'm confused).

Do I have to change the order of my time vector accordingly to the frequency vector ? What is the reason for this ? Or is it something else that I am missing ?

My goal is to retrieve a Ricker wavelet centered in zero (or even lagged), but I don't know why my results are like these and how to justify flipping or slicing my time vector.

Please find below my code which also generate the plots. Please let me know if you need any further information.

Thanks in advance, Luis

import matplotlib.pyplot as plt
import numpy as np
# Dummy signal length
nsamples = 338
    
dt = 1.6199375667655787e-10

freq = np.fft.fftfreq(len(trace),d=dt)

# peak angular frequency
omega_p = 2*np.pi*250e6

#Using only the positive frequencies for the Ricker Wavelet calculation
omega = 2*np.pi*freq[0:169]

# Ricker Wavelet in Frequency Domain
S_desired = (2/np.sqrt(np.pi))*((omega**2)/(omega_p**3))*np.exp(-(omega**2)/(omega_p**2))

# Appending the Ricker Wavelet values
S_flip = np.flip(S_desired).copy()

S = np.concatenate((S_desired,S_flip))


S_desired_time = np.fft.ifft(S) 
    
time = time=np.arange(0,nsamples*dt,dt)

plt.plot(freq,np.abs(S),'r',label='Power spectrum Ricker wavelet')
plt.xlabel('Frequency [Hz]')
plt.figure()
plt.plot(time,S_desired_time,label='IFFT of the Ricker Wavelet')
plt.xlabel('time [s]')

Answer 1

For the centring on t=0: the frequency domain doesn't have any absolute time information, so it's up to you to define where your t=0 is (after all, when you do an fft of a time signal, you just supply the samples, and the dt is only used to set the frequency scale).

AS for the shape: you have to remember that the discrete FT of a signal that has a finite number of samples is the same one as with the signal repeated:

And so, you have a phase shift between the signal you get and the one you want of 1/2 your sample length.

The easiest solution to visualise your signal is just to shift the first and last part of your signal, by using for example:

S_desired_time = np.fft.ifftshift((np.fft.ifft(S))
plt.plot(time-time[-1]/2,S_desired_time.real) # to centre on 0

(note that the output of an ifft is complex, so if you want to just should plot the real part use signal.real -- plt.plot() should do it by default but it doesn't hurt to specify it)

Another way is to remember that phase shifting a time signal is the same a multiplying the Fourier transform with an exponential:

$$\mathscr{F}\big\{x(t-t_0)\big\}=X(f)e^{-j2\pi f t_0}$$

You can add this lines to you code:

timeshift=1/2
f = np.linspace(0, nsamples-1, nsamples)
phaseshift  = (np.exp((-2*np.pi*1j*f*timeshift+2*np.pi*1j*timeshift*nsamples)))
S=S*phaseshift
S_desired_time = (np.fft.ifft(S))

As abs(phaseshift)=1, the magnitude of your frequency graph is the same (plt.plot(freq,np.abs(S)), but the real part is different:

plt.plot(freq,S.real):

blue: no delay. green: with delay.

and in the time domain (I didn't center on 0):

more information:

How to do FFT fractional time delay (SOLVED)

https://antongrin.github.io/Seismika/Chapters/SignalProcessing/Ricker_wavelet.html

Answer 2

No time for full answer but in short, the wavelet positioning in time looks correct and is centered about $t=0$: the 0 is at index 0, and DFT is circular about it, so right half of the plot is actually negatives. ifftshift will visually center it, but this should not be done with FFT convolution.

Finding correct correspondence between continuous and discrete is the subject of periodization (accounting for aliasing, frame shift and size). Responding to another answer:

note that the output of an ifft is complex, so you should just plot the real part

Only do this if you've confirmed the imaginary part is zero; it's very easy to make discretization errors, and dropping the imaginary part will hide them.