Contents

Ricker wavelet and its spectrum

Introduction 🎍

A Ricker wavelet is a signal commonly used in seismic interpretation, modelling and inversion. It has a relatively simple visual form with a positive peak and two negative side-lobes.

Sometimes it is referred to as Mexican hat wavelet.

Sombrero from Collins Dictionary

Fig. 1 Ricker wavelet / Mexican hat / sombrero

Analytically, it has a following representation in time domain:

\[A(t) = (1-2 \pi^2 f^2 t^2) e^{-\pi^2 f^2 t^2}\]

See also

Objectives 🎯

  1. Generate a Ricker wavelet given its dominant frequency

  2. Plot its time response and frequency spectrum

  3. Apply phase rotation to the wavelet

Generate a Ricker wavelet 🆕

First, we need to set the peak dominant frequency of the wavelet (\(f\) term in the formula), let’s say 30 Hz. Since typical frequencies of seismic signals fall between 10 and 100 Hz, that would be appropriate.

freq = 30

Second, let’s define a time vector (array). As the wavelet reaches highest amplitude when time is zero, it is convenient to have a zero-centered time vector.

Normally, in seismic exploration we work with time sampling rate of 0.5, 1, or 2 milliseconds.

Here, let the time-vector span from -0.05 to 0.05 seconds with a step of 1 ms.

import numpy as np

dt = 0.001
tmin = -0.05
tmax = 0.05
t = np.arange(tmin, tmax, dt)
NSamples = len(t)

Time and frequency having been set up, we can proceed to computation of amplitudes.

from numpy import pi, exp

w = (1 - 2*(pi*freq*t)**2) * exp(-(pi*freq*t)**2)

import matplotlib.pyplot as plt

plt.plot(t, w)
plt.xlabel('Time, s')
plt.ylabel('Amplitude')
plt.show()
../../_images/Ricker_wavelet_16_0.png

Compute its spectrum 🌈

If the spectrum of the wavelet is \(Z\) (array of complex numbers), then its amplitude spectrum is \(|Z|\) and phase spectrum is \(arg(Z)\):

\[ Z = |Z| \cdot e^{j \cdot arg(Z)} \]

Fourier spectrum can be computed by Discrete Fourier Transform (fft function in scipy.fftpack module). We can also easily determine its amplitude and phase by applying abs and angle functions of numpy.

from scipy.fftpack import fft, ifft, fftfreq

S = fft(w)
AS = np.abs(S)
PS = np.angle(S)

plt.plot(AS, label = 'Phase spectrum')
plt.plot(PS, label = 'Amplitude spectrum')
../../_images/Ricker_wavelet_22_0.png

What do we see here?

First, the amplitude and phase spectra are symmetric with respect to the center. We covered this in another page.

Second, the amplitude spectrum has a clear maximum at a certain point. We can expect this point to correspond to the dominant frequency of the wavelet. This can be verified by computing FFT frequencies (via fftfreq, for example) and making another plot. We omit the “mirrored” part of the spectra here, zooming-in on range from 0 to 250 Hz.

f = fftfreq(NSamples, dt)

plt.plot(f, AS, label = 'Amplitude spectrum')
plt.plot(f, PS, label = 'Phase spectrum')
plt.xlim(0, 250)
plt.axvline(freq, lw = 1, ls = '--', label = 'Dominant frequency')
../../_images/Ricker_wavelet_26_0.png

Third, the phase spectrum is non-zero, which is somewhat surprising! We did not expect this to happen since the amplitude at zero-time is maximum and the wavelet is perfectly symmetric.

Then why does the phase differ from zero? How does the true zero-phase wavelet look like?

Zero-Phase Ricker wavelet ⚖️

The apparent contradiction mentioned above comes from the fact that the time-vector did not count in FFT.

Indeed, to compute the spectrum we simply passed variable w to the fft function, ignoring time-vector t.

Thus, the above spectra were the spectra of the following signal (pay attention to the x-axis):

../../_images/Ricker_wavelet_30_0.png

This wavelet is actually not zero-phase! It has a maximum at the middle-sample (at time 0.05 s) rather than at 0.

So let’s try and create the new Ricker wavelet, with phase equal to 0.

In order to achieve this, we have to:

  1. Leave its amplitude unchanged;

  2. Replace phase spectrum by zeros;

  3. Take the inverse Fourier transform of the new spectrum;

For that, we will multiply the phase spectrum PS of wavelet w by zero.

AS0 = AS
PS0 = np.zeros_like(PS)

