PROGRAM FOR MULTIDIMENSIONAL DE/CONVOLUTION

Let us unroll the filter helix seen in Figure 2 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 as in the figure.

Notice a ``1'' at the filter's onset (a common convention) and other values damping to zero a little distance away from the ``1''. Unrolling this filter into a 1-D strip, we find the ``1'' on the end of the strip, a few other values nearby, a long sequence of zeros, zeros that were on the back side of the coil, and eventually we come upon some more small values, values that are close to the ``1'' in the next column of the two-dimensional space, but are far away in one-dimensional space. Thus, in the helical approach, a 2-D filter is a 1-D filter containing some long intervals of zeros. The intervals are the length of a 1-D seismogram.

My 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 I do. I do not multiply by the backside zeros, nor do I even store them in memory. Subroutine helicon() did the convolution job for Figure 1. It looks like a one-dimensional convolution program except where the input is shifted from the output. Let a counter ia range over the nonzero filter coefficients. Whereas an ordinary convolution program would do time shifting by a code line like ix=iy-lag, subroutine helicon() ignores the many zero filter values on backside of the tube by using the code ix=iy-lag(ia). Before subroutine helicon() is invoked, we need to prepare two lists, one list containing nonzero filter coefficients aa(ia), and the other list containing the corresponding lags lag(ia) measured to include multiple wraps around the helix.

You can ignore the second half of subroutine helicon() which does adjoint filtering. Because I solve so many least-squares problems, in my textbooks I always accompany every linear operator by its adjoint (matrix transpose). As you might notice, the adjoint of filtering is filtering backwards. The ends of the 1-D axis are a little trickier than you might imagine, but there is no need to examine them now. (One end is transient and the other end is truncated.)

# Polynomial multiplication and its adjoint.
#      Designed for gapped filters (helical 2-D, 3-D, etc).
#      Requires the filter be causal with an implicit "1." at the onset.
#
subroutine helicon( adjoint, aa,na, lag,     xx,n123,  yy   )
integer iy, ix, ia, adjoint,    na, lag(na),    n123
real                         aa(na),         xx(n123), yy(n123)
do ia = 1, na
        if( lag(ia) <= 0 ) call erexit('Filter is non causal')
if( adjoint == 0)
        do iy = 1, n123 {
                        yy(iy) = xx(iy)
                do ia = 1, na {                        ix = iy - lag(ia)
                                                   if( ix < 1 ) break
                        yy(iy) = yy(iy) + aa(ia) * xx( ix)
                        }
                }
else
        do ix = n123, 1, -1  {
                        xx(ix) = yy(ix)
                do ia = 1, na {                        iy = ix + lag(ia)
                                                   if( iy > n123 ) break
                        xx(ix) = xx(ix) + aa(ia) * yy( iy)
                        }
                }
return; end