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.

The multiscale helix filter is defined in module mshelix , analogous to the single-scale module helix . A new function onescale extracts our usual helix filter of one particular scale from the multiscale filter.

 

module mshelix {                                # multiscale helix filter type
  use helix
  type msfilter {             
    real,    dimension( :),    pointer :: flt   # (nh) filter coefficients
    integer, dimension( :, :), pointer :: lag   # (nh,ns) filter (lags,scales)
    logical, dimension( :, :), pointer :: mis   # (nd,ns) boundary conditions
  }
contains
  subroutine msallocate( msaa, nh, ns) {
    type( msfilter)   :: msaa
    integer           :: nh, ns
    allocate( msaa%flt( nh), msaa%lag( nh, ns))
    msaa%flt = 0.;  nullify( msaa%mis)
    }
  subroutine msdeallocate( msaa) {
    type( msfilter) :: msaa
    deallocate( msaa%flt, msaa%lag)
    if( associated( msaa%mis)) deallocate( msaa%mis)
    }
  subroutine onescale( i, msaa, aa) {      # Extract single-scale filter.
    integer, intent (in) :: i
    type( filter)        :: aa
    type( msfilter)      :: msaa
    aa%flt => msaa%flt
    aa%lag => msaa%lag( :, i)
    if( associated( msaa%mis)) 
	aa%mis => msaa%mis( :, i)
    else
        nullify( aa%mis)
    }
}

We create a multscale helix with module createmshelixmod . An expanded scale helix filter is like an ordinary helix filter except that the lags are scaled according to a jump.

 

module createmshelixmod {         # Create multiscale helix filter lags and mis
use mshelix
use createhelixmod
use bound
contains
  function createmshelix( nd, center, gap, jump, na)  result( msaa) {
    type( msfilter)                   :: msaa    # needed by mshelicon.
    integer, dimension(:), intent(in) :: nd, na  # data and filter axes
    integer, dimension(:), intent(in) :: center  # normally (na1/2,na2/2,...,1)
    integer, dimension(:), intent(in) :: gap     # normally ( 0,   0,  0,...,0)
    integer, dimension(:), intent(in) :: jump	 # jump(ns) stretch scales
    type( filter)                     :: aa
    integer                           :: is, ns, nh, n123 
    aa = createhelix( nd, center, gap, na)
    ns = size( jump);  nh = size( aa%lag);  n123 = product( nd)
    call msallocate( msaa, nh, ns)
    do is = 1, ns	
        msaa%lag(:,is) = aa%lag(:)*jump(is)	 # set lags for expanded scale
    call deallocatehelix( aa)
    allocate( msaa%mis( n123, ns))
    do is = 1, ns {				 # for all scales
     call onescale( is, msaa, aa);  nullify( aa%mis)	# extract a filter
     call boundn( nd, nd, na*jump(is), aa)		# set up its boundaries
     msaa%mis( :, is) = aa%mis;  deallocate( aa%mis)	# save them
    }
  }
}

First we examine code for estimating a prediction-error filter that is applicable at many scales. We simply invoke the usual filter operator hconest for each scale.

 

module mshconest {      # multi-scale helix convolution, adjoint is the filter.
  use mshelix
  use hconest
  use helix
  integer, private                   :: nx, ns
  real,    dimension(:),     pointer :: x
  type( msfilter)                    :: msaa
#% _init( x, msaa)
   nx = size( x);   ns = size( msaa%lag, 2)
#% _lop( a(:), y(nx,ns))
    integer  :: is, stat1
    type (filter) :: aa
    do is = 1, ns {  
  	call onescale (is, msaa, aa)
       	call hconest_init( x, aa)
       	stat1 = hconest_lop( adj, .true., a, y(:,is))
    }
}

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 helix filter module hconest by a call to mshconest.  

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

The purpose of pack(dd,.true.) is to produce the one-dimensional array expected by solver() .

Similar code applies to the operator in (4) which is needed for missing data problems. This is like mshconest except the adjoint is not the filter but the input.

 

module mshelicon {                                  # Multi-scale convolution
  use mshelix
  use helicon
  integer           :: nx, ns
  type( msfilter)   :: msaa
#%  _init (nx, ns, msaa)
#%  _lop ( xx( nx),  yy( nx, ns))
  integer :: is, stat1
  type (filter) :: aa
  do is = 1, ns { 
     call onescale( is, msaa, aa)
     call helicon_init( aa)
     stat1 = helicon_lop( adj, .true., xx, yy(:,is))
     }
}

The multiscale missing-data module msmis2 is just like the usual missing-data module mis2 except that the filtering is done with the multiscale filter mshelicon .

 

module msmis2 {                      # multi-scale missing data interpolation
  use mshelicon
  use cgstep_mod
  use solver_mod
contains
  subroutine mis1( niter, nx, ns, xx, aa, known) {
    integer,                intent( in)     :: niter, nx, ns
    logical, dimension( :), intent( in)     :: known
    type( msfilter),        intent( in)     :: aa
    real,    dimension( :), intent( in out) :: xx
    real,    dimension( nx*ns)              :: dd
    dd = 0.
    call mshelicon_init( nx,ns, aa)
    call solver( mshelicon_lop, cgstep, niter= niter, x = xx, dat = dd, 
                    known = known, x0 = xx)
    call cgstep_close()
  }
}