# Source code for bruges.filters.wavelets

```
"""
Seismic wavelets.
:copyright: 2021 Agile Geoscience
:license: Apache 2.0
"""
import warnings
from collections import namedtuple
import numpy as np
import scipy.signal
import matplotlib.pyplot as plt
from bruges.util import deprecated
def _get_time(duration, dt, sym=True):
"""
Make a time vector.
If `sym` is `True`, the time vector will have an odd number of samples,
and will be symmetric about 0. If it's False, and the number of samples
is even (e.g. duration = 0.016, dt = 0.004), then 0 will bot be center.
"""
# This business is to avoid some of the issues with `np.arange`:
# (1) unpredictable length and (2) floating point weirdness, like
# 1.234e-17 instead of 0. Not using `linspace` because figuring out
# the length and offset gave me even more of a headache than this.
n = int(duration / dt)
odd = n % 2
k = int(10**-np.floor(np.log10(dt)))
dti = int(k * dt) # integer dt
if (odd and sym):
t = np.arange(n)
elif (not odd and sym):
t = np.arange(n + 1)
elif (odd and not sym):
t = np.arange(n)
elif (not odd and not sym):
t = np.arange(n) - 1
t -= t[-1] // 2
return dti * t / k
def _generic(func, duration, dt, f, t=None, return_t=True, taper='blackman', sym=True):
"""
Generic wavelet generator: applies a window to a continuous function.
Args:
func (function): The continuous function, taking t, f as arguments.
duration (float): The length in seconds of the wavelet.
dt (float): The sample interval in seconds (often one of 0.001, 0.002,
or 0.004).
f (array-like): Dominant frequency of the wavelet in Hz. If a sequence is
passed, you will get a 2D array in return, one row per frequency.
t (array-like): The time series to evaluate at, if you don't want one
to be computed. If you pass `t` then `duration` and `dt` will be
ignored, so we recommend passing `None` for those arguments.
return_t (bool): If True, then the function returns a tuple of
wavelet, time-basis.
taper (str or function): The window or tapering function to apply.
To use one of NumPy's functions, pass 'bartlett', 'blackman' (the
default), 'hamming', or 'hanning'; to apply no tapering, pass
'none'. To apply your own function, pass a function taking only
the length of the window and returning the window function.
sym (bool): If True (default behaviour before v0.5) then the wavelet
is forced to have an odd number of samples and the central sample
is at 0 time.
Returns:
ndarray. wavelet(s) with centre frequency f sampled on t. If you
passed `return_t=True` then a tuple of (wavelet, t) is returned.
"""
if not return_t:
m = "return_t is deprecated. In future releases, return_t will always be True."
warnings.warn(m, DeprecationWarning, stacklevel=2)
f = np.asanyarray(f).reshape(-1, 1)
# Compute time domain response.
if t is None:
t = _get_time(duration, dt, sym=sym)
else:
if (duration is not None) or (dt is not None):
m = "`duration` and `dt` are ignored when `t` is passed."
warnings.warn(m, UserWarning, stacklevel=2)
t = np.array(t)
t[t == 0] = 1e-12 # Avoid division by zero.
f[f == 0] = 1e-12 # Avoid division by zero.
w = np.squeeze(func(t, f))
if taper is not None:
tapers = {
'bartlett': np.bartlett,
'blackman': np.blackman,
'hamming': np.hamming,
'hanning': np.hanning,
'none': lambda _: 1,
}
taper = tapers.get(taper, taper)
w *= taper(t.size)
if return_t:
Wavelet = namedtuple('Wavelet', ['amplitude', 'time'])
return Wavelet(w, t)
else:
return w
[docs]def sinc(duration, dt, f, t=None, return_t=True, taper='blackman', sym=True):
"""
sinc function centered on t=0, with a dominant frequency of f Hz.
If you pass a 1D array of frequencies, you get a wavelet bank in return.
.. plot::
import matplotlib.pyplot as plt
import bruges
w, t = bruges.filters.sinc(0.256, 0.002, 40)
plt.plot(t, w)
Args:
duration (float): The length in seconds of the wavelet.
dt (float): The sample interval in seconds (often one of 0.001, 0.002,
or 0.004).
f (array-like): Dominant frequency of the wavelet in Hz. If a sequence is
passed, you will get a 2D array in return, one row per frequency.
t (array-like): The time series to evaluate at, if you don't want one
to be computed. If you pass `t` then `duration` and `dt` will be
ignored, so we recommend passing `None` for those arguments.
return_t (bool): If True, then the function returns a tuple of
wavelet, time-basis.
taper (str or function): The window or tapering function to apply.
To use one of NumPy's functions, pass 'bartlett', 'blackman' (the
default), 'hamming', or 'hanning'; to apply no tapering, pass
'none'. To apply your own function, pass a function taking only
the length of the window and returning the window function.
Returns:
ndarray. sinc wavelet(s) with centre frequency f sampled on t. If
you passed `return_t = True` then a tuple of (wavelet, t) is returned.
"""
def func(t_, f_):
return np.sin(2*np.pi*f_*t_) / (2*np.pi*f_*t_)
return _generic(func, duration, dt, f, t, return_t, taper, sym=sym)
[docs]def cosine(duration, dt, f, t=None, return_t=True, taper='gaussian', sigma=None, sym=True):
"""
With the default Gaussian window, equivalent to a 'modified Morlet'
also sometimes called a 'Gabor' wavelet. The `bruges.filters.gabor`
function returns a similar shape, but with a higher mean frequancy,
somewhere between a Ricker and a cosine (pure tone).
If you pass a 1D array of frequencies, you get a wavelet bank in return.
.. plot::
import matplotlib.pyplot as plt
import bruges
w, t = bruges.filters.cosine(0.256, 0.002, 40)
plt.plot(t, w)
Args:
duration (float): The length in seconds of the wavelet.
dt (float): The sample interval in seconds (often one of 0.001, 0.002,
or 0.004).
f (array-like): Dominant frequency of the wavelet in Hz. If a sequence is
passed, you will get a 2D array in return, one row per frequency.
t (array-like): The time series to evaluate at, if you don't want one
to be computed. If you pass `t` then `duration` and `dt` will be
ignored, so we recommend passing `None` for those arguments.
return_t (bool): If True, then the function returns a tuple of
wavelet, time-basis.
taper (str or function): The window or tapering function to apply.
To use one of NumPy's functions, pass 'bartlett', 'blackman' (the
default), 'hamming', or 'hanning'; to apply no tapering, pass
'none'. To apply your own function, pass a function taking only
the length of the window and returning the window function.
sigma (float): Width of the default Gaussian window, in seconds.
Defaults to 1/8 of the duration.
Returns:
ndarray. sinc wavelet(s) with centre frequency f sampled on t. If
you passed `return_t=True` then a tuple of (wavelet, t) is returned.
"""
if sigma is None:
sigma = duration / 8
def func(t_, f_):
return np.cos(2 * np.pi * f_ * t_)
def taper(length):
return scipy.signal.gaussian(length, sigma/dt)
return _generic(func, duration, dt, f, t, return_t, taper, sym=sym)
[docs]def gabor(duration, dt, f, t=None, return_t=True, sym=True):
"""
Generates a Gabor wavelet with a peak frequency f0 at time t.
https://en.wikipedia.org/wiki/Gabor_wavelet
If you pass a 1D array of frequencies, you get a wavelet bank in return.
.. plot::
import matplotlib.pyplot as plt
import bruges
w, t = bruges.filters.gabor(0.256, 0.002, 40)
plt.plot(t, w)
Args:
duration (float): The length in seconds of the wavelet.
dt (float): The sample interval in seconds (often one of 0.001, 0.002,
or 0.004).
f (array-like): Centre frequency of the wavelet in Hz. If a sequence is
passed, you will get a 2D array in return, one row per frequency.
t (array-like): The time series to evaluate at, if you don't want one
to be computed. If you pass `t` then `duration` and `dt` will be
ignored, so we recommend passing `None` for those arguments.
return_t (bool): If True, then the function returns a tuple of
wavelet, time-basis.
Returns:
ndarray. Gabor wavelet(s) with centre frequency f sampled on t. If
you passed `return_t=True` then a tuple of (wavelet, t) is returned.
"""
def func(t_, f_):
return np.exp(-2 * f_**2 * t_**2) * np.cos(2 * np.pi * f_ * t_)
return _generic(func, duration, dt, f, t, return_t, sym=sym)
[docs]def ricker(duration, dt, f, t=None, return_t=True, sym=True):
"""
Also known as the mexican hat wavelet, models the function:
.. math::
A = (1 - 2 \pi^2 f^2 t^2) e^{-\pi^2 f^2 t^2}
If you pass a 1D array of frequencies, you get a wavelet bank in return.
.. plot::
import matplotlib.pyplot as plt
import bruges
w, t = bruges.filters.ricker(0.256, 0.002, 40)
plt.plot(t, w)
Args:
duration (float): The length in seconds of the wavelet.
dt (float): The sample interval in seconds (often one of 0.001, 0.002,
or 0.004).
f (array-like): Centre frequency of the wavelet in Hz. If a sequence is
passed, you will get a 2D array in return, one row per frequency.
t (array-like): The time series to evaluate at, if you don't want one
to be computed. If you pass `t` then `duration` and `dt` will be
ignored, so we recommend passing `None` for those arguments.
return_t (bool): If True, then the function returns a tuple of
wavelet, time-basis.
sym (bool): If True (default behaviour before v0.5) then the wavelet
is forced to have an odd number of samples and the central sample
is at 0 time.
Returns:
ndarray. Ricker wavelet(s) with centre frequency f sampled on t. If
you passed `return_t=True` then a tuple of (wavelet, t) is returned.
"""
if not return_t:
m = "return_t is deprecated. In future releases, return_t will always be True."
warnings.warn(m, DeprecationWarning, stacklevel=2)
f = np.asanyarray(f).reshape(-1, 1)
if t is None:
t = _get_time(duration, dt, sym=sym)
else:
if (duration is not None) or (dt is not None):
m = "`duration` and `dt` are ignored when `t` is passed."
warnings.warn(m, UserWarning, stacklevel=2)
pft2 = (np.pi * f * t)**2
w = np.squeeze((1 - (2 * pft2)) * np.exp(-pft2))
if return_t:
RickerWavelet = namedtuple('RickerWavelet', ['amplitude', 'time'])
return RickerWavelet(w, t)
else:
return w
[docs]def klauder(duration, dt, f,
autocorrelate=True,
t=None,
return_t=True,
taper='blackman',
sym=True,
**kwargs):
"""
By default, gives the autocorrelation of a linear frequency modulated
wavelet (sweep). Uses scipy.signal.chirp, adding dimensions as necessary.
.. plot::
import matplotlib.pyplot as plt
import bruges
w, t = bruges.filters.klauder(0.256, 0.002, [12, 48])
plt.plot(t, w)
Args:
duration (float): The length in seconds of the wavelet.
dt (float): is the sample interval in seconds (usually 0.001, 0.002,
or 0.004)
f (array-like): Upper and lower frequencies. Any sequence like (f1, f2).
A list of lists will create a wavelet bank.
autocorrelate (bool): Whether to autocorrelate the sweep(s) to create
a wavelet. Default is `True`.
t (array-like): The time series to evaluate at, if you don't want one
to be computed. If you pass `t` then `duration` and `dt` will be
ignored, so we recommend passing `None` for those arguments.
return_t (bool): If True, then the function returns a tuple of
wavelet, time-basis.
taper (str or function): The window or tapering function to apply.
To use one of NumPy's functions, pass 'bartlett', 'blackman' (the
default), 'hamming', or 'hanning'; to apply no tapering, pass
'none'. To apply your own function, pass a function taking only
the length of the window and returning the window function.
sym (bool): If True (default behaviour before v0.5) then the wavelet
is forced to have an odd number of samples and the central sample
is at 0 time.
**kwargs: Further arguments are passed to scipy.signal.chirp. They are
`method` ('linear','quadratic','logarithmic'), `phi` (phase offset
in degrees), and `vertex_zero`.
Returns:
ndarray: The waveform. If you passed `return_t=True` then a tuple of
(wavelet, t) is returned.
"""
if not return_t:
m = "return_t is deprecated. In future releases, return_t will always be True."
warnings.warn(m, DeprecationWarning, stacklevel=2)
if t is None:
t = _get_time(duration, dt, sym=sym)
else:
if (duration is not None) or (dt is not None):
m = "`duration` and `dt` are ignored when `t` is passed. "
m += "Pass None to suppress this warning."
warnings.warn(m, UserWarning, stacklevel=2)
t0, t1 = -duration/2, duration/2
f = np.asanyarray(f).reshape(2, -1)
f1, f2 = f
c = [scipy.signal.chirp(t, f1_+(f2_-f1_)/2., t1, f2_, **kwargs)
for f1_, f2_
in zip(f1, f2)]
if autocorrelate:
w = [np.correlate(c_, c_, mode='same') for c_ in c]
else:
w = c
w = np.squeeze(w) / np.amax(w)
if taper:
funcs = {
'bartlett': np.bartlett,
'blackman': np.blackman,
'hamming': np.hamming,
'hanning': np.hanning,
'none': lambda x: x,
}
func = funcs.get(taper, taper)
w *= func(t.size)
if return_t:
Sweep = namedtuple('Sweep', ['amplitude', 'time'])
return Sweep(w, t)
else:
return w
sweep = klauder
def _ormsby(t, f1, f2, f3, f4):
"""
Compute a single Ormsby wavelet. Private function.
"""
def numerator(f, t):
"""The numerator in the Ormsby equation."""
return (np.sinc(f * t)**2) * ((np.pi * f) ** 2)
pf43 = (np.pi * f4) - (np.pi * f3)
pf21 = (np.pi * f2) - (np.pi * f1)
w = ((numerator(f4, t)/pf43) - (numerator(f3, t)/pf43) -
(numerator(f2, t)/pf21) + (numerator(f1, t)/pf21))
return np.squeeze(w) / np.amax(w)
[docs]def ormsby(duration, dt, f, t=None, return_t=True, sym=True):
"""
The Ormsby wavelet requires four frequencies which together define a
trapezoid shape in the spectrum. The Ormsby wavelet has several sidelobes,
unlike Ricker wavelets.
.. plot::
import matplotlib.pyplot as plt
import bruges
w, t = bruges.filters.ormsby(0.256, 0.002, [5, 10, 40, 80])
plt.plot(t, w)
Args:
duration (float): The length in seconds of the wavelet.
dt (float): The sample interval in seconds (usually 0.001, 0.002,
or 0.004).
f (array-like): Sequence of form (f1, f2, f3, f4), or list of lists of
frequencies, which will return a 2D wavelet bank.
t (array-like): The time series to evaluate at, if you don't want one
to be computed. If you pass `t` then `duration` and `dt` will be
ignored, so we recommend passing `None` for those arguments.
return_t (bool): If True, then the function returns a tuple of
wavelet, time-basis.
sym (bool): If True (default behaviour before v0.5) then the wavelet
is forced to have an odd number of samples and the central sample
is at 0 time.
Returns:
ndarray: A vector containing the Ormsby wavelet, or a bank of them. If
you passed `return_t=True` then a tuple of (wavelet, t) is returned.
"""
if not return_t:
m = "return_t is deprecated. In future releases, return_t will always be True."
warnings.warn(m, DeprecationWarning, stacklevel=2)
try:
f = np.atleast_2d(f).reshape(-1, 4)
except ValueError as e:
raise ValueError("The last dimension of the frequency array must be of size 4.")
if t is None:
t = _get_time(duration, dt, sym=sym)
else:
if (duration is not None) or (dt is not None):
m = "`duration` and `dt` are ignored when `t` is passed."
warnings.warn(m, UserWarning, stacklevel=2)
w = np.squeeze([_ormsby(t, *fs) for fs in f])
if return_t:
OrmsbyWavelet = namedtuple('OrmsbyWavelet', ['amplitude', 'time'])
return OrmsbyWavelet(w, t)
else:
return w
def _ormsby_fft(duration, dt, f, P, sym):
fs = 1 / dt
fN = fs // 2
n = int(duration / dt)
a = map(lambda p: 10**(p/20), P)
# Linear interpolation of points.
x = np.linspace(0, int(fN), int(10*n))
xp = [ 0.] + list(f) + [fN]
fp = [0., 0.] + list(a) + [0., 0.]
W = np.interp(x, xp, fp)
# Compute inverse FFT.
w_ = np.fft.fftshift(np.fft.irfft(W))
L = int(w_.size // 2)
normalize = lambda d: d / np.max(abs(d))
return normalize(w_[L-n//2:L+n//2+sym])
[docs]def ormsby_fft(duration, dt, f, P=(0, 0), return_t=True, sym=True):
"""
Non-white Ormsby, with arbitary amplitudes.
Can use as many points as you like. The power of f1 and f4 is assumed to be 0,
so you only need to provide p2 and p3 (the corners). (You can actually provide
as many f points as you like, as long as there are n - 2 matching p points.)
.. plot::
import matplotlib.pyplot as plt
import bruges
w, t = bruges.filters.ormsby(0.256, 0.002, [5, 10, 40, 80])
plt.plot(t, w)
Args:
duration (float): The length in seconds of the wavelet.
dt (float): The sample interval in seconds (usually 0.001, 0.002,
or 0.004).
f (array-like): Sequence of form (f1, f2, f3, f4), or list of lists of
frequencies, which will return a 2D wavelet bank.
P (tuple): The power of the f2 and f3 frequencies, in relative dB.
(The magnitudes of f1 and f4 are assumed to be -∞ dB, i.e. a
magnitude of 0.) The default power values of (0, 0) results in a
trapezoidal spectrum and a conventional Ormsby wavelet. Pass, e.g.
(0, -15) for a 'pink' wavelet, with more energy in the lower
frequencies.
return_t (bool): If True, then the function returns a tuple of
wavelet, time-basis.
sym (bool): If True (default behaviour before v0.5) then the wavelet
is forced to have an odd number of samples and the central sample
is at 0 time.
Returns:
ndarray: A vector containing the Ormsby wavelet, or a bank of them. If
you passed `return_t=True` then a tuple of (wavelet, t) is returned.
"""
if not return_t:
m = "return_t is deprecated. In future releases, return_t will always be True."
warnings.warn(m, DeprecationWarning, stacklevel=2)
try:
f = np.atleast_2d(f).reshape(-1, 4)
except ValueError as e:
raise ValueError("The last dimension of the frequency array must be of size 4.")
w = np.squeeze([_ormsby_fft(duration, dt, fs, P, sym) for fs in f])
t = _get_time(duration, dt, sym=sym)
if return_t:
OrmsbyWavelet = namedtuple('OrmsbyWavelet', ['amplitude', 'time'])
return OrmsbyWavelet(w, t)
else:
return w
[docs]def berlage(duration, dt, f, n=2, alpha=180, phi=-np.pi/2, t=None, return_t=True, sym=True):
r"""
Generates a Berlage wavelet with a peak frequency f. Implements
.. math::
w(t) = AH(t) t^n \mathrm{e}^{- \alpha t} \cos(2 \pi f_0 t + \phi_0)
as described in Aldridge, DF (1990), The Berlage wavelet, GEOPHYSICS
55 (11), p 1508-1511. Berlage wavelets are causal, minimum phase and
useful for modeling marine airgun sources.
If you pass a 1D array of frequencies, you get a wavelet bank in return.
.. plot::
import matplotlib.pyplot as plt
import bruges
w, t = bruges.filters.berlage(0.256, 0.002, 40)
plt.plot(t, w)
Args:
duration (float): The length in seconds of the wavelet.
dt (float): The sample interval in seconds (often one of 0.001, 0.002,
or 0.004).
f (array-like): Centre frequency of the wavelet in Hz. If a sequence is
passed, you will get a 2D array in return, one row per frequency.
n (float): The time exponent; non-negative and real.
alpha(float): The exponential decay factor; non-negative and real.
phi (float): The phase.
t (array-like): The time series to evaluate at, if you don't want one
to be computed. If you pass `t` then `duration` and `dt` will be
ignored, so we recommend passing `None` for those arguments.
return_t (bool): If True, then the function returns a tuple of
wavelet, time-basis.
sym (bool): If True (default behaviour before v0.5) then the wavelet
is forced to have an odd number of samples and the central sample
is at 0 time.
Returns:
ndarray. Berlage wavelet(s) with centre frequency f sampled on t. If
you passed `return_t=True` then a tuple of (wavelet, t) is returned.
"""
if not return_t:
m = "return_t is deprecated. In future releases, return_t will always be True."
warnings.warn(m, DeprecationWarning, stacklevel=2)
f = np.asanyarray(f).reshape(-1, 1)
if t is None:
t = _get_time(duration, dt, sym=sym)
else:
if (duration is not None) or (dt is not None):
m = "`duration` and `dt` are ignored when `t` is passed."
warnings.warn(m, UserWarning, stacklevel=2)
H = np.heaviside(t, 0)
w = H * t**n * np.exp(-alpha * t) * np.cos(2 * np.pi * f * t + phi)
w = np.squeeze(w) / np.max(np.abs(w))
if return_t:
BerlageWavelet = namedtuple('BerlageWavelet', ['amplitude', 'time'])
return BerlageWavelet(w, t)
else:
return w
[docs]def generalized(duration, dt, f, u=2, t=None, return_t=True, imag=False, sym=True):
"""
Wang's generalized wavelet, of which the Ricker is a special case where
u = 2. The parameter u is the order of the time-domain derivative, which
can be a fractional derivative.
As given by Wang (2015), Generalized seismic wavelets. GJI 203, p 1172-78.
DOI: https://doi.org/10.1093/gji/ggv346. I am using the (more accurate)
frequency domain method (eq 4 in that paper).
.. plot::
import matplotlib.pyplot as plt
import bruges
w, t = bruges.filters.generalized(0.256, 0.002, 40, u=1.0)
plt.plot(t, w)
Args:
duration (float): The length of the wavelet, in s.
dt (float): The time sample interval in s.
f (float or array-like): The frequency or frequencies, in Hertz.
u (float or array-like): The fractional derivative parameter u.
t (array-like): The time series to evaluate at, if you don't want one
to be computed. If you pass `t` then `duration` and `dt` will be
ignored, so we recommend passing `None` for those arguments.
return_t (bool): Whether to return the time basis array.
center (bool): Whether to center the wavelet on time 0.
imag (bool): Whether to return the imaginary component as well.
sym (bool): If True (default behaviour before v0.5) then the wavelet
is forced to have an odd number of samples and the central sample
is at 0 time.
Returns:
ndarray. If f and u are floats, the resulting wavelet has duration/dt
= A samples. If you give f as an array of length M and u as an
array of length N, then the resulting wavelet bank will have shape
(M, N, A). If f or u are floats, their size will be 1, and they
will be squeezed out: the bank is always squeezed to its minimum
number of dimensions. If you passed `return_t=True` then a tuple
of (wavelet, t) is returned.
"""
if not return_t:
m = "return_t is deprecated. In future releases, return_t will always be True."
warnings.warn(m, DeprecationWarning, stacklevel=2)
# Make sure we can do banks.
f = np.asanyarray(f).reshape(-1, 1)
u = np.asanyarray(u).reshape(-1, 1, 1)
# Compute time domain response.
if t is None:
t = _get_time(duration, dt, sym=sym)
else:
if (duration is not None) or (dt is not None):
m = "`duration` and `dt` are ignored when `t` is passed."
warnings.warn(m, UserWarning, stacklevel=2)
dt = t[1] - t[0]
duration = len(t) * dt
# Basics.
om0 = f * 2 * np.pi
u2 = u / 2
df = 1 / duration
nyquist = (1 / dt) / 2
nf = 1 + nyquist / df
t0 = duration / 2
om = 2 * np.pi * np.arange(0, nyquist, df)
# Compute the spectrum from Wang's eq 4.
exp1 = np.exp((-om**2 / om0**2) + u2)
exp2 = np.exp(-1j*om*t0 + 1j*np.pi * (1 + u2))
W = (u2**(-u2)) * (om**u / om0**u) * exp1 * exp2
w = np.fft.ifft(W, t.size)
if not imag:
w = w.real
# At this point the wavelet bank has the shape (u, f, a),
# where u is the size of u, f is the size of f, and a is
# the number of amplitude samples we generated.
w_max = np.max(np.abs(w), axis=-1)[:, :, None]
w = np.squeeze(w / w_max)
if return_t:
GeneralizedWavelet = namedtuple('GeneralizedWavelet', ['amplitude', 'time'])
return GeneralizedWavelet(w, t)
else:
return w
```