Coding multidimensional de/convolution

Let us unroll the filter helix seen in Figure and see what we have. Start from the idea that a 2-D filter is generally made from a cluster of values near one another in two dimensions similar to the Laplacian operator in the figure. We see that in the helical approach, a 2-D filter is a 1-D filter containing some long intervals of zeros. The intervals are about the length of a 1-D seismogram.

Our program for 2-D convolution with a 1-D convolution program, could convolve with the somewhat long 1-D strip, but it is much more cost effective to ignore the many zeros, which is what we do. We do not multiply by the backside zeros, nor do we even store them in memory. Whereas an ordinary convolution program would do time shifting by a code line like iy=ix+lag, Module helicon ignores the many zero filter values on backside of the tube by using the code iy=ix+lag(ia) where a counter ia ranges over the nonzero filter coefficients. Before operator helicon is invoked, we need to prepare two lists, one list containing nonzero filter coefficients flt(ia), and the other list containing the corresponding lags lag(ia) measured to include multiple wraps around the helix. For example, the 2-D Laplace operator can be thought of as the 1-D filter  

$$\begin{array} {\vert r\vert r\vert r\vert r\vert r\vert r\ve... ... & 1 & \\ \hline 1 & -4 & 1\\ \hline & 1 & \\ \hline\end{array}$$ (6)

The first filter coefficient in equation (6) is +1 as implicit to module helicon. To apply the Laplacian on a $1000\times 1000$ mesh requires the filter inputs:

                i    lag(i)   flt(i)
               ---   ------   -----
                1      999      1
                2     1000     -4
                3     1001      1
                4     2000      1

Here we choose to use ``declaration of a type'', a modern computer language feature that is absent from Fortran 77. Fortran 77 has the built in complex arithmetic type. In module helix we define a type filter, actually, a helix filter. After making this definition, it will be used by many programs. The helix filter consists of three vectors, a real valued vector of filter coefficients, an integer valued vector of filter lags, and an optional vector that has logical values ``.TRUE.'' for output locations that will not be computed (either because of boundary conditions or because of missing inputs). The filter vectors are the size of the nonzero filter coefficents (excluding the leading 1.) while the logical vector is long and relates to the data size. The helix module allocates and frees memory for a helix filter. By default, the logical vector is not allocated but is set to null with the nullify operator and ignored. helixdefinition for helix-type filters

For those of you with no Fortran 90 experience, the ``%'' appearing in the helix module denotes a pointer. Fortran 77 has no pointers (or everything is a pointer). The C, C++, and Java languages use ``.'' to denote pointers. C and C++ also have a second type of pointer denoted by ``->''. The behavior of pointers is somewhat different in each language. Never-the-less, the idea is simple. In module helicon you see the expression aa%flt(ia). It refers to the filter named aa. Any filter defined by the helix module contains three vectors, one of which is named flt. The second component of the flt vector in the aa filter is referred to as aa%flt(2) which in the example above refers to the value 4.0 in the center of the laplacian operator. For data sets like above with 1000 points on the 1-axis, this value 4.0 occurs after 1000 lags, thus aa%lag(2)=1000.

Our first convolution operator tcai1 was limited to one dimension and a particular choice of end conditions. With the helix and Fortran 90 pointers, the operator helicon is a multidimensional filter with considerable flexibility (because of the mis vector) to work around boundaries and missing data. heliconhelical convolution The code fragment aa%lag(ia) corresponds to b-1 in tcai1 .

Operator helicon did the convolution job for Figure . As with tcai1 the adjoint of filtering is filtering backwards which means unscrewing the helix.

The companion to convolution is deconvolution. The module polydiv is essentially the same as polydiv1 , but here it was coded using our new filter type in module helix which will simplify our many future uses of convolution and deconvolution. Although convolution allows us to work around missing input values, deconvolution does not (any input affects all subsequent outputs), so polydiv never references aa%mis(ia). polydivhelical deconvolution