FFT3F

Computes Fourier coefficients of a complex periodic three-dimensional array.

Required Arguments

A — Three-dimensional complex matrix containing the data to be transformed. (Input)

B — Three-dimensional complex matrix containing the Fourier coefficients of A. (Output)
The matrices A and B may be the same.

Optional Arguments

N1 — Limit on the first subscript of matrices A and B. (Input)
Default: N1 = size(A, 1)

N2 — Limit on the second subscript of matrices A and B. (Input)
Default: N2 = size(A, 2)

N3 — Limit on the third subscript of matrices A and B. (Input)
Default: N3 = size(A, 3)

LDA — Leading dimension of A exactly as specified in the dimension statement of the calling program. (Input)
Default: LDA = size (A,1).

MDA — Middle dimension of A exactly as specified in the dimension statement of the calling program. (Input)
Default: MDA = size (A,2).

LDB — Leading dimension of B exactly as specified in the dimension statement of the calling program. (Input)
Default: LDB = size (B,1).

MDB — Middle dimension of B exactly as specified in the dimension statement of the calling program. (Input)
Default: MDB = size (B,2).

FORTRAN 90 Interface

Generic: CALL FFT3F (A, B [,…])

Specific: The specific interface names are S_FFT3F and D_FFT3F.

FORTRAN 77 Interface

Single: CALL FFT3F (N1, N2, N3, A, LDA, MDA, B, LDB, MDB)

Double: The double precision name is DFFT3F.

Description

The routine FFT3F computes the forward discrete complex Fourier transform of a complex three-dimensional array of size (N1 = N) ´ (N2 = M) ´ (N3 = L). The method used is a variant of the Cooley-Tukey algorithm , which is most efficient when N, M, and L are each products of small prime factors. If N, M, and L satisfy this condition, then the computational effort is proportional to N M L log N M L. This considerable savings has historically led people to refer to this algorithm as the “fast Fourier transform” or FFT.

Specifically, given an N ´ M ´ L array a, FFT3F returns in c = COEF

Furthermore, a vector of Euclidean norm S is mapped into a vector of norm

Finally, note that an unnormalized inverse is implemented in FFT3B. The routine FFT3F is based on the complex FFT in FFTPACK. The package FFTPACK was developed by Paul Swarztrauber at the National Center for Atmospheric Research.

Comments

1. Workspace may be explicitly provided, if desired, by use of F2T3F/DF2T3F. The reference is:

CALL F2T3F (N1, N2, N3, A, LDA, MDA, B, LDB, MDB, WFF1, WFF2, WFF3, CPY)

The additional arguments are as follows:

WFF1 — Real array of length 4 * N1 + 15 initialized by FFTCI. The initialization depends on N1. (Input)

WFF2 — Real array of length 4 * N2 + 15 initialized by FFTCI. The initialization depends on N2. (Input)

WFF3 — Real array of length 4 * N3 + 15 initialized by FFTCI. The initialization depends on N3. (Input)

CPY — Real array of size 2 * MAX(N1, N2, N3). (Workspace)

2. The routine FFT3F is most efficient when N1, N2, and N3 are the product of small primes.

3. If FFT3F/FFT3B is used repeatedly with the same values for N1, N2 and N3, then use FFTCI to fill WFF1(N = N1), WFF2(N = N2), and WFF3(N = N3). Follow this with repeated calls to F2T3F/F2T3B. This is more efficient than repeated calls to FFT3F/FFT3B.

Example

In this example, we compute the Fourier transform of the pure frequency input for a 2 ´ 3 ´ 4 array

for 1 ≤ n ≤ 2, 1 ≤ m ≤ 3, and 1 ≤ l ≤ 4 using the IMSL routine FFT3F. The result

has all zeros except in the (2, 3, 3) position.

USE FFT3F_INT

USE UMACH_INT

USE CONST_INT

IMPLICIT NONE

INTEGER LDA, LDB, MDA, MDB, NDA, NDB

PARAMETER (LDA=2, LDB=2, MDA=3, MDB=3, NDA=4, NDB=4)

! SPECIFICATIONS FOR LOCAL VARIABLES

INTEGER I, J, K, L, M, N, N1, N2, N3, NOUT

REAL PI

COMPLEX A(LDA,MDA,NDA), B(LDB,MDB,NDB), C, H

! SPECIFICATIONS FOR INTRINSICS

INTRINSIC CEXP, CMPLX

COMPLEX CEXP, CMPLX

! SPECIFICATIONS FOR SUBROUTINES

! SPECIFICATIONS FOR FUNCTIONS

! Get output unit number

CALL UMACH (2, NOUT)

PI = CONST('PI')

C = CMPLX(0.0,2.0*PI)

! Set array A

DO 30 N=1, 2

DO 20 M=1, 3

DO 10 L=1, 4

H = C*(N-1)*1/2 + C*(M-1)*2/3 + C*(L-1)*2/4

A(N,M,L) = CEXP(H)

10 CONTINUE

20 CONTINUE

30 CONTINUE

CALL FFT3F (A, B)

WRITE (NOUT,99996)

DO 50 I=1, 2

WRITE (NOUT,99998) I

DO 40 J=1, 3

WRITE (NOUT,99999) (A(I,J,K),K=1,4)

40 CONTINUE

50 CONTINUE

WRITE (NOUT,99997)

DO 70 I=1, 2

WRITE (NOUT,99998) I

DO 60 J=1, 3

WRITE (NOUT,99999) (B(I,J,K),K=1,4)

60 CONTINUE

70 CONTINUE

99996 FORMAT (13X, 'The input for FFT3F is')

99997 FORMAT (/, 13X, 'The results from FFT3F are')

99998 FORMAT (/, ' Face no. ', I1)

99999 FORMAT (1X, 4('(',F6.2,',',F6.2,')',3X))

END

Output

            The input for FFT3F is

Face no. 1
( 1.00, 0.00)   ( -1.00, 0.00)   ( 1.00, 0.00)   ( -1.00, 0.00)
( -0.50, -0.87)   ( 0.50, 0.87)   ( -0.50, -0.87)   ( 0.50, 0.87)
( -0.50, 0.87)   ( 0.50, -0.87)   ( -0.50, 0.87)   ( 0.50, -0.87)

Face no. 2
( -1.00, 0.00)   ( 1.00, 0.00)   ( -1.00, 0.00)   ( 1.00, 0.00)
( 0.50, 0.87)   ( -0.50, -0.87)   ( 0.50, 0.87)   ( -0.50, -0.87)
( 0.50, -0.87)   ( -0.50, 0.87)   ( 0.50, -0.87)   ( -0.50, 0.87)

The results from FFT3F are

Face no. 1
( 0.00, 0.00)   ( 0.00, 0.00)   ( 0.00, 0.00)   ( 0.00, 0.00)
( 0.00, 0.00)   ( 0.00, 0.00)   ( 0.00, 0.00)   ( 0.00, 0.00)
( 0.00, 0.00)   ( 0.00, 0.00)   ( 0.00, 0.00)   ( 0.00, 0.00)

Face no. 2
( 0.00, 0.00)   ( 0.00, 0.00)   ( 0.00, 0.00)   ( 0.00, 0.00)
( 0.00, 0.00)   ( 0.00, 0.00)   ( 0.00, 0.00)   ( 0.00, 0.00)
( 0.00, 0.00)   ( 0.00, 0.00)   ( 24.00, 0.00)   ( 0.00, 0.00)

Visual Numerics, Inc.
Visual Numerics - Developers of IMSL and PV-WAVE
http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260