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.

 

#$
#$=head1 NAME
#$
#$createmshelix - create a multi-scale helix filter
#$
#$
#$=head1 SYNOPSIS
#$
#$C<aa=createmshelix(nd,center,gap,jump,na)>
#$
#$=head1 INPUT PARAMETERS
#$
#$=over 4
#$
#$=item nd      - C<integer(:)>
#$
#$      size of data
#$
#$=item center  - C<integer(:)>
#$
#$      location of the 1 of the filter within na box
#$
#$=item jump     - C<integer(:)>
#$
#$      stretches 
#$
#$=item gap     - C<integer(:)> 
#$
#$      distance along each axis before filter coef
#$
#$=item na      - C<integer(:)>
#$
#$      size of box arround filter
#$
#$=back
#$
#$=head1 RETURN VALUE
#$
#$=over 4
#$
#$=item aa       - C<type(nfilter)>
#$
#$      Helix filter
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$Create a multi-scale helix filter
#$
#$=head1 SEE ALSO
#$
#$L<createhelix>,L<mshelix>,L<createnhelix>,L<createmshelix>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut
#$
#$

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.

 

#$
#$=head1 NAME
#$
#$mshconest - convolution using multi-scale helix filters, adjoint is filter
#$  
#$=head1 SYNOPSIS
#$  
#$Initializer - C<call mshconest_init(x,aa)>
#$  
#$Operator    - C<ierr=mshconest_lop(adj,add,xx,yy)>
#$  
#$=head1 PARAMETERS
#$     
#$=over 4
#$  
#$=item x  - C<real(:)>
#$
#$      data
#$
#$=item aa - type(msfilter)
#$
#$      helix filter to perform convolution with
#$
#$=item adj,add,xx,yy -
#$
#$      standard operators parameters
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$  Masked multi-scale helix convolution, adjoint is the filter
#$
#$=head1 SEE ALSO
#$
#$L<mshelix>,L<hconest>,L<mshelicon>,L<nhconest>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$

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.  

#$=head1 NAME
#$
#$mspef  - find prediction error filter
#$
#$=head1 SYNOPSIS
#$
#$C<call find_pef(dd,aa,niter)>
#$
#$=head1 INPUT PARAMETERS
#$
#$=over 4
#$
#$=item dd  -  C<real(:)>
#$
#$      Input data 
#$
#$=item niter - integer 
#$
#$      Number of itterations
#$
#$=back  
#$
#$=head1 OUTPUT PARAMETERS
#$
#$=over 4
#$
#$=item  aa - type(msfilter)
#$ 
#$       output filter
#$
#$=back  
#$
#$=head1 DESCRIPTION
#$
#$ Find multi-scale prediction-error filter (helix magic)
#$
#$=head1 SEE ALSO
#$
#$L<npef>,L<hconest>,L<solver>,L<cgstep>,L<mshelix>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut

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.

 

#$
#$=head1 NAME
#$
#$mshelicon - convolution using mutli scale helix filters
#$  
#$=head1 SYNOPSIS
#$  
#$Initializer - C<call mshelicon_init(aa)>
#$  
#$Operator    - C<ierr=mshelicon_lop(adj,add,xx,yy)>
#$  
#$=head1 PARAMETERS
#$     
#$=over 4
#$  
#$=item aa - type(msfilter)
#$
#$      helix filter to perform convolution with
#$
#$=item adj,add,xx,yy - 
#$
#$      standard operators parameters
#$
#$=back 
#$
#$=head1 DESCRIPTION
#$
#$ Multi-scale onvolution, inverse to deconvolution.
#$ Requires the filter be causal with an implicit "1." at the onset.
#$
#$
#$=head1 SEE ALSO
#$
#$L<mshelix>,L<mshconest>,L<helicon>,L<nhelicon>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut
#$  
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 .

 

#$=head1 NAME
#$
#$msmis - Fill in missing data using a multi scale
#$
#$=head1 SYNOPSIS
#$
#$C<call mis1(niter,nx,ns,xx,aa,known,doprec)>
#$
#$=head1 PARAMETERS
#$
#$=over 4
#$
#$=item  niter - integer
#$
#$       number of iterations
#$
#$=item  nx    - integer
#$
#$       number of points in the data
#$
#$=item  ns    - integer  
#$
#$       number of scales
#$
#$=item  xx    -real        
#$
#$       fitting variable
#$ 
#$=item  aa    -type(filter) 
#$
#$       filter to apply
#$
#$=item  known -C<logical(:)>
#$
#$       Known data
#$
#$=item  doprec -logical
#$
#$       Whether or not to run preconditioning
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$fill in missing data by minimizing power out of a given multi scale filter
#$ by helix magic works in any number of dimensions
#$
#$=head1 SEE ALSO
#$  
#$L<nmis2>,L<mis2>,L<mask1>,L<helicon>,L<solver>,L<mshelicon>
#$
#$=head1 LIBRARY
#$   
#$B<geef90>
#$  
#$=cut

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()
  }
}