Convolution from scratch

Processing Prerequisites

Learning objectives

  • Define discrete convolution and compute it by the flip-slide-multiply-sum recipe
  • Recognize the convolutional model of a seismic trace: s(t)=w(t)∗r(t)+n(t)s(t) = w(t) \ast r(t) + n(t)
  • See why convolving a spike with a wavelet reproduces the wavelet, and why close spikes tune
  • Name the three algebraic properties of convolution (commutative, associative, distributive) and why they matter in practice

A seismic trace is not what the earth reflects; it is what the earth reflects convolved with the wavelet your source emitted. Every processing step that tries to recover the earth’s reflectivity (deconvolution and impedance inversion above all) has to undo this convolution; even a migrated image is still reflectivity convolved with the wavelet. So we start by understanding what convolution actually does.

1. The convolutional model

s(t)  =  w(t) ∗ r(t)  +  n(t)s(t) \;=\; w(t) \,\ast\, r(t) \;+\; n(t)

where

  • s(t)s(t) is the recorded seismic trace,
  • w(t)w(t) is the source wavelet (plus receiver response, coupling, everything band-limited by acquisition),
  • r(t)r(t) is the earth’s reflectivity, a spike train at every impedance contrast,
  • n(t)n(t) is noise (ambient and coherent),
  • the ast\\ast symbol is convolution.

This one line is the forward model for almost everything we do. Deconvolution is the inverse problem: given ss, recover rr.

2. Discrete convolution, the recipe

In the digital world, signals are sequences. The discrete convolution of an input x\[n\] with a kernel h\[n\] is defined as

(x \\ast h)\[\\tau\] \\;=\\; \\sum\_{n=-\\infty}^{\\infty}\\, x\[n\]\\,\\cdot\\,h\[\\tau - n\]

That looks intimidating. In practice it is a four-step mechanical procedure you can always fall back on:

  1. Flip the kernel: h\[k\] becomes h\[-k\].
  2. Slide the flipped kernel to position tau\\tau.
  3. Multiply the flipped-and-shifted kernel, sample by sample, against xx.
  4. Sum the products. That number is the output at time tau\\tau.

Advance tau\\tau by one sample and repeat: the whole output trace is built one sample at a time. In Figure 0.3 the input is a reflectivity r\[n\] sampled every 2 ms and the kernel is a wavelet ww. You slide the output time tau\\tau, see which spikes the flipped wavelet covers, check the sum of the products against the trace, and then bring two spikes together to see what their copies of the wavelet do to each other.

Convolution: input ⊛ filter = outputInput pulse⊛Filter response=Output (smoothed)Interactive figure, enable JavaScript to design custom filters and watch the output shape change.

3. What to notice in the figure

The figure opens on two spikes of opposite sign, +0.20 at 60 ms and −0.20 at 72 ms, with a 30 Hz Ricker wavelet and tau\\tau = 60 ms. The flipped wavelet covers both spikes, and the sum of its two products, 0.20times1.000+(−0.20)times(−0.434)=0.2870.20 \\times 1.000 + (-0.20) \\times (-0.434) = 0.287, is exactly the sample of the trace at 60 ms.

  • One spike: at most one product is ever nonzero, so as tau\\tau moves the trace spells out the wavelet sample by sample: a copy of the wavelet placed at the spike, centred on it for the zero-phase Ricker and starting at it for a causal (minimum-phase) wavelet, whose largest lobe arrives 28 ms later. Convolution with a unit impulse returns the kernel itself; this is why the kernel is called the impulse response.
  • The flip: the Ricker is symmetric, so its flip in (d) lies exactly on it and the flip step is invisible. Switch to the minimum-phase wavelet and the flipped copy in (a) leans the other way; forget the flip and the trace would start before the spike.
  • Two spikes of opposite sign are the top and base of a thin bed. Their copies overlap, and each copy’s side lobe adds to the other’s main lobe: at 12 ms apart the pair peaks at 1.4 times one spike, and the peak is largest at 14 ms, close to the rule of thumb 1/(2.6f_mathrmp)1/(2.6 f\_{\\mathrm p}) = 12.8 ms for a Ricker. This is tuning: near a quarter of a wavelength the bed looks brighter than either reflection. Thinner than that the copies cancel (at 2 ms apart the pair peaks at 36 % of one spike, and at 0 ms it is gone), so a very thin bed looks weaker than it is.
  • Two spikes of the same sign first partly cancel, down to 56 % of one spike at 14 ms, where each copy’s trough falls on the other’s peak, then merge into one broader, stronger event (1.8 times at 4 ms, 2.0 times at 0 ms) that no trace can split into two interfaces.
  • Blocky log: six interfaces of an impedance log, each smeared into a copy of the wavelet; at 30 Hz the copies interfere and the trace peaks at 1.3 times its largest interface alone. Deconvolution is trying to undo this smearing.

