Nfft

NFFT: Algorithm for irregular sampling
Akshay Gulati and Robert J. Ferguson

Summary
The nonuniform discrete Fourier transform (NDFT), used in many processing schemes, can be computed using a fast algorithm known as the non uniform fast Fourier transform (NFFT). The NFFT is not a new algorithm, but it is an approximation scheme that can be use to calculate an approximate spectrum. In onedimension, computational complexity of the NFFT is O(N log N) which is a dramatic improvement from the O(N2) complexity of the NDFT. This algorithm can be easily extended to higher dimensions. The approximate spectrum is calculated using a simple algorithm scheme which involves convolution of an irregularly sampled signal with a truncated Gaussian in the spatial domain. A new empirical expressionbased on numerical experiment for the analytical Gaussian width is proposed. Synthetic data examples, some with analytical solutions, demonstrate the utility and validity of this approach. The approximate spectrum obtained can be use further in a reconstruction algorithm. This algorithm removes the bottle neck from forward process by replacing NDFT with NFFT in many conventional processing algorithms.Introduction
The problem of analyzing a signal P (tj) having irregularly spaced measurements is common in geophysics (Ferguson, 2006). Faulty equipment, errors in positioning, obstacles, and noise sources can be reasons for irregularly spaced measurements. Several approaches can be found which are based on some kind of interpolation techniques, but most of these approaches do not handle thedata optimally. Transformation methods which include Radon transformations have been used to handle the problem of irregular sampling (Clarebout,1992).In the parabolic Radon transform, two CMP gathers are combined to improve offset sampling, and thus differences between mid point positions are ignored. Similarly, hyperbolic and linear Radon transforms (Clarebout,1992) as well as the parabolic Radontransform are suitable for estimating frequencies at irregular nodes, but they suffer aliasing problems due to sparse sampling. Prediction error filtering is another method which is used to interpolate missing traces, but it works only on regularly sampled data, and it only corrects for aliasing. if sampling is irregular, the result will be erroneous. Ronen et al (1991) suggests a method for DMOstacks which can handle irregular sampling, but it is still not efficient for very large seismic data volumes. Approximate regularization/datuming (Ferguson, 2006) allows extrapolation of data recorded on an irregular grid onto a regular grid, but it requires a velocity model. My aim in this paper is to solve the forward case with a faster algorithm, as done by Duijndam and Schonewille (1999) andGreengard and Lee (2004). I will begin with the theory behind the algorithm structure for NFFT. After reviewing the basic algorithms, improvements are suggested in the NFFT algorithm. Then, using synthetic examples to validate the algorithm and I show that the NFFT is approximately 100 times faster than the NDFT, and that the approximate spectrum obtained using NFFT gives a good approximation tothe original spectrum.

Theory
The NFFT algorithm can be numerically expressed in following steps: convolution, FFT, and deconvolution. Two parameters are significant in our algorithm, one is numerical width q of the truncated filter, and other is analytical width b for Gaussian filter. Convolution with the short Gaussian filter g(x) is carried out to make the signal approximately band-limitedaccording to

GeoCanada 2010 – Working with the Earth

1

p g ( m ) = g ( x ) * p( x )

(1)

where pg(m) is the result is the result of spatial convolution. Equation 1 can be written as multiplication in the Fourier domain as

Pg ( m ) = G( m ) * P ( m )

(2)

where pg(m) is the Fourier spectrum of pg(m) in Fourier domain. For efficiency Gaussian need to be truncated, thus...