-
Notifications
You must be signed in to change notification settings - Fork 97
/
plot_shift.py
executable file
·92 lines (78 loc) · 2.87 KB
/
plot_shift.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
r"""
Shift
=====
This example shows how to use the :py:class:`pylops.signalprocessing.Shift`
operator to apply fractional delay to an input signal. Whilst this operator
acts on 1D signals it can also be applied on any multi-dimensional signal on
a specific direction of it.
"""
import numpy as np
import matplotlib.pyplot as plt
import pylops
plt.close('all')
###############################################################################
# Let's start with a 1D example. Define the input parameters: number of samples
# of input signal (``nt``), sampling step (``dt``) as well as the input
# signal which will be equal to a ricker wavelet:
nt = 127
dt = 0.004
t = np.arange(nt) * dt
ntwav = 41
wav = pylops.utils.wavelets.ricker(t[:ntwav], f0=20)[0]
wav = np.pad(wav, [0, nt-len(wav)])
WAV = np.fft.rfft(wav, n=nt)
###############################################################################
# We can shift this wavelet by :math:`5.5*dt`:
shift = 5.5 * dt
Op = pylops.signalprocessing.Shift(nt, shift, sampling=dt,
real=True, dtype=np.float64)
wavshift = Op * wav
wavshiftback = Op.H * wavshift
plt.figure(figsize=(10, 3))
plt.plot(t, wav, 'k', lw=2, label='Original')
plt.plot(t, wavshift, 'r', lw=2, label='Shifted')
plt.plot(t, wavshiftback, '--b', lw=2, label='Adjoint')
plt.axvline(t[ntwav-1], color='k')
plt.axvline(t[ntwav-1] + shift, color='r')
plt.xlim(0, .3)
plt.legend()
plt.title('1D Shift')
plt.tight_layout()
###############################################################################
# We can repeat the same exercise for a 2D signal and perform the shift
# along the first and second dimensions.
shift = 10.5 * dt
# 1st dir
wav2d = np.outer(wav, np.ones(10))
Op = pylops.signalprocessing.Shift((nt, 10), shift, dir=0, sampling=dt,
real=True, dtype=np.float64)
wav2dshift = (Op * wav2d.ravel()).reshape(nt, 10)
wav2dshiftback = (Op.H * wav2dshift.ravel()).reshape(nt, 10)
fig, axs = plt.subplots(1, 3, figsize=(10, 3))
axs[0].imshow(wav2d, cmap='gray')
axs[0].axis('tight')
axs[0].set_title('Original')
axs[1].imshow(wav2dshift, cmap='gray')
axs[1].set_title('Shifted')
axs[1].axis('tight')
axs[2].imshow(wav2dshiftback, cmap='gray')
axs[2].set_title('Adjoint')
axs[2].axis('tight')
fig.tight_layout()
# 2nd dir
wav2d = np.outer(wav, np.ones(10)).T
Op = pylops.signalprocessing.Shift((10, nt), shift, dir=1, sampling=dt,
real=True, dtype=np.float64)
wav2dshift = (Op * wav2d.ravel()).reshape(10, nt)
wav2dshiftback = (Op.H * wav2dshift.ravel()).reshape(10, nt)
fig, axs = plt.subplots(1, 3, figsize=(10, 3))
axs[0].imshow(wav2d, cmap='gray')
axs[0].axis('tight')
axs[0].set_title('Original')
axs[1].imshow(wav2dshift, cmap='gray')
axs[1].set_title('Shifted')
axs[1].axis('tight')
axs[2].imshow(wav2dshiftback, cmap='gray')
axs[2].set_title('Adjoint')
axs[2].axis('tight')
fig.tight_layout();