D702: Complex Fast Fourier Transform

Author(s): R.C. Singeton (Stanford)	Library: MATHLIB
Submitter: B. Fornberg	Submitted: 03.05.1968
Language: Fortran	Revised: 01.10.1974

A discrete Fourier transform is defined by:

and the inverse

satisfying . CFT evaluates these sums using fast Fourier technique. It is not required that N is a power of 2. One-, two- and three-dimensional transforms can be performed.

Structure:

SUBROUTINE subprogram
User Entry Names: CFT
Files Referenced: Printer

Usage:

    CALL CFT(A,B,NTOT,N,NSPAN,ISN)

Arrays A and B originally hold the real and imaginary components of the data, and return the real and imaginary components of the resulting Fourier coefficients.

Multivariate data is indexed according to the Fortran array element successor function, without limit on the number of implied multiple subscripts. The SUBROUTINE is called once for each variate. The calls for a multivariate transform may be in any order. NTOT is the total number of complex data values. N is the dimension of the current variable. NSPAN/N is the spacing of consecutive data values while indexing the current variable. The sign of ISN determines the sign of the complex exponential, and the magnitude of ISN is normally one.

For a single-variate transform, (number of complex data values), e.g.

    CALL CFT(A,B,N,N,N,1)

A tri-variate transform with A(N1,N2,N3), B(N1,N2,N3) is computed by

CALL CFT(A,B,N1*N2*N3,N1,N1,1)

CALL CFT(A,B,N1*N2*N3,N2,N1*N2,1) and

CALL CFT(A,B,N1*N2*N3,N3,N1*N2*N3,1)

The data may alternatively be stored in a single COMPLEX array A, then the magnitude of ISN changed to two to give the correct indexing increment and the second parameter used to pass the initial address for the sequence of imaginary values, for example:

    REAL
    EQUIVALENCE (A,S)
    ...
    CALL CFT (A,S(2),NTOT,N,NSPAN,2)

Arrays AT(MAXF), CK(MAXF), BT(MAXF), SK(MAXF), and NP(MAXP) are used for temporary storage. If the available storage is insufficient the program is terminated by a STOP.

MAXF must be the maximum prime factor of N. MAXP must be > the number of prime factors of N. In addition, if the square-free portion K of N has two or more prime factors, then MAXP must be . Storage in NFAC allows for a maximum of 11 factors of N. If N has more than one square-free factor, the product of the square-free factors must be 210.

Notes:

CFT is very general since the number of points is not restricted to powers of two, as is the case for RFT (D700) and FFTRC (D701). For the routines in FFTRC (D701) are considerably faster.

References:

1. R.C. Singleton, An Algorithm for Computing the Mixed Radix F.F.T., IEEE Trans. Audio Electroacoust., AU-1(1969) 93-107.
2. Reprinted in: L.R. Rabiner and C.M. Rader: Digital Signal Processing, IEEE Press New York (1972) 294.