Interlacing a filter

The filter below can be designed despite alternate missing traces. This filter destroys plane waves. If the plane wave should happen to pass halfway between the ``d'' and the ``e'', those two points could interpolate the halfway point, at least for well-sampled temporal frequencies, and the time axis should always be well sampled. For example, d=e=-.5 would almost destroy the plane wave and it is an aliased planewave for its higher frequencies.  

$$\begin{array} {ccccccccc} a &\cdot &b &\cdot &c &\cdot &d &... ...&\cdot &\cdot &\cdot &1 &\cdot &\cdot &\cdot &\cdot \end{array}$$ (1)

We could use module pef to find the filter (1), if we set up the lag table lag appropriately. Then we could throw away alternate zeroed rows and columns (rescale the lag) to get the filter  

$$\begin{array} {ccccc} a &b &c &d &e \\ \cdot &\cdot &1 &\cdot &\cdot \end{array}$$ (2)

which could be used with subroutine mis1() , to find the interleaved data because both the filters (1) and (2) have the same dip characteristics.

Figure 1 shows three plane waves recorded on five channels and the interpolated data.

 

lace390
Figure 1 Left is five signals, each showing three arrivals. With the data shown on the left (and no more), the signals have been interpolated. Three new traces appear between each given trace, as shown on the right.

Both the original data and the interpolated data can be described as ``beyond aliasing,'' because on the input data the signal shifts exceed the signal duration. The calculation requires only a few seconds of a two-stage least-squares method, in which the first stage estimates a PEF (inverse spectrum) of the known data, and the second uses the PEF to estimate the missing traces. Figure 1 comes from PVI which introduces the clever method described above. We will review how that was done and examine the F90 codes that generalize it to N-dimensions. Then we'll go on to more general methods that allow missing data in any location. Before the methods of this section are applied to field data for migration, data must be broken into many overlapping tiles of size about like those shown here and the results from each tile pieced together. That is described later in chapter .

A PEF is like a differential equation. The more plane-wave solutions you expect, the more lags you need on the data. Returning to Figure 1, the filter must cover four traces (or more) to enable it to predict three plane waves. In this case, na=(9,4). As usual, the spike on the 2-D PEF is at center=(5,1). We see the filter is expanded by a factor of jump=4. The data size is nd=(75,5) and gap=0. Before looking at the code lace for estimating the PEF, it might be helpful to recall the basic utilities line2cart and cart2line for conversion between a multidimensional space and the helix filter lag. We need to sweep across the whole filter and ``stretch'' its lags on the 1-axis. We do not need to stretch its lags on the 2-axis because the data has not yet been interlaced by zero traces.  

module lace {                                     # find PEF on interlaced data
  use createhelixmod 
  use bound
  use pef
  use cartesian
contains
  function lace_pef( dd, jump, nd, center, gap, na)  result( aa) {
    type( filter)                      :: aa
    integer,               intent( in) :: jump
    integer, dimension(:), intent( in) :: nd, center, gap, na
    real,    dimension(:), pointer     :: dd          # input data
    integer, dimension(:), pointer     :: savelags    # holding place
    integer, dimension( size( nd))     :: ii
    integer                            :: ih, nh, lag0, lag
    aa = createhelix( nd, center, gap, na);  nh = size( aa%lag)
    savelags => aa%lag;   allocate( aa%lag( nh))      # prepare interlaced helix
    call cart2line(  na, center, lag0)
    do ih = 1, nh {				 	# Sweep thru the filter.
       call line2cart( na, ih+lag0, ii)
       ii = ii - center;            ii(1) = ii(1)*jump  # Interlace on 1-axis.
       call cart2line( nd,  ii+1, lag)
       aa%lag( ih) = lag - 1
       }
    call boundn( nd, nd, (/ na(1)*jump, na(2:) /), aa)  # Define  aa.mis
    call find_pef( dd, aa, nh*2)                        # Estimate aa coefs
    deallocate( aa%lag);   aa%lag => savelags		# Restore filter lags
  }
}

The line ii(1)=ii(1)*jump means we interlace the 1-axis but not the 2-axis because the data has not yet been interlaced with zero traces. For a 3-D filter aa(na1,na2,na3), the somewhat obtuse expression (/na(1)*jump, na(2:)/) is a three component vector containing (na1*jump, na2, na3).

After the PEF has been found, we can get missing data in the usual way with with module mis2 .