Finite differences and spectral factorization

In the one-dimensional case, a finite-difference representation of the squared Laplacian can be defined as a centered 5-point filter with coefficients (1,-4,6,-4,1). On the same grid, the Laplacian operator can be approximated to the same order of accuracy with the filter (1/12,-4/3,5/2,-4/3,1/12). Combining the two filters in accordance with equation (3) and performing a spectral factorization with one of the standard methods Claerbout (1976, 1992), we can obtain a 3-point minimum-phase filter, suitable for inverse filtering. Figure 1 shows a family of one-dimensional minimum-phase filters for different values of the parameter t. Figure 2 demonstrates the interpolation results obtained with these filters on a simple one-dimensional synthetic. As expected, a small tension value (t=0.01) produces a smooth interpolation, but creates artificial oscillations in the unconstrained regions around sharp changes in the gradient. The value of t=1 leads to linear interpolation with no extraneous inflections, but with discontinuous derivatives. Intermediate values of t allow us to achieve a compromise: a smooth surface with constrained oscillations.

 

otens Figure 1 One-dimensional minimum-phase filters for different values of the tension parameter t. The filters range from the second derivative for t=0 to the first derivative for t=1.

 

int
Figure 2 Interpolating a simple one-dimensional synthetic with recursive filter preconditioning for different values of the tension parameter t. The input data is shown on the top. The interpolation results range from a natural cubic spline interpolation for t=0 to linear interpolation for t=1.

To design the corresponding filters in two dimensions, I define the finite-difference representation of operator (3) on a 5-by-5 stencil. The filters coefficients are chosen with the help of the Taylor expansion to match the desired spectrum of the operator around the zero spatial frequency. The matching conditions lead to the following set of coefficients for the squared Laplacian:

-1/60 2/5 7/30 2/5 -1/60

2/5 -14/15 -44/15 -14/15 2/5

7/30 -44/15 57/5 -44/15 7/30

2/5 -14/15 -44/15 -14/15 2/5

-1/60 2/5 7/30 2/5 -1/60

= 1/60

-1 24 14 24 -1

24 -56 -176 -56 24

14 -176 684 -176 14

24 -56 -176 -56 24

-1 24 14 24 -1

Laplacian representation with the same order of accuracy has the coefficients

-1/360 2/45 2/45 -1/360

2/45 -14/45 -4/5 -14/45 2/45

-4/5 41/10 -4/5

2/45 -14/45 -4/5 -14/45 2/45

-1/360 2/45 2/45 -1/360

= 1/360

-1 16 16 -1

16 -112 -288 -112 16

-288 1476 -288

16 -112 -288 -112 16

-1 16 16 -1

For the sake of simplicity, I assumed an equal physical spacing in x and y directions. The coefficients can be easily adjusted for anisotropic spacing. Figures 3 and 4 show the spectra of the finite-difference representations of operator (3) for the different values of the tension parameter. The finite-different spectra appear as fairly isotropic. They match the exact expressions at small frequencies.

 

specc
Figure 3 Spectra of the finite-difference splines-in-tension schemes for different values of the tension parameter (contour plots).

 

specp
Figure 4 Spectra of the finite-difference splines-in-tension schemes for different values of the tension parameter (cross-section plots). Dashed lines show the exact spectra for continuous operators.

Regarding the finite-difference operators as two-dimensional auto-correlations and applying the efficient Wilson-Burg method of spectral factorization Claerbout (1999); Sava et al. (1998), I obtain two-dimensional minimum-phase filters suitable for inverse filtering. The exact filters contain many coefficients, which rapidly decrease in magnitude at a distance from the first coefficient. For reasons of efficiency, it is advisable to restrict the shape of the filter so that it contains only the valuable coefficients. Keeping all the coefficients that are 1000 times smaller in magnitude than the leading coefficient creates a 53-point filter for t=0 and a 35-point filter for t=1, with intermediate filter lengths for intermediate values of t. When the ratio is changed to 200, we obtain 25- and 16-point filters, respectively. The restricted filters don't factor the autocorrelation exactly, but provide an effective approximation of the exact factors. As outputs of the Wilson-Burg spectral factorization process, they obey the minimum-phase condition.

 

splin
Figure 5 Inverse filtering with the tension filters. The left plots show the inputs composed of filters and spikes. Inverse filtering turns filters into impulses and turns spikes into inverse filter responses (middle plots). Adjoint filtering creates smooth isotropic shapes (right plots). The tension parameter takes values 0, 0.3, 0.7, and 1 (from top to bottom).

Figure 5 shows the two-dimensional filters for different values of t and illustrates inverse recursive filtering, which is the essence of the helix method Claerbout (1998a,b, 1999). The case of t=1 leads to the filter known as helix derivative Claerbout (1999); Zhao (1999). The filter values are spread mostly on two columns. The other boundary case of t=0 leads to a three-column filter, which serves as the minimum-phase version of the Laplacian. As expected from the theory, the inverse impulse response of this filter is noticeably smoother and wider than the inverse response of the helix derivative. Filters corresponding to intermediate values of t exhibit intermediate properties. Theoretically, the inverse impulse response of the filter corresponds to the Green function of equation (3). The theoretical Green function for the case of t=1 is   where r is the distance from the impulse: $r=\sqrt{\left(x-x_k\right)^2 + \left(y-y_k\right)}$. In the case of t=0, the Green function is smoother at the origin:  

$$G = \frac{1}{8\pi}r^2\ln{r}\;.$$ (6)

The theoretical Green function expression for an arbitrary value of t is not known, but we can assume that its smoothness lies between the two boundary conditions.

In the next section, I illustrate an application of helical inverse filtering to a two-dimensional interpolation problem.