Fourier methods of seismic data regularisation

Anti-Leakage Fourier Transform

The first documented solution to the leakage problem was performed by Hogbom (1974) as a solution within Astrophysics. The situation was similar to that of seismic data regularisation - radio interferometric baselines after aperture synthesis did not fall on a regular grid as desired, and the adaptive process known as `The CLEAN Algorithm' was devised.

An adaptation of this can give a robust solution to Fourier regularisation, as discussed by Xu and Pham (2004a), Xu and Pham (2004b) and Schonewille (2009), and has been named the `anti-leakage Fourier transform' due to its ability to suppress spectral leakage. This is an iterative algorithm, which exploits the fact that leakage is caused by high energy components contaminating other Fourier components. The maximum component of a given spectrum is calculated and added to an `estimated' spectrum, then the contribution of this component is removed from the input data in data space, by inverse transforming the single component (taking care to reconstruct the same irregular axis) and subtracting. The new spectrum will see this peak removed, and then the next highest energetic component is found, added to the estimated spectrum and removed from the input data, and then iterations continue. The number of iterations is problem dependent, but the maximum required to rebuild a leakage free spectrum will be $N$ , the number of Fourier components, and often far fewer will be required. The final spectrum is free of leakage, however some attention must be paid to Fourier artifacts; in this case cosine tapering is applied to the edges of the data

This algorithm can be summarised as follows

1. Compute the irregular, weighted DFT of the input data
2. Select the strongest Fourier component
3. Add this value to the estimated spectrum
4. Inverse irregular DFT this component (onto the irregular input axis)
5. Subtract this component from the input data
6. Iterate until the desired number of Fourier components have been estimated

It is clear from this set up that this algorithm is very limited by problem size. If the input axis is of length $N$ and $N_k$ components are to be estimated then $O(N x N x N_k)$ operations are required. Weights are only axis dependent so these are calculated and stored outside of the main loop.

Fourier methods of seismic data regularisation