The routines for fast Fourier transform (and their origin) are:
This section describes the differences between the way NAG and FFTPACK store arrays of Fourier coefficients. In the following, the Fourier transform of an array of N data values is represented by a sequence of N complex values [A0+i*B0], [A1+i*B1], ..., [A(N-1)+i*B(N-1)].
The differences are basically in the organisation of the Fourier coefficients within the returned array, and also in the normalisation. The normalisation of the FFTPACK values is such that doing a forward transform followed by a backward transform will result in the original array values being multiplied by a factor of N.
Routines to do conversions between FFTPACK and NAG formats have been added to the library.
The relevant NAG routines are C06FAF and C06FBF (the “Hermitian” routines), and the FFTPACK routines are PDA_DRFFTI, PDA_DRFFTF and PDA_DRFFTB. These routines take advantage of the symmetries present in the Fourier transform of a purely real sequence. Only half of the real (A) and imaginary (B) terms need to be calculated and stored because the other halves are just the same. This means that only half the space is required to store the Fourier transform (i.e. N elements rather than 2*N), and it takes roughly half the time to evaluate. The disadvantage is that the resulting Fourier transform array can be rather more difficult to use than if all the real and imaginary parts are stored explicitly. There are routines PDA_DNAG2R and PDA_DR2NAG to do in-situ conversions between NAG and FFTPACK format. Note, each of these routines divides the supplied values by SQRT(N), so successive calls to PDA_DR2NAG and PDA_DNAG2R do not leave the original data unaffected (they are divided by N). This is done to cancel the effect of successive calls of PDA_DRFFTF and PDA_DRFFTB which multiplies the original data by N.
The real and imaginary coefficients produced by PDA_DRFFTF are numerically larger than the corresponding C06FAF coefficients by a factor of SQRT(N), and are ordered differently in the returned arrays. Both routines return A0 (i.e. the DC level in the array) in element 1. PDA_DRFFTF then has corresponding real and imaginary terms in adjacent elements, whereas C06FAF has all the real terms together, followed by all the imaginary terms (in reverse order):
The zeroth imaginary term (B0) always has the value zero and so is not stored in the array. Care has to be taken about the parity of the array size. If it is even, then there is one more real term than there are imaginary terms (excluding A0), i.e. if N = 10, then the coefficients are stored as follows:
If N = 9, then the coefficients are stored as follows:
The relevant NAG routine is C06FCF and the FFTPACK routines are PDA_DCFFTI, PDA_DCFFTF and PDA_DCFFTB. These routines take the Fourier transform of a general complex sequence of N values (i.e. 2*N real values), also returning the Fourier transform in a sequence of N complex values. FFTPACK and NAG differ in that FFTPACK stores the real and imaginary parts of each complex value in adjacent elements of the array, whereas NAG has two separate arrays, one for the real terms and one for the imaginary terms. There is also a difference in the normalisation of the routines in that the real and imaginary Fourier coefficients produced by PDA_DRFFTF are numerically larger than the corresponding C06FAF coefficients by a factor of SQRT(N). There are subroutines PDA_DNAG2C and PDA_DC2NAG to convert between NAG and FFTPACK format. Successive calls to PDA_DC2NAG and PDA_DNAG2C will result in the original data being divided by N. This is done to cancel the multiplication by N which occurs when successive calls to PDA_DCFFTF and PDA_DCFFTB are made.
C06FAF is the NAG routine for finding the FFT of a one-dimensional sequence of real data values. The routine performs a forward transform, storing the FFT as a “Hermitian” sequence in which only half of the real and imaginary terms are kept. The inverse transform is obtained by calling C06FBF, which finds the FFT of a one-dimensional Hermitian sequence.
The following steps are involved in replacing C06FAF calls with equivalent FFTPACK calls:
The work array passed to the FFT routine needs to be increased in size from N elements (for NAG) to 2*N+15 (for FFTPACK).
Replace the call
with
The work array supplied to PDA_DRFFTF needs initialising before calling PDA_DRFFTF. This is done by calling PDA_DRFFTI:
There is no need to re-initialise WORK if the value of N has not changed since the previous call to PDA_DRFFTI (and if the contents of the work array have not been altered). No harm will occur (except for significant slowing down of execution) if the WORK array is unnecessarily re-initialised, but it is a good idea to include some logic to prevent this.
Compared to the Fourier coefficients created by NAG, those created by FFTPACK are stored in a different order in the output array and are normalised differently. You can either modify your application to use the FFTPACK format throughout, or call the PDA_DR2NAG routine to convert the FFTPACK results into NAG format.
where X is the output from PDA_DRFFTF. On return, X holds a NAG-style Hermitian sequence.
C06FBF is the NAG routine for finding the FFT of a one-dimensional Hermitian sequence such as created by C06FAF. The routine performs a forward transform, but it is usually used to perform an inverse transform by preceeding it with a call to C06GBF to form the complex conjugates of the input (frequency domain) data.
The following steps are involved in replacing C06FBF calls with equivalent FFTPACK calls:
Compared to the Hermitian sequences created by NAG, those created by FFTPACK are stored in a different order and are normalised differently. You can either modify your application to use the FFTPACK format throughout, or call the PDA_DNAG2R routine to convert the supplied NAG format data into the equivalent FFTPACK format data:
where X is the supplied NAG-style data. On return, X holds the FFTPACK-style data, ready for use by PDA_DRFFTB. If this call is made, the values returned by PDA_DRFFTB will have the same normalisation as the original data supplied to PDA_DRFFTF.
The work array passed to the FFT routine needs to be increased in size from N elements (for NAG) to 2*N+15 (for FFTPACK).
Replace the two calls:
where X is in NAG format, with
where X is in FFTPACK format.
The work array supplied to PDA_DRFFTB needs initialising before calling PDA_DRFFTB. This is done by calling PDA_DRFFTI:
There is no need to re-initialise WORK if the value of N has not changed since the previous call to PDA_DRFFTI (and if the contents of the work array have not been altered).
C06FCF is the NAG routine for finding the FFT of a one-dimensional sequence of complex data values. The routine performs a forward transform. To do an inverse transform the complex conjugate of the input data is taken before calling C06FCF (using C06GCF), and the complex conjugate of the output data is taken on return from C06FCF.
The steps involved in replacing C06FCF calls with equivalent FFTPACK calls are listed separately for forward and inverse transforms.
The NAG routine expects real and imaginary parts in separate arrays, whereas FFTPACK expects them in the same array, with corresponding real and imaginary values in adjacent elements. If the application can be changed to supply the input data in this format, so well and good. Otherwise you will have to have an extra work array in which to hold the input (and output) data in FFTPACK format. You would convert the supplied input data using code such as:
or
where the X and Y arrays hold the supplied data, C is a work array, and N is the number of data points.
The work array passed to the FFT routine needs to be increased in size from N elements (for NAG) to 4*N+15 (for FFTPACK).
Replace the call
with
The work array supplied to PDA_DCFFTF needs initialising before calling PDA_DCFFTF. This is done by calling PDA_DCFFTI:
There is no need to re-initialise WORK if the value of N has not changed since the previous call to PDA_DCFFTI (and if the contents of the work array have not been altered). No harm will occur (except for significant slowing down of execution) if the WORK array is unnecessarily re-initialised, but it is a good idea to include some logic to prevent this.
The Fourier coefficients created by FFTPACK are stored in a single array and are not normalised, whereas NAG stores them in two arrays and normalises them. You can either modify the way your application to use the FFTPACK format instead of the NAG format, or call the PDA_DC2NAG routine to convert the FFTPACK results into NAG format.
where C is the output from PDA_DCFFTF, and X and Y hold the corresponding real and imaginary coefficients as returned by C06FCF.
If you choose not to modify your application to use FFTPACK data format throughout, you can instead do all the conversions just before (and after) calling the FFTPACK routines. So, if your application supplied frequency domain data in NAG format, first convert it to FFTPACK format using the PDA_DNAG2C routine:
where C is an additional work array used to hold the FFTPACK format data, ready for use by PDA_DCFFTB. X and Y are the supplied frequency domain data in NAG format. If this call to PDA_DC2NAG is made, the values returned by PDA_DCFFTB will have the same normalisation as the original data supplied to PDA_DCFFTF.
The work array passed to the FFT routine needs to be increased in size from N elements (for NAG) to 4*N+15 (for FFTPACK).
Using NAG, the inverse transform is usually done by the three calls:
These three calls should be replaced by the single call:
where C is the array into which the X and Y arrays have been converted using the method of the previous section.
The WORK array passed to PDA_DCFFTB should be initialised using PDA_DCFFTI before calling PDA_DCFFTB. Once the array has been initialised it can be used in multiple calls to PDA_DCFFTF and PDA_DCFFTB so long as they all have the same value for N.
NAG puts the spatial domain results into two arrays (one real, one imaginary), whereas FFTPACK puts them into one. You can either modify your application to use the FFTPACK format or convert the FFTPACK results into NAG-style results using code such as:
or
C06FJF is the NAG routine for finding the FFT of an N-dimensional array of complex data values. The routine performs a forward transformation. To do an inverse transform the complex conjugate of the input data is taken before calling C06FJF (using C06GCF), and the complex conjugate of the output data is taken on return from C06FJF.
There are no equivalent routines in FFTPACK as found in NETLIB. PDA_DNFFTF and PDA_DNFFTB have been written, which do the equivalent of C06FJF. These routines are a bit different to genuine FFTPACK routines in that they do not need any initialisation, and use NAG format rather than native FFTPACK format for complex data arrays and Fourier coefficient arrays. Consequently, replacing C06FJF is a bit easier than replacing the one-dimensional routines.
The steps involved in replacing C06FJF calls with equivalent FFTPACK calls are listed separately for forward and inverse transforms.
The work array passed to the FFT routine needs to be increased in size from 3*MAXDIM elements (for NAG) to 6*MAXDIM+15 (for FFTPACK). Here, MAXDIM is the size of the largest array dimension.
Replace the call
with
The work array passed to the FFT routine needs to be increased in size from 3*MAXDIM elements (for NAG) to 6*MAXDIM+15 (for FFTPACK). Here, MAXDIM is the size of the largest array dimension.
Using NAG, the inverse transform is usually done by the three calls:
These three calls should be replaced by the single call:
C06FUF is the NAG routine for finding the FFT of a two-dimensional sequence of complex data values. There is no direct equivalent. Use the N-dimensional routines instead (with N = 2). See C06FJF.
The complex conjugation NAG routines C06GBF and C06GCF should no longer be needed since separate routines are provided within FFTPACK for doing inverse transformation.