Editing Deconvolution (section)

===Seismology===

The concept of deconvolution had an early application in [[reflection seismology]].  In 1950, [[Enders Robinson]] was a graduate student at [[MIT]].  He worked with others at MIT, such as [[Norbert Wiener]], [[Norman Levinson]], and economist [[Paul Samuelson]], to develop the "convolutional model" of a reflection [[seismogram]].  This model assumes that the recorded seismogram ''s''(''t'') is the convolution of an Earth-reflectivity function ''e''(''t'') and a [[seismic]] [[wavelet]] ''w''(''t'') from a [[point source]], where ''t'' represents recording time. Thus, our convolution equation is

:<math>s(t) = (e * w)(t). \, </math>

The seismologist is interested in ''e'', which contains information about the Earth's structure.  By the [[convolution theorem]], this equation may be [[Fourier transform]]ed to

: <math>S(\omega) = E(\omega)W(\omega) \, </math>

in the [[frequency domain]], where <math>\omega</math> is the frequency variable. By assuming that the reflectivity is white, we can assume that the [[Spectral density|power spectrum]] of the reflectivity is constant, and that the power spectrum of the seismogram is the spectrum of the wavelet multiplied by that constant.  Thus,

: <math>|S(\omega)| \approx k|W(\omega)|. \, </math>

If we assume that the wavelet is [[minimum phase]], we can recover it by calculating the minimum phase equivalent of the power spectrum we just found.  The reflectivity may be recovered by designing and applying a [[Wiener filter]] that shapes the estimated wavelet to a [[Dirac delta function]] (i.e., a spike).  The result may be seen as a series of scaled, shifted delta functions (although this is not mathematically rigorous):

: <math>e(t)=\sum_{i=1}^N r_i\delta(t-\tau_i),</math>

where ''N'' is the number of reflection events, <math>r_i</math> are the [[reflection coefficient]]s, <math>t-\tau_i</math> are the reflection times of each event, and <math>\delta</math> is the [[Dirac delta function]].

In practice, since we are dealing with noisy, finite [[Bandwidth (computing)|bandwidth]], finite length, [[Sampling (signal processing)|discretely sampled]] datasets, the above procedure only yields an approximation of the filter required to deconvolve the data.  However, by formulating the problem as the solution of a [[Toeplitz matrix]] and using [[Levinson recursion]], we can relatively quickly estimate a filter with the smallest [[mean squared error]] possible.  We can also do deconvolution directly in the frequency domain and get similar results.  The technique is closely related to [[linear prediction]].