(13) | ||

(14) |

# vv = (aa Dx - pp Dt) uu # subroutine wavekill1( aa, pp, m1,n1,n2 , uu, vv) integer conj, i1,i2, m1,n1,n2 real aa, pp, s11, s12, s21, s22, uu(m1,n2), vv(m1,n2-1) s11 = -aa+pp; s12 = aa+pp s21 = -aa-pp; s22 = aa-pp call zero( m1*(n2-1), vv) do i1=1, n1-1 { # vv is one point shorter than uu on both axes. do i2=1, n2-1 { vv(i1,i2) = vv(i1,i2) + uu(i1 ,i2) * s11 + uu(i1 ,i2+1) * s12 + uu(i1+1,i2) * s21 + uu(i1+1,i2+1) * s22 }} return; end

In a finite difference representation,
the difference operator is a two-dimensional filter
that fits on a differencing star.
As a matrix operation,
this two dimensional convolution is denoted .(More details can be found in my previous book ``''.)
The program `wavekill1()` applies the operator
which can be specialized the operators
,
,.I carefully arranged the side boundaries so the filter never
runs off the sides of the data.
Thus the output is shorter than the input by one point
both the *t*-axis and the *x*-axis.
The reason for using these side boundaries is so that large data
sets can be chopped into independent sections
without the boundaries themselves
affecting the result.
By chopping a large data set into sections,
we can handle curved events as piecewise linear.

When only one wave is present and the data is adequately sampled,
then finding the best value of *p* is a single-parameter,
linear least-squares problem.
Let be a vector whose components are values of
taken on a mesh in (*t*,*x*).
Likewise let contain .Since we wish ,we minimize the quadratic function of *p*

(15) |

(16) |

(17) |

(18) |

# measure coherency and dip, and compute residual res = (Dx - p Dt) uu # subroutine puck ( m1, n1, n2, uu, coh, pp, res ) integer i1, i2, m1, n1,n2 real uu(m1,n2), res(m1,n2), xx, xt, tt, coh, pp temporary real dx(m1,n2-1), dt(m1,n2-1) call wavekill1( 1., 0., m1,n1,n2 , uu, dx) call wavekill1( 0.,-1., m1,n1,n2 , uu, dt) xx = 1.e-30; tt = 1.e-30; xt = 0. do i1=1,n1-1 { do i2=1,n2-1 { xt = xt + dt(i1,i2) * dx(i1,i2) tt = tt + dt(i1,i2) * dt(i1,i2) xx = xx + dx(i1,i2) * dx(i1,i2) }} coh = sqrt((xt/tt)*(xt/xx)) pp = xt/tt call wavekill1( 1., pp, m1,n1,n2 , uu, res) return; end

Finally, an undesirable feature of equation
(13)
is that the derivatives roughen the data so the
spectrum of and is not that of .What I am doing now is putting the data through
leaky integration on the time axis before I start.
I am planning to switch this from *t* to *x*
so the operator can be thought of as .Then the residual should be in phase with the input signal
so I can subtract it from the data leaving ``cleaned up'' data.

1/13/1998