plt.plot(f, AS0, label = 'Amplitude spectrum')
plt.plot(f, PS0, label = 'Phase spectrum')
plt.xlim(0, 250)
../../_images/Ricker_wavelet_35_0.png 
w0 = np.real_if_close(ifft(AS0 * exp(1j*PS0)))

w0 is the new wavelet obtained by zeroing the phase of original Ricker wavelet.

Therefore, denoting zero-phase spectrum by Z0:

\[\begin{split} \begin{eqnarray*} & arg(Z_0) = 0, \\ & Z_0 = |Z_0| \cdot e^{j \cdot arg(Z_0)} = |Z_0| \end{eqnarray*} \end{split}\]

We expect Z0 to be real.

What is np.real_if_close and why we need it here

Generally, the wavelet resulting from inverse FFT is complex. In our particular case, we expect the wavelet to be real (because its spectrum is symmetric). However, the imaginary part of new wavelet is not perfectly zero due to computational issues. But it is so small (order of 1*10-17) that can be neglected with no negative effect. For such sutuations, numpy has an implemented function numpy.real_if_close. It converts complex-but-almost-real numbers to real.

plt.plot(t-np.min(t), w0)
plt.xlabel('Time, s')
plt.ylabel('Amplitude')
../../_images/Ricker_wavelet_42_0.png

In fact, this operation circularly shifted (or rolled) the wavelet by NSamples//2 samples. This signal appears to be ugly at first glance, but!

Recalling the periodicity implied in Discrete Fourier Transform, the new wavelet actually looks like this:

../../_images/Ricker_wavelet_44_0.png

And despite its strange appearance, it is exactly what we need! It is a wavelet with a maximum at 0 sample, symmetric with respect to y-axis. And it is the periodicity of Fourier transforms that allows us to construct the zero-phase wavelet.

This effect resembles the periodicity of discrete Fourier spectrum. As with spectra, the ‘mirrored’ right part of the wavelet can be virtually shifted to the negative part of time axis.

../../_images/Ricker_wavelet_46_0.png

In fact, this is only a matter of visualization. Of course, for most purposes (convolution, for instance) it is preferable to use a normal-looking signal. By the way, the ‘True zero-phase’ and ‘Centered maximum’ wavelets can be related using the shift theorem, a property of Fourier transform.

However, in terms of discrete Fourier transform, always remember how the zero-phase wavelet looks like:

../../_images/Ricker_wavelet_48_0.png

Constant-phase Ricker wavelet ➖

As an experiment, why not introduce a constant phase shift to the wavelet? Let’s say, 60 degrees (which is \(\pi/3\)).

Here, we should not forget about the symmetry of the phase spectrum which must be an odd function for real signals. Hence we have to add 60° to ‘positive’ frequencies and subtract 60° from so-called ‘negative’ frequencies. We start with a ‘normal’ good-looking wavelet (with maximum in the middle).

AS60 = AS
PS60 = np.zeros_like(PS)
PS60[1:NSamples//2] = PS[1:NSamples//2] + pi/3
PS60[NSamples//2+1:] = PS[NSamples//2+1:] - pi/3
w60 = np.real_if_close(ifft(AS0 * exp(1j*PS60)))
../../_images/Ricker_wavelet_52_0.png

A little bit of playing with phase yields the following picture:

PS_sh = np.zeros_like(PS)
AS_sh = np.zeros_like(AS)
phase_step = 45
colors = plt.cm.rainbow(np.linspace(0, 1, 360//phase_step))


for i, phshift in enumerate(range(-180, 180, phase_step)):
    AS_sh = AS0
    PS_sh[1:NSamples//2] = PS[1:NSamples//2] + pi*phshift/180
    PS_sh[NSamples//2+1:] = PS[NSamples//2+1:] - pi*phshift/180
    w_sh = np.real_if_close(ifft(AS_sh * exp(1j*PS_sh)))
    plt.plot(t, w_sh, c = colors[i], label = '{0:.0f}'.format(phshift))
    
plt.axhline(0, lw=0.5, c = '#333333', zorder = 0)
plt.axvline(0, lw=0.5, c = '#333333', zorder = 0)
plt.xlabel('Time, s')
plt.ylabel('Amplitude')
plt.legend(loc = 'upper right')

plt.show()
../../_images/Ricker_wavelet_54_0.png

The imaginary enveloping line that is tangent to each wavelet of this set is called Instantaneous Amplitude or Amplitude Envelope.

previous

Seismika | Intro

next

Complex trace attributes