next up previous print clean
Next: Matrix view of the Up: Multidimensional recursive filters via Previous: FEATURES OF 1-D THAT

FINITE DIFFERENCES ON A HELIX

The function  
 \begin{displaymath}
\bold r = (-1,0,\cdots,0,-1,4,-1,0, \cdots,0,-1)\end{displaymath} (7)
is an autocorrelation function. It is symmetrical about the ``4'' and its Fourier transform is positive for all frequencies. Digging out an old textbook Claerbout (1976), we discover how to compute a causal wavelet with this autocorrelation. I used the ``Kolmogoroff spectral-factorization method'' to find this wavelet $\bold h$: 
 \begin{displaymath}
\bold h \quad = \quad(1.791,-.651,-.044,-.024,\cdots,\cdots,-.044,-.087,-.200,-.558 )\end{displaymath} (8)
Wind the signal $\bold r$ around a vertical-axis helix to see its two-dimensional shape $\bold R$: 
 \begin{displaymath}
\bold R \quad = \quad
\left[
 \begin{array}
{rrr}
 & -1 & \\  -1 & 4 & -1 \\  & -1 &
 \end{array}\right]\end{displaymath} (9)
This 2-D filter is is the negative of the finite-difference representation of the negative of the Laplacian operator, generally denoted $\nabla^2= {\partial^2\over \partial x^2} + {\partial^2\over \partial y^2} $.Now wind the signal $\bold h$ around the same helix to see its two-dimensional shape $\bold H$: 
 \begin{displaymath}
\bold H \quad = \quad
 \left[
 \begin{array}
{rrrrrrrrr}
 & ...
 ...ots &-.044 & -.087 & -.200 & -.558 & & & &
 \end{array} \right]\end{displaymath} (10)
In the 2-D representation (10) we see the coefficients diminishing rapidly away from maximum value 1.791. My claim is that the 2-D autocorrelation of (10) is (9). You verified this idea at the beginning of this paper where the numbers were all ones. You can check it again in a few moments if you drop the small values, say 0.2 and smaller.

Since the autocorrelation of $\bold H$ is $\bold H'\bold H =\bold R=-\nabla^2$is a second derivative, the operator $\bold H$ must be something like a first derivative. As a geophysicist, I found it natural to compare the operator ${\partial\over\partial y}$to $\bold H$ by applying them to a local topographic map. The result shown in Figure 5 is that $\bold H$ enhances drainage patterns whereas ${\partial\over\partial y}$ enhances mountain ridges.

 
helocut
helocut
Figure 5
(a) Topography, (b) helical derivative, (c) slope south.


view

Our construction makes $\bold H$ have the energy spectrum kx2+ky2, so the magnitude of its Fourier transform is $\sqrt{\vert k_x^2+k_y^2\vert}$.This rotationally invariant function in the Fourier domain contrasts sharply with the nonrotationally invariant function shape in (x,y)-space. The difference must arise from the phase spectrum. The factorization (10) is nonunique in that causality associated with the helix mapping can be defined along either x or y axes; thus the operator (10) can be rotated or reflected.

The operator $\bold H$ has curious similarities and differences with the familiar gradient and divergence operators. In two-dimensional physical space, the gradient maps one field to two fields (north slope and east slope). The factorization of $-\nabla^2$ with the helix gives us the operator $\bold H$that maps one field to one field. Being a one-to-one transformation (unlike gradient and divergence) the operator $\bold H$ is potentially invertible by deconvolution (recursive filtering).

I have chosen the name[*] ``helix derivative'' or ``helical derivative'' for the operator $\bold H$.A telephone pole has a narrow shadow behind it. The helix integral (middle frame of Figure 6) and the helix derivative (left frame) show shadows with an angular bandwidth approaching $180^\circ$.Thus, $\bold H$ is much less directional than either ${\partial\over\partial x}$ or ${\partial\over\partial y}$.

 
lapfac
lapfac
Figure 6
Deconvolution by a filter whose autocorrelation is the two-dimensional Laplacian operator. Amounts to solving the Poisson equation. Left is $\bold q$; Middle is $\bold q/\bold H$; Right is $(\bold q/\bold H)/\bold H'$.


view burn build edit restore

This is where the story all comes together. One-dimensional theory, such as the Kolmogoroff spectral factorization, produces not merely a causal wavelet with the required autocorrelation. It produces a wavelet that is stable in deconvolution. Using $\bold H$ in one-dimensional polynomial division, we can solve many formerly difficult problems very rapidly. Consider the Laplace equation with sources (Poisson's equation). Polynomial division and its reverse (adjoint) gives us $\bold p =(\bold q/\bold H)/\bold H'$which means that we have solved $\ \ \nabla^2 \bold p = -\bold q\ \ $ in Figure 6 by using polynomial division on a helix. Using the seven coefficients shown, the cost is fourteen multiplications (because we need to run both ways) per mesh point.



 
next up previous print clean
Next: Matrix view of the Up: Multidimensional recursive filters via Previous: FEATURES OF 1-D THAT
Stanford Exploration Project
10/23/1998