Tuning (plate (e)) is set by the wavelet, not the geology: raise f_mathrmpf\_{\\mathrm p} and the tuning spacing shrinks with 1/f_mathrmp1/f\_{\\mathrm p}, from 16 ms at 25 Hz to 6 ms at 60 Hz, so a bed 20 ms thick that reads 1.3 times too bright at 25 Hz reads at its true amplitude at 60 Hz.

4. Properties you will use constantly

Three properties make convolution pleasant to reason about:

  • Commutative: xasth=hastxx \\ast h = h \\ast x. You can think of either signal as “the kernel.”
  • Associative: (xasth_1)asth_2=xast(h_1asth_2)(x \\ast h\_1) \\ast h\_2 = x \\ast (h\_1 \\ast h\_2). Two filters in series equal one combined filter: the source signature, ghosts, attenuation and the recording instrument cascade into the single effective wavelet w(t)w(t).
  • Distributive: xast(h_1+h_2)=xasth_1+xasth_2x \\ast (h\_1 + h\_2) = x \\ast h\_1 + x \\ast h\_2. This is linearity: each reflection is smeared independently and the trace is the sum of the copies, which is why the pair of spikes in Figure 0.3 is the sum of two single-spike traces.

5. Length arithmetic

If xx has length NN and hh has length MM, then xasthx \\ast h has length N+M−1N + M - 1. In Figure 0.3 the reflectivity has 80 samples and the wavelet 41, so the trace has 120: the wavelet smears each input sample across MM output samples, past both ends of the input. This matters when you design filters: a long operator produces long transients at the trace ends, which is why traces are tapered before filtering.

6. Special cases worth memorizing

  • Convolution with a delta \\delta\[n\] is the identity: xastdelta=xx \\ast \\delta = x.
  • Convolution with a shifted delta shifts: (x \\ast \\delta\[\\cdot - k\])\[n\] = x\[n-k\].
  • Convolution with a rectangular pulse of length LL and height 1/L1/L is a moving average of width LL.
  • Convolution in the time domain is multiplication in the frequency domain (proved in Section 0.4). This is the single fact that makes FFT-based processing fast enough to be practical.
The one sentence to remember

Convolution smears each input sample into the shape of the kernel. A seismic trace is the earth’s reflectivity smeared by the source wavelet, and every inverse operation in processing is some flavor of unsmearing.

Where this goes next

Section 0.4 proves that this “smearing” operation has a staggeringly simpler form once we pass to the frequency domain: convolution in time becomes multiplication in frequency. That one theorem is why the FFT sits at the heart of virtually every processing flow.

References

  • Bracewell, R. N. (1999). The Fourier Transform and Its Applications (3rd ed.). McGraw-Hill.
  • Oppenheim, A. V., Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
  • Yilmaz, Ö. (2001). Seismic Data Analysis (2 vols.). SEG.
  • Widess, M. B. (1973). How thin is a thin bed? Geophysics, 38(6), 1176–1180.
  • Claerbout, J. F. (1976). Fundamentals of Geophysical Data Processing. McGraw-Hill.

This page is prerendered for SEO and accessibility. The interactive widgets above hydrate on JavaScript load.