(17) |

(18) |

(19) |

We can solve the Helmholtz equation on a regular grid by approximating the differential operator with a finite-difference stencil. A simple (all-zero) convolutional approximation to the Laplacian, , produces the matrix equation:

(20) |

Unfortunately the direct solution of equation () requires the inversion of a multi-dimensional convolutional matrix. As with Poisson's equation above, the application of helical boundary conditions allows me to factor the convolutional filter into a minimum-phase causal and anti-causal pair that can be inverted rapidly by polynomial division.

Factoring the Helmholtz operator is not as simple as factoring the Poisson operator, since its spectrum is not positive definite. The spectrum of the differential Helmholtz operator can be obtained by taking the spatial Fourier transform of equation (), to give

(21) |

Fortunately replacing by , where is a small positive number, successfully stabilizes the spectrum, by pushing the function off the negative-real axis. Equation (), therefore becomes

(22) |

Generalizing the concept of spectral factorization to cross-spectral factorization, Claerbout (1998c) showed that any function whose Fourier transform does not wrap around the origin in the complex plane can be factored into the crosscorrelation of two minimum-phase factors. Since for this class of function, the phase of the Fourier component at the positive Nyquist equals the phase at the negative Nyquist (with no bulk phase shift), he termed this class of function `level-phase'. Although the Helmholtz operator is not strictly an autocorrelation, it does satisfy the `level-phase' criterion, and so it can still be factored into a pair of minimum-phase factors.

Rather than considering a simple convolutional approximation to the Laplacian, we can form a rational approximation,

(23) |

Inserting equation () into equation () yields a matrix equation of similar form, but with increased accuracy at high spatial wavenumbers:

(24) | ||

(25) |

The operator on the left-hand-side of equation ()
represents a multi-dimensional convolution matrix, that can be mapped
to an equivalent one-dimensional convolution by applying helical
boundary conditions.
Although the complex coefficients on the main diagonal cause
the matrix not to be Hermitian, the spectrum of the matrix
is of level-phase.
Therefore, for constant *v*, it can be
factored into causal and anti-causal (triangular) components with any
spectral factorization algorithm that has been adapted for
cross-spectra Claerbout (1998c).

(26) |

The challenge of extrapolation is to find that satisfies both the above equation and our initial conditions, . Starting from , we can invert recursively to obtain a function that satisfies both the initial conditions, and

(27) |

5/27/2001