Chapter formulates the basic principle of missing data interpolation as follows:
A method for restoring missing data is to ensure that the restored data, after specified filtering, has minimum energy.Mathematically, this principle can be expressed by the simple equation
(26) |
(27) |
(28) |
(29) |
If D is represented by a discrete convolution, the natural choice for P is the corresponding deconvolution operator:
(30) |
The following module, mis2 , does the actual interpolation.
module mis2 {
use helicon
use polydiv
use mask1
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, known, precon, trans) {
logical, intent( in) :: precon, trans
integer, intent( in) :: niter
logical, dimension( :), intent( in) :: known
real, dimension( :), intent( in out) :: xx
real, dimension( :), pointer :: aa
integer, dimension( :), pointer :: lag
logical, dimension( :), pointer :: mm
real, dimension( :), allocatable :: yy, dd
integer :: nx, ny, pad, n1, n2
nx = size( xx); pad = maxval( lag); ny = nx + 2*pad
allocate( yy( ny), mm( ny)); yy = 0.; mm = .false.
if (trans) { n1 = 1; n2 = nx # emulate transient boundaries
mm( ny-pad+1:ny) = .true. # zero values beyond end.
} else { n1 = 2*pad+1; n2 = ny}
yy (n1:n2) = xx; mm (n1:n2) = known
if( precon) { # KP p = K y, m = P p
call mask1_init( mm)
call polydiv_init( aa, lag)
call solver_prec( mask1_lop, cgstep, niter= niter, x= yy, dat= yy,
prec= polydiv_lop, nprec= ny, eps= 0.)
} else { # KA m = 0
allocate( dd( ny)); dd = 0.
call helicon_init( aa, lag)
call solver( helicon_lop, cgstep, niter= niter, x = yy, dat = dd,
known = mm, x0 = yy)
deallocate( dd)
}
call cgstep_close( )
where( xx == 0.) xx = yy(n1:n2)
deallocate( yy, mm)
}
}
The operator mask1, apearing in mis2 . is a simple masked assignment (copy under mask.)
module mask1 { # masking operator
logical, dimension( :), pointer :: m
#% _init( m)
#% _lop( x, y)
if( adj)
where( m) x += y
else #
where( m) y += x
}