Finite differences on a helix

The function  

$$\bold r = (-1,0,\cdots,0,-1,4,-1,0, \cdots,0,-1)$$ (4)

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$: 

$$\bold h \quad = \quad(1.791,-.651,-.044,-.024,\cdots,\cdots,-.044,-.087,-.200,-.558 )$$ (5)

Wind the signal $\bold r$ around a vertical-axis helix to see its two-dimensional shape $\bold R$: 

$$\bold R \quad = \quad \left[ \begin{array} {rrr} & -1 & \\ -1 & 4 & -1 \\ & -1 & \end{array}\right]$$ (6)

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$: 

$$\bold H \quad = \quad \left[ \begin{array} {rrrrrrrrr} & ... ...ots &-.044 & -.087 & -.200 & -.558 & & & & \end{array} \right]$$ (7)

In the 2-D representation (7) we see the coefficients diminishing rapidly away from maximum value 1.791. My claim is that the 2-D autocorrelation of (7) is (6). 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 9 is that $\bold H$ enhances drainage patterns whereas ${\partial\over\partial y}$ enhances mountain ridges.

 

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

Our construction makes $\bold H$ have the energy spectrum k_x²+k_y², 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 (7) is nonunique in that causality associated with the helix mapping can be defined along either x or y axes; thus the operator (7) 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 10) 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
Figure 10 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'$.

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 10 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.