REFERENCES

Claerbout, J., 1976, Fundamentals of Geophysical Data Processing: http://sepwww.stanford.edu/sep/prof/.

Claerbout, J., 1997, Multidimensional recursive filters via a helix: SEP-95, 1-13.

Kolmogoroff, A. N., 1939, Sur l'interpolation et l'extrapolation des suites stationnairres: C.R. Acad.Sci., 208, 2043-2045.

Rickett, J., and Claerbout, J., 1998, Helical factorization of the Helmholtz equation: SEP-97, 353-362.

Wilson, G., 1969, Factorization of the covariance generating function of a pure moving average process: SIAM J. Numer. Anal., 6, no. 1, 1-7.

A We use the code below to factor cross-correlation functions into two minimum phase filters. The module contains three subroutines, one for initialization, one for the actual factorization and a third one for deallocating the memory.

Type cfilter is a complex filter defined in the helical coordinate system. Such a filter is composed of a ``zero lag'' coefficient (e.g. s0) and two vectors, one for the lags (e.g. ss%lag) and one for the complex values of the filter (e.g. ss%flt).

We loop for the desired number of iterations. At each step, we take the auto-correlation (ss), divide it by A and $\bar B$, select the filter coefficients of the positive/negative lags, and then convolve by A/B.

The result is a pair of filters (aa and bb) of type cfilter.

module ccwil { 
# Factorization of complex cross-spectra using the Wilson-Burg algorithm
  use chelicon
  use cpolydiv
  use cprint
#----------------------------------
  integer,                             private :: n         
  complex, dimension( :), allocatable, private :: cross, ab, cc, b, c  
contains
  subroutine ccwil_init( nmax) {
    integer, intent(in) :: nmax; n = nmax
    allocate( cross(2*n-1), ab(2*n-1), cc(2*n-1), b(n), c(n))
    }
  subroutine ccwil_factor( niter, s0, ss, aa, bb) { 
    integer,                intent( in)     :: niter
    complex,                intent( in)     :: s0
    type( cfilter),         intent( in)     :: ss
    type( cfilter),         intent(out)     :: aa, bb    
    integer                                 :: i,stat
    real                                    :: eps
    logical                                 :: T=.true.,F=.false.
    cross = 0.; cross( n) = s0; cross( n+ss%lag) = ss%flt
    ab=0.0; cc=0.0; b=0.0; b(1)=1.0; c=0.0

do i = 1, niter {
       call cpolydiv_init( aa)
       stat = cpolydiv_clop( F, F, cross, ab)  # ab = S/A
       call cpolydiv_init( bb)
       stat = cpolydiv_clop( T, F, cc   , ab)  # cc = S/(AB*')

b( 2:n) = cc( n+1:2*n-1)/cc(n)
       eps = maxval(abs(b(2:n)))
       call chelicon_init( aa)
       stat = chelicon_clop( F, F,     b,  c)  # c = A b
       aa%flt = c(1+aa%lag)

b( 2:n) = cc( n-1:1:-1 )/cc(n)
       eps = max(abs(maxval(abs(b(2:n)))),abs(eps))
       call chelicon_init( bb)
       stat = chelicon_clop( F, F,     b,  c)  # c = B* b
       bb%flt = c(1+bb%lag)
       write( 0,*) i, eps
       }
       bb%flt=conjg(bb%flt)
  }
  subroutine ccwil_close () {
    deallocate( cross, ab, cc, b, c) }
}