Getting Started

Documentation

Glish

Learn More

Programming

Contact Us

Version 1.9 Build 803

News	FAQ
Search	Home

FFT and Convolutions

Profiling of imager showed that the application spent maximum fraction of time in the Copy::objcopy() method. This method is used to copy data between memory buffers and all these calls were triggered either in the Table system or in the FFTServer (via the LatticeConvolver::convolve() -> LatticeFFT::fft() -> FFTServer::fft() -> FFTServer::fft0() -> FFTServer::flip() -> Copy::objcopy() chain). We therefore examined FFTServer with the goal of minimizing memory copies if possible.

AIPS++ uses the FFTPACK library for the various versions of the FFT algorithm (Real->Complex, Complex->Complex, Complex->Real). FFTServer is the C++ wrapper class which provides these FFT services. Images are represented using the TempLattice classes, which in turn are derived from the Lattice base class and uses the ArrayLattice if the desired image fits in the available memory or PagedArray if the available memory is not sufficient. The decision about the available memory is controlled by the .aipsrc variable system.resources.memory.

FFT and FFT based convolution operations on images are done via the LatticeFFT and LatticeConvolver classes respectively. As mentioned in Section 2, FFT is used once for making the PSF, and then in every Clark Clean major cycle via the convolution for making the model visibilities. Convolution of the PSF with the CC model is done via the following operations:

1. Flip the PSF, take the FFT to make the transfer function (XFR) and flip the XFR back again

2. Flip the CC image, take the FFT (V), and flip the result again

3. multiply the XFR with V, flip the result, do an inverse FFT and flip the result again

The flipping operation is a circular shift of the functions along each axis by N/2 pixels where N is the size of the image along an axis. Flipping an image via a call to FFTServer::flip() costs O(3N) calls per axis to Copy::objcopy(). The FFT algorithm expects the origin of the function along each axis to be located at the first pixel along each axis - often referred to as the flipped order. Functions (images in our case) represented in the natural order have their original located at N/2. Without the flipping operation, the Fourier transform of a function I is given by

$$\rightleftharpoons \pi \iota \left[\vphantom{u N_x/2 + v N_y/2}\right. \left.\vphantom{u N_x/2 + v N_y/2}\right]$$ (1)

where N_x and N_y are image sizes along the two axis and V will be in the flipped order.

The flipping operation on I eliminates the extra phase term in the Fourier domain. For the purpose of convolution then, it is sufficient to flip the PSF only. Without the flips in Step 2, V will carry the phase term. Complex point-by-point multiplication of XFR with V will ensure that after taking an inverse Fourier transform, the resultant convolution will be in the natural order. This reduces the number of flips to only O(N) per axis of the PSF, and is independent of the number of convolutions (provided, of course, the transfer function does not change!)⁸.