Theoretically, a double integration of the second derivative
gives the desired pulse.
A representation in the discrete time domain is
the product of (1-*Z*)^{2} with
, which is 1.
Double integration
amounts to spectral division
by .Mathematically the problem is that vanishes at .In practice the problem is that dividing by where it is small
amplifies noises at those low frequencies.
(Inversion theorists are even more frustrated
because they are trying to create something like
a velocity profile, roughly a step function,
and they need to do something like a third integration.)
Old nuts like this illustrate
the dichotomy between theory and practice.

Chapter provides a theoretical solution to this problem in the Fourier domain. Here we will express the same concepts in the time domain. Define as follows:

Withy_{t}Given data.b_{t}Known filter.x_{t}Excitation (to be found). Noise: data minus filtered excitation.

(50) | ||

(51) |

(52) |

(53) |

subroutine ident( adj, add, epsilon, n, pp, qq ) integer i, adj, add, n real epsilon, pp(n), qq(n) # equivalence (pp,qq) OK if( adj == 0 ) { if( add == 0 ) { do i=1,n { qq(i) = epsilon * pp(i) } } else { do i=1,n { qq(i) = qq(i) + epsilon * pp(i) } } } else { if( add == 0 ) { do i=1,n { pp(i) = epsilon * qq(i) } } else { do i=1,n { pp(i) = pp(i) + epsilon * qq(i) } } } return; end

We can use any convolution routine we like,
but for simplicity, I selected `contrunc()` so the output
would be the same length as the input.
The two operators `ident()` and `contrunc()`
could be built into a new operator.
I found it easier to simply cascade them
in the deghosting subroutine `deghost()` below.

# deghost: min |rrtop| = | y - bb (contrunc) xx | # x |rrbot| | 0 - epsilon I xx | subroutine deghost( eps, nb,bb, n, yy, xx, rr, niter) integer iter, nb, n, niter real bb(nb), yy(n), eps # inputs. typically bb=(1,-2,1) real xx(n), rr(n+n) # outputs. temporary real dx(n), sx(n), dr(n+n), sr(n+n) call zero( n, xx) call copy( n, yy, rr(1 )) # top half of residual call zero( n , rr(1+n)) # bottom of residual do iter= 0, niter { call contrunc(1,0,1,nb,bb, n,dx,n,rr); call ident(1,1,eps, n,dx,rr(1+n)) call contrunc(0,0,1,nb,bb, n,dx,n,dr); call ident(0,0,eps, n,dx,dr(1+n)) call cgstep( iter, n,xx,dx,sx, _ n+n,rr,dr,sr) } return; end

10/21/1998