## 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 , 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) }
}