Introduction

Many algorithms in seismic data processing require that input data have uniform sampling intervals both in time and in space. However, in reality, the difficulties in positioning the sources and receivers may cause the spatial sampling intervals to vary from place to place. If the variations of the sampling intervals are too large to be acceptable, an interpolation process must be done. Because seismic data are almost always sampled uniformly in time, the spatial interpolation can be done along each constant time slice or each constant frequency slice. Therefore, the problem can be reduced to one dimension.

Although the interpolation of seismic data beyond aliasing is still under development, the theories and techniques for interpolating unaliased data are well understood and routinely used in industry. In the 1950s, Yen (1956 and 1957) developed an algorithm for interpolating or resampling irregularly sampled, one-dimensional data. In his algorithm, the irregularly sampled data are modeled as the interpolation of unknown uniformly sampled data:  

$$s(t_i) = \sum_j s(j\Delta t)\hbox{sinc}(t_i-j\Delta t) \ \ \ \ i=1,\ldots,N.$$ (1)

Yen solved this linear system for $s(j\Delta t)$, and derived an exact, explicit formula for evaluating the continuous signal s(t).

Hale (1980) tested Yen's algorithm and pointed out that, although the method is exact for noise-free data, it could give erroneous results for noisy data. With Yen's method, a single noisy sample, even one far away from the correct sampling positions, may have an undesirable influence on the values of the interpolated samples. Moreover, Yen's algorithm is computationally expensive. Resampling N points requires an order of N³ operations. To solve these two problems, Hale replaced the sinc function in equation (1) with a truncated, tapered sinc function. Although this procedure sacrifices accuracy at high frequencies, it performs better and at much less computational cost than Yen's method for data contaminated with truncations, aliasing data or noisy data.

This paper describes a new, one-dimensional interpolation method. In this method, a continuous signal is modeled with its Fourier series expansion. The coefficients of this expansion are discrete spectra of the signal, and can be estimated by fitting the modeled data to the irregularly sampled data using the standard least squares technique. It turns out that the normal matrix in this least squares problem is a Hermitian Toeplitz matrix, and its inverse requires only an order of N² operations. Therefore, the algorithm is quite efficient. Assuming that the irregularly sampled signal is not aliased and that it has a finite duration, the continuous signal can be recovered from the irregular samples exactly. This algorithm can also be used to compute the discrete Fourier transform of the irregularly sampled data.

In the following sections, I first formulate the problem and observe the Toeplitz structure of the matrix to be inverted. Then, I show two examples of one-dimensional interpolation and one example of computing the discrete Fourier transform. Finally, I offer conclusions and discuss some potential problems of the method.