### 8 Fast Fourier transform (FFT)

The routines for fast Fourier transform (and their origin) are:

#### 8.1 Differences between NAG and FFTPACK

• FFTPACK expects and returns data in a different format to NAG.
• FFTPACK includes initialisation routines (PDA_RFFTI and PDA_CFFTI) which should be called prior to the other routines, but which don’t need to be called again until the size of the data array changes. There are no equivalent NAG initialisation routines. The main NAG FFT routines (e.g. C06FAF, etc.) do this initialisation each time they are called, irrespective of the array size.
• FFTPACK has separate forward and backward transform routines, whereas NAG only has forward routines (backward transforms are performed by using complex conjugation with the forward transform). This means that there is probably no need to supply equivalents to the complex conjugation NAG routines (a trivial operation anyway).
• FFTPACK can accept arrays of any length, whereas NAG puts some restrictions on the array length (no prime factor larger than 19 allowed in the array size, and the total number of prime factors must be less than 21).
• FFTPACK routines require a differently sized work space array.
• FFTPACK has no error checking.

#### 8.2 Data formats for FFTPACK and NAG

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.

##### 8.2.1 Fourier transforms of sequences of purely real values

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):

PDA_DRFFTF:  A0,    A1, B1,     A2, B2,     A3, B3,   ...
C06FAF:      A0,    A1, A2, A3, ...,        ..., B3, B2, B1

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:

PDA_DRFFTF:  A0, A1, B1, A2, B2, A3, B3, A4, B4, A5
C06FAF:      A0, A1, A2, A3, A4, A5, B4, B3, B2, B1

If N = 9, then the coefficients are stored as follows:

PDA_DRFFTF:  A0, A1, B1, A2, B2, A3, B3, A4, B4
C06FAF:      A0, A1, A2, A3, A4, B4, B3, B2, B1
##### 8.2.2 Fourier transforms of sequences of complex values

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.

#### 8.3 Replacing calls to C06FAF

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:

• Increase the size of the work array

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 call to C06FAF with PDA_DRFFTF

Replace the call

DOUBLE PRECISION X(N), WORK(N)
CALL C06FAF( X, N, WORK, IFAIL )

with

DOUBLE PRECISION X(N), WORK(2*N+15)
CALL PDA_DRFFTF( N, X, WORK )
• Add calls to PDA_DRFFTI if necessary

The work array supplied to PDA_DRFFTF needs initialising before calling PDA_DRFFTF. This is done by calling PDA_DRFFTI:

DOUBLE PRECISION WORK(2*N+15)
CALL PDA_DRFFTI( N, WORK )

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.

• Convert output (frequency domain) data to NAG format

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.

DOUBLE PRECISION X(N)
CALL PDA_DR2NAG( N, X )

where X is the output from PDA_DRFFTF. On return, X holds a NAG-style Hermitian sequence.

#### 8.4 Replacing calls to C06FBF

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:

• Convert input (frequency domain) data to FFTPACK format

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:

DOUBLE PRECISION X(N)
CALL PDA_DNAG2R( N, X )

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.

• Increase the size of the work array

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 call to C06FBF and C06GBF with PDA_DRFFTB

Replace the two calls:

DOUBLE PRECISION X(N), WORK(N)
CALL C06GBF( X, N, IFAIL )
CALL C06FBF( X, N, WORK, IFAIL )

where X is in NAG format, with

DOUBLE PRECISION X(N), WORK(2*N+15)
CALL PDA_DRFFTB( N, X, WORK )

where X is in FFTPACK format.

• Add calls to PDA_DRFFTI if necessary

The work array supplied to PDA_DRFFTB needs initialising before calling PDA_DRFFTB. This is done by calling PDA_DRFFTI:

DOUBLE PRECISION WORK(2*N+15)
CALL PDA_DRFFTI( N, WORK )

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).

#### 8.5 Replacing calls to C06FCF

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.

##### 8.5.1 Forward transforms
• Re-organise the input (spatial domain) data

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:

DOUBLE PRECISION X( N ), Y( N ), C( 2*N )
DO J = 1, N
I = 2*J
C( I - 1 ) = X( J )
C( I ) = Y( J )
END DO

