Coding 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 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  

$$\left[ \begin{array} {ccc} & 1 & \\ 1 & -4 & 1\\ & 1 &\end... ...itions} (1,\,0, \; ... \; 0,\,1,\,-4,\,1,\,0, \; ... \; 0,\,1).$$ (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.  

#$
#$=head1 NAME
#$
#$helix - module containing allocate and deallocate of a  helix filter
#$
#$=head1 SYNOPSIS
#$  
#$C<call allocatehelix(aa,nh)>
#$  
#$C<call deallocatehelix(aa)>
#$  
#$=head1 INPUT PARAMETERS
#$  
#$=over 4
#$  
#$=item aa - type(helix)   
#$
#$      Filter
#$
#$=item nh - integer
#$
#$      Number of coefs in filter
#$
#$=back  
#$
#$=head1 DESCRIPTION
#$
#$Allocate and deallocation of helix filter.
#$
#$=head1 COMMENTS
#$
#$type(filter) :
#$
#$=over 4
#$
#$=item flt - C<real(:)> (nh)
#$
#$      filter coefficients
#$
#$=item lag - C<real(:)> (nh)
#$
#$      filter lags
#$    
#$=item mis - C<real(:)>  (nd)
#$
#$      boundary conditions
#$  
#$=back 
#$
#$=head1 SEE ALSO
#$
#$L<nhelix>,L<mshelix>
#$ 
#$=head1 LIBRARY 
#$    
#$B<geef90>
#$
#$=cut

module helix {                                  # DEFINE helix filter type
  type filter {             
    real,    dimension( :), pointer :: flt      # (nh) filter coefficients
    integer, dimension( :), pointer :: lag      # (nh) filter lags
    logical, dimension( :), pointer :: mis      # (nd) boundary conditions
  }
contains
  subroutine allocatehelix( aa, nh ) {		# allocate a filter
    type( filter) :: aa
    integer       :: nh				# count of filter coefs (excl 1)
    allocate( aa%flt( nh), aa%lag( nh))		# allocate filter and lags.
    nullify( aa%mis)				# set null pointer for "mis".
    aa%flt = 0.					# zero filter coef values
  }
  subroutine deallocatehelix( aa) {		# destroy a filter
    type( filter) :: aa
    deallocate( aa%flt, aa%lag)			# free memory
    if( associated( aa%mis))			# if logicals were allocated
        deallocate( aa%mis)			# free them
  }
}

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.  

#$
#$=head1 NAME
#$
#$helicon - convolution using helix filters
#$
#$=head1 SYNOPSIS
#$  
#$Initializer - C<call helicon_init(aa)>
#$  
#$Operator    - C<ierr=helicon_lop(adj,add,xx,yy)> 
#$  
#$=head1 PARAMETERS
#$  
#$=over 4
#$  
#$=item aa - type(filter)  
#$
#$      helix filter to perform convolution with
#$
#$=item adj,add,xx,yy - 
#$
#$      standard operators parameters
#$
#$=back  
#$
#$=head1 DESCRIPTION
#$
#$ Convolution, inverse to deconvolution.
#$ Requires the filter be causal with an implicit "1." at the onset.
#$
#$
#$=head1 SEE ALSO
#$
#$L<helix>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut  
#$
module helicon {                       # Convolution, inverse to deconvolution.
#             Requires the filter be causal with an implicit "1." at the onset.
use helix
type( filter) :: aa
#%  _init( aa)
#%  _lop ( xx,  yy)
integer iy, ix, ia
if( adj)              # zero lag
	xx += yy
else
	yy += xx
do ia = 1, size( aa%lag) {
      do iy = 1  + aa%lag( ia), size( yy) {    	
	    if( associated( aa%mis)) { if( aa%mis( iy)) cycle}
            ix = iy - aa%lag( ia)
            if( adj)
                        xx(ix) += yy(iy) * aa%flt(ia)
	    else	
                        yy(iy) += xx(ix) * aa%flt(ia)
	    }
      }
}

The code fragment aa%lag(ia) corresponds to ib-1 in tcai1 .

Operator helicon did the convolution job for Figure 1. 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).  

#$
#$=head1 NAME
#$
#$polydiv - polynomial division
#$ 
#$=head1 SYNOPSIS 
#$
#$Initializer - C<call polydiv_init(nd,aa)>
#$
#$Operator    - C<ierr=polydiv_lop(adj,add,xx,yy)>
#$
#$=head1 PARAMETERS 
#$
#$=over 4
#$
#$=item nd - integer
#$
#$      number of data points
#$
#$=item aa - type(filter)
#$
#$      helix filter to perform convolution with
#$
#$=item adj,add,xx,yy -
#$
#$      standard operators parameters
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$ Polynomial division (deconvolution), inverse to convolution.
#$ Requires the filter be causal with an implicit "1." at the onset.
#$
#$
#$=head1 SEE ALSO
#$
#$L<helix>,L<hconest>,L<helicon>,L<npolydiv>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut
module polydiv {			  # Helix polynomial division
use helix
integer                           :: nd
type( filter)                     :: aa
real, dimension (nd), allocatable :: tt
#%  _init ( nd, aa)
#%  _lop  ( xx, yy)
integer  ia, ix, iy
tt = 0.
if( adj) {
	do ix= nd, 1, -1 {
		tt( ix) = yy( ix)
		do ia = 1, size( aa%lag) {
			iy = ix + aa%lag( ia);     if( iy > nd)  next
			tt( ix) -=  aa%flt( ia) * tt( iy)
			} 
		}
	xx += tt
 } else { 
	do iy= 1, nd {
		tt( iy) = xx( iy)
		do ia = 1, size( aa%lag) {
			ix = iy - aa%lag( ia);      if( ix < 1)  next
			tt( iy) -=  aa%flt( ia) * tt( ix)
			} 
		}
	yy += tt
        }
}

EXERCISES:

1. Observe the matrix matrix () which corresponds to subroutine tcai1 . What is the matrix corresponding to helicon ?