Reducing the filter length

Although mostly zeros, the 9-point filter we 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, we 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 we 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 we still encounter numerical problems depending on the factorization algorithm, and the two factors are barely minimum-phase.

At the time of this report, we have been using the Fourier-domain Kolmogoroff method (, ) to factor the Helmholtz equation. The Kolmogoroff 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, we need to pad the filters to over 4,000 times their original length. Secondly, the Kolmogoroff 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 (, ) may eventually overcome these problems. 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, so even the best result before divergence starts may be unusable.