or

DOUBLE PRECISION X( N ), Y( N ), C( 2, N )
DO J = 1, N
C( 1, J ) = X( J )
C( 2, J ) = Y( J )
END DO

where the X and Y arrays hold the supplied data, C is a work array, and N is the number of data points.

• Increase the size of the work array

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 call to C06FCF with PDA_DCFFTF

Replace the call

DOUBLE PRECISION X(N), Y(N), WORK(N)
CALL C06FCF( X, Y, N, WORK, IFAIL )

with

DOUBLE PRECISION C(2*N), WORK(4*N+15)
CALL PDA_DCFFTF( N, C, WORK )
• Add calls to PDA_DCFFTI if necessary

The work array supplied to PDA_DCFFTF needs initialising before calling PDA_DCFFTF. This is done by calling PDA_DCFFTI:

DOUBLE PRECISION WORK( 4*N+15 )
CALL PDA_DCFFTI( N, WORK )

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.

• Convert output (frequency domain) data to NAG format

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.

DOUBLE PRECISION X(N), Y(N), C(2*N)
CALL PDA_DC2NAG( N, C, X, Y )

where C is the output from PDA_DCFFTF, and X and Y hold the corresponding real and imaginary coefficients as returned by C06FCF.

##### 8.5.2 Inverse transforms
• Convert input (frequency domain) data to FFTPACK format

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:

DOUBLE PRECISION X(N), Y(N), C(2*N)
CALL PDA_DNAG2C( N, X, Y, C )

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.

• Increase the size of the work array

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 call to C06FCF and C06GCF with PDA_DCFFTB

Using NAG, the inverse transform is usually done by the three calls:

CALL C06GCF( Y, N, IFAIL )
CALL C06FCF( X, Y, N, WORK, IFAIL )
CALL C06GCF( Y, N, IFAIL )

These three calls should be replaced by the single call:

CALL PDA_DCFFTB( N, C, WORK )

where C is the array into which the X and Y arrays have been converted using the method of the previous section.

• Add calls to PDA_DCFFTI if necessary

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.

• Re-organise the output (spatial domain) data

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:

DOUBLE PRECISION X( N ), Y( N ), C( 2*N )
DO J = 1, N
I = 2*J
X( J ) = C( I - 1 )
Y( J ) = C( I )
END DO

or

DOUBLE PRECISION X( N ), Y( N ), C( 2, N )
DO J = 1, N
X( J ) = C( 1, J )
Y( J ) = C( 2, J )
END DO

#### 8.6 Replacing calls to C06FJF

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.

##### 8.6.1 Forward transforms
• Increase the size of the work array

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 call to C06FJF with PDA_DNFFTF

Replace the call

DOUBLE PRECISION X(N), Y(N), WORK( 3*MAXDIM )
INTEGER ND( NDIM )
CALL C06FJF( NDIM, ND, N, X, Y, WORK, LWORK, IFAIL )

with

DOUBLE PRECISION X(N), Y(N), WORK( 6*MAXDIM + 15 )
INTEGER ND( NDIM )
CALL PDA_DNFFTF( NDIM, ND, X ,Y, WORK, ISTAT )
IF( ISTAT .NE. 0 ) THEN
This means that NDIM was either less than 1 or greater than 20.
Report a programming error!
END IF
##### 8.6.2 Inverse transforms
• Increase the size of the work array

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 call to C06FJF and C06GCF with PDA_DNFFTB

Using NAG, the inverse transform is usually done by the three calls:

CALL C06GCF( Y, N, IFAIL )
CALL C06FJF( NDIM, ND, N, X, Y, WORK, LWORK, IFAIL )
CALL C06GCF( Y, N, IFAIL )

These three calls should be replaced by the single call:

CALL PDA_DNFFTB( NDIM, ND, X ,Y, WORK, ISTAT )

#### 8.7 Replacing calls to C06FUF

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.

#### 8.8 Replacing calls to C06GBF and C06GCF

The complex conjugation NAG routines C06GBF and C06GCF should no longer be needed since separate routines are provided within FFTPACK for doing inverse transformation.