Fast Spherical Transforms: SpharmonicKit: Modified FFTPACK

SpharmonicKit: Using the Modified FFTPACK

This file is basically the file HOWTO_FFTPACK found in SpharmonicKit.

A considerably speed-up in Legendre transform runtime is possible by using a slightly modified version of FFTPACK, the freely available collection of Fortran routines developed at NCAR which compute 1-dimensional Fourier (and related) transforms. The original, unmodified version of the library is available here.

The modified version of FFTPACK, the version that SpharmonicKit uses, is available in the following formats:

.tar.Z (around 41 kb)
.tar.gz (around 24 kb)
.zip (around 46 kb)

We make the same statement concerning our modified FFTPACK as NCAR does about the original: the modified FFTPACK is provided free of charge, but there is no guarantee.

Directions to make and install the modified FFTPACK library are included in the distribution of FFTPACK we provide.

The specific FFTPACK routines used by SpharmonicKit are:

rffti ( initialize arrays necessary for doing ffts )
cosqi ( initialize arrays necessary for doing dcts )
rfftf ( forward fft: sample to coefficient )
rfftb ( inverse fft: coefficient to sample )
cosqb ( our forward dct: sample to coefficient )
cosqf ( our inverse dct: coefficient to sample )

The functions rffti, cosqi, rfftf and rfftb were not modified in anyway. The functions cosqb and cosqf were modified. Exact details are given later in this document.

To allow easy use of the FFTPACK-based routines (which are written in Fortran), small C-"interface" functions were written that call the appropriate FFTPACK-library routines. This should make substituting your own optimized FFT and DCT routines relatively straightfoward - just modify the definitions of the "interface" functions. If you use your own FFT and DCT routines, make sure that all scaling of coefficients and samples is done properly.

IMPORTANT: In order to use the FFTPACK-based (or your own) routines, define the symbol FFTPACK in the Makefile before compiling, and make sure the symbol FFTFLAGS in the Makefile is defined and is pointing to where the FFTPACK (or your own) library lives.

The FFT routines

The interface functions for the fft routines are

precomp_fft ( calls the fftpack routine rffti )
grid_fourier_FFTPACK ( calls the fftpack routine rfftf )
grid_invfourier_FFTPACK ( calls the fftpack routine rfftb )

and they are defined in fft_grids.c . If you want to use your own FFT, these are the three functions you want to modify.

The DCT routines

The interface functions for the dct routines are

precomp_dct ( calls the fftpack routine cosqi )
DCTf ( "f" for forward: calls the modified fftpack routine cosqb )
DCTb ( "b" for backward: calls the modified fftpack routine cosqf )

and they are defined in newFCT.c . If you want to use your own dct, these are the three functions you want to modify.

The DCT Details

DCTf and DCTb use slightly modified versions of the fftpack functions cosqb ( in DCTf ) and cosqf (in DCTb ). Before compiling the fftpack library, these functions were modified in order to scale the coefficients the same way as they are in kFCT, kFCTX, and ExpIFCT (i.e. the routines we provide). Originally, the transforms cosqb and cosqf were not orthogonal. THEY ARE NOW.

Exactly how was the scaling modified in the original, unmodified fftpack routines cosqb, cosqf ? Let X be an input array of length N. Let Y be the output of kFCT (or kFCTX - they're identical), Z be the output of the original, unmodified fftpack routine cosqb. Let SCALE = 1 / (2 * N) .Then (using C and not Fortran indexing)

Y[0] = SCALE * Z[0] / 2
Y[i] = SCALE * Z[i] for i = 1, 2, ..., N-1

The original FFTPACK library routine cosqb was further modified to return only the first p <= N many coefficients. This required modifying the original, unmodified fftpack-source file cosqb.f .

Let X be an input array of length N. Let Y be the output of ExpIFCT. Now let me multiply the first coefficient of X by 2 ( so X[0] *= 2 ) and plug this modified X into the original, unmodified fftpack-routine cosqf. Let Z be the resulting output of cosqf. Let SCALE = 1/2. Then (using C and not Fortran indexing)

Y[i] = SCALE * Z[i] for i = 0, 1, 2, ..., N-1

All this scaling is now done within the source files cosqf1.f and cosqb1.f . It's a little hokey, but it works.

Return to SpharmonicKit