Coding the multiscale filter operator

Equation (3) shows an example where the first output signal is the ordinary one and the second output signal used a filter interlaced with zeros. We prepare subroutines that allow for more output signals, each with its own filter interlace parameter given in the table jump(ns). Each entry in the jump table corresponds to a particular scaling of a filter axis. The number of output signals is ns and the number of zeros interlaced between filter points for the j-th signal is jump(j)-1.

First we examine code for estimating a prediction-error filter that is applicable at many scales. The coding is much like that in hconest , the new ingredient being the loop over scales and the jumping along axes.  

module mshconest {      # multi-scale helix convolution, adjoint is the filter.
  integer, private                   :: nx, ns
  real,    dimension( :),    pointer :: x
  logical, dimension( :, :), pointer :: bad	# output depends on bad inputs
  integer, dimension( :),    pointer :: lag, jump
#% _init( bad, x, lag, jump)
   nx = size( x); ns = size( jump)
#% _lop( a (:), y (nx,ns))
    integer  ia, ix, iy, is
    do is = 1, ns {
    do ia = 1, size( a) {
         do iy = 1  + lag( ia)*jump( is), nx {    if( bad( iy, is)) cycle  
            ix = iy - lag( ia)*jump( is)
             if( adj) 
                    a( ia) += y( iy, is) * x( ix)
             else
                    y( iy, is) += a( ia) * x( ix)
         }
    }}
}

The multiscale prediction-error filter finding subroutine is nearly identical to the usual subroutine find_pef() . (That routine cleverly ignores missing data while estimating a PEF.) To easily extend pef to multiscale filters we replace its call to the ordinary 2-D filter subroutine by a call to mshconest. Of course, it also needs to know that multiscale filtering has a larger output space (with more panels) than single-scale filtering does.  

module mspef {  # Find multi-scale prediction-error filter (helix magic)
  use mshconest
  use cgstep_mod
  use solver_mod
contains
  subroutine find_pef( yy, aa, lag, jump, mx, niter) {
    integer,                   intent( in)  :: niter
    real,    dimension( :),    pointer      :: yy
    logical, dimension( :, :), pointer      :: mx
    integer, dimension( :),    pointer      :: lag, jump
    real,    dimension( :),    intent( out) :: aa
    integer                                 :: is, nx
    real,    dimension( size( yy), size( jump)) :: dd
    nx = size( yy)
    do is = 1, size( jump)
       dd (:,is) = -yy
    call mshconest_init( mx, yy, lag, jump)
    call solver( mshconest_lop, cgstep, aa, pack(dd, .true.), niter)
    call cgstep_close()
  }
}

Likewise, I coded up the operator in (4).  

module mshelicon {                                  # Multi-scale Convolution
#                          Designed for gapped filters (helical 2-D, 3-D, etc).
#            Requires the filter be causal with an implicit "1." at the onset.
integer                         :: nx, ns
real,    dimension (:), pointer :: aa
integer, dimension (:), pointer :: lag, jump
#%  _init (nx, ns, aa, lag, jump)
if( any( lag <= 0)) call erexit ("Filter is not causal")
#%  _lop ( xx (nx),  yy (nx,ns))
integer iy, ix, ia, is
do is= 1, ns {
if( adj)              # zero lag
	xx += yy (:,is)
else
	yy (:,is) += xx
do ia= 1, size(aa) {
do ix= 1, nx-lag(ia)*jump(is) {  iy = ix + lag(ia)*jump(is)
        if( adj)
                        xx(ix) += yy(iy,is) * aa(ia)
	else	
                        yy(iy,is) += xx(ix) * aa(ia)
        }}}
}

The multiscale missing-data module msmis2 is just like the usual missing-data module mis2 except that the filtering is done with a multiscale subroutine and the size of the residual must be precomputed externally.

 

module msmis2 { # multi-scale missing data interpolation
  use mshelicon
  use cgstep_mod
  use solver_mod
contains
# fill in missing data by minimizing power out of a given filter.
# by helix magic works in any number of dimensions
  subroutine mis1( niter, xx, aa, lag, jump, known) {
    integer,                intent( in)     :: niter
    logical, dimension( :), intent( in)     :: known
    real,    dimension( :), intent( in out) :: xx
    real,    dimension( :), pointer         :: aa
    integer, dimension( :), pointer         :: lag, jump
    logical, dimension( :), allocatable     :: mm
    real,    dimension( :), allocatable     :: yy, dd
    integer                                 :: nx, ny, pad, ns
    nx = size( xx); ns = size( jump)
    pad = maxval( lag)*maxval( jump); ny = nx + 2*pad
    allocate( yy( ny), mm( ny), dd( ny*ns)); dd = 0.
    yy( :2*pad) = 0.;      yy( 2*pad+1:) = xx
    mm( :2*pad) = .false.; mm( 2*pad+1:) = known
    call mshelicon_init( ny, ns, aa, lag, jump)
    call solver( mshelicon_lop, cgstep, niter= niter, x = yy, dat = dd, 
                    known = mm, x0 = yy)
    call cgstep_close( )
    where( xx == 0.)  xx = yy( 2*pad+1:)
    deallocate( yy, mm, dd)
  }
}