Reducing the filter length

Although mostly zeros, the 9-point filter I factorize in 2-D is actually 2 N_x+3 points long on a helix, and so its two factors each contain N_x+2 points. For Figures ,  and  above, I did not discard any of these filter coefficients: our recursive filters contained the full 128 coefficients. For this extrapolator to compete effectively with other methods (especially in 3-D problems), the number of filter coefficients has to be reduced significantly.

The factor coefficients themselves should be independent of the diameter of the helix. This leads us to expect many zero value coefficients on the helix `backside'. Fortunately this is observed in practice, and the filter coefficient amplitude drops rapidly from either end.

Unfortunately, however, the operator that I factor has roots very close to the unit circle in the complex plane. The Laplacian already has a pair of roots at Z=0, and the effect of the extra term on the main diagonal is to destabilize the operator further. The damping rescues us in theory, but in practice I still encounter numerical problems depending on the factorization algorithm.

The results in this thesis have been generated with the Fourier-domain Kolmogorov method Claerbout (1998c); Kolmogorov (1939) to factor the Helmholtz equation. The Kolmogorov method has two main problems, both due to the proximity of the roots to the unit circle. Firstly, circular boundary conditions require us to pad the cross-correlation function before transforming it to the Fourier domain. With the roots close to the unit circle, extreme amounts of padding are needed: in the 2-D examples above, I need to pad the filters to over 4,000 times their original length. Secondly, the Kolmogorov method simultaneously computes all the filter coefficients. With roots so close to the unit circle, truncating filter coefficients, even in a reasonable manner, often leads to non-minimum-phase filters and divergent results.

The Wilson-Burg algorithm Sava et al. (1998); Wilson (1969) does not suffer from the same difficulties. By working in the time domain, the algorithm avoids circular boundary conditions, and the number of filter coefficients can be defined at each iteration, providing a best-fit filter with a given number of coefficients. However, the Wilson-Burg algorithm also encounters problems with roots close to the unit circle. Specifically, numerical problems cause filters to lose their minimum-phase nature, causing the algorithm to diverge. With roots very close to the unit circle, this can happen within the first couple of iterations.