2
$\begingroup$

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]')

enter image description here

enter image description here

$\endgroup$
2
  • $\begingroup$ you can just define your time array with time=np.arange(0,nsamples*dt,dt) and not type out each value, it would make your code look more compact $\endgroup$
    – XaC
    Jan 11, 2023 at 14:34
  • $\begingroup$ You're right, just edited the post. $\endgroup$
    – Luis Fraga
    Jan 12, 2023 at 8:47

2 Answers 2

0
$\begingroup$

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:

enter image description here

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

real FFT

blue: no delay. green: with delay.

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

more information:

How to do FFT fractional time delay (SOLVED)

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

$\endgroup$
1
  • $\begingroup$ Many thanks for your reponse and for all the graph visualisations, that helps a lot ! I definetely havo to find a good book on continuous and discrete TFT. As in my comment for @OverLordGoldDragon, my final goal will be to carry out a deconvolution on ameasured and then a convultion with a desired source signal. $\endgroup$
    – Luis Fraga
    Jan 12, 2023 at 9:05
1
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Thanks @OverLordGoldDragon, specially for the topic on periodization. My main goal is exatcly to do a convolution of the Ricker signal with an impulse response (like in this article : earthdoc.org/content/journals/10.3997/1873-0604.2003015). Maybe I will have more questions popping out soon. $\endgroup$
    – Luis Fraga
    Jan 12, 2023 at 8:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.