LMO by nearest-neighbor interpolation

To do linear moveout (LMO) correction, we need to time-shift data. Shifting data requires us to interpolate it. The easiest interpolation method is the nearest-neighbor method. We begin with a signal given at times t = t0+dt*(it-1) where it is an integer. Then we can use equation (3), namely $\tau=t-px$.Given the location tau of the desired value we backsolve for an integer, say itau. In Fortran, conversion of a real value to an integer is done by truncating the fractional part of the real value. To get rounding up as well as down, we add 0.5 before conversion to an integer, namely itau=int(1.5+(tau-tau0)/dt). This gives the nearest neighbor. The way the program works is to identify two points, one in (t,x)-space and one in $(\tau,x)$-space. Then the data value at one point in one space is carried to the other. The adjoint operation copies $\tau$ space back to t space. The subroutine used in the illustrations above is lmo() with adj=1.  

# linear moveout
#
subroutine lmo( adj,add, slow, tau0, t0,dt, x0,dx, modl,nt,nx,  data       )
integer         adj,add,                                nt,nx,       it,ix,iu
real t, x, tau,          slow, tau0, t0,dt, x0,dx, modl(nt,nx), data(nt,nx)
call adjnull(   adj,add,                           modl,nt*nx,  data,nt*nx)
do ix= 1, nx {  x= x0 + dx * (ix-1)
do it= 1, nt {  t= t0 + dt * (it-1)
                tau =  t - x * slow
                 iu = 1.5001 + (tau-tau0)/dt
                if( 0 < iu  &&  iu <= nt)
                        if( adj == 0 )
                                data(it,ix) = data(it,ix) + modl(iu,ix)
                        else
                                modl(iu,ix) = modl(iu,ix) + data(it,ix)
                }}
return; end

Nearest neighbor rounding is crude but ordinarily very reliable. I discovered a very rare numerical roundoff problem peculiar to signal time-shifting, a problem which arises in the linear moveout application when the water velocity, about 1.48km/sec is approximated by 1.5=3/2. The problem arises only where the amount of the time shift is a numerical value (like 12.5000001 or 12.499999) and the fractional part should be exactly 1/2 but numerical rounding pushes it randomly in either direction. We would not care if an entire signal was shifted by either 12 units or by 13 units. What is troublesome, however, is if some random portion of the signal shifts 12 units while the rest of it shifts 13 units. Then the output signal has places which are empty while adjacent places contain the sum of two values. Linear moveout is the only application where I have ever encountered this difficulty. A simple fix here was to modify the lmo() subroutine changing the ``1.5'' to ``1.5001''. The problem disappears if we use a more accurate sound velocity or if we switch from nearest-neighbor interpolation to linear interpolation.