next up previous index
Next: D703 Real Fast Up: CERNLIB Previous: D601 Solution of

D700 Real Fast Fourier Transform

Routine ID: D700
Author(s): C. IselinLibrary: MATHLIB
Submitter: Submitted: 04.09.1972
Language: FortranRevised: 15.01.1977

Let the discrete Fourier transform be defined by

yj = {1N}∑k=0N-1exp({2πijkN}) xk, (j=0,1,...,N).

The subroutines of package RFT compute this transform or its inverse

xk = {1N}∑j=0N-1exp({-2πijkN}) yj, (k=0,1,...,N)

for real functions, with the restriction that N is a power of 2.

Structure:

SUBROUTINE subprograms
User Entry Names: RFT, RCA, RPA, RPS, RSA
Internal Entry Names: D700SU
Files Referenced: Printer
COMMMON Block Names and Lengths: /D700DT/ 6, /FWORK/ 321

Usage:

CALL RFT(M,X,IX,Y,IY,MODE) or
CALL RCA(M,X,IX,Y,IY) or
CALL RPA(M,X,IX,Y,IY) or
CALL RPS(M,X,IX,Y,IY) or
CALL RSA(M,X,IX,Y,IY)
M
( INTEGER) Number m (such that n=2m ) of input values (full period or half period).
X
( REAL) Input array. The input values are taken from X(k*IX+1) for k=0,1,...,n .
Y
( REAL) Output array. The results are stored in Y(k*IY+1) for j=0,1,...,n .
MODE
( INTEGER) Selects the mode of operation for RFT as follows:
MODE = 1: Analysis of a general real function.
CALL RFT(M,X,IX,Y,IY,1) or
CALL RPA(M,X,IX,Y,IY)
assumes xk=X(k*IX+1) (k=0,1,...,n-1); n=2m=N to define a full period of the function to be analysed. The value xn is ignored. The results are returned in the following order:
y0=yn=Y(1)

yj=yn-j=Y(j*IY+1)+iY((j+n/2)*IY+1), &quad;(j=1,2,...,n/2) .
The other values in Y are not changed.

MODE = 4: Synthesis of a general real function.

CALL RFT(M,X,IX,Y,IY,4) or
CALL RPS(M,X,IX,Y,IY)
is exactly the inverse of MODE=1 as described above. The value xn is set equal to x0 .
MODE=2/5: Analysis/Synthesis of a real even function.
For an even function, the transform is identical to its inverse.
CALL RFT(M,X,IX,Y,IY,2) or
CALL RFT(M,X,IX,Y,IY,5) or
CALL RCA(M,X,IX,Y,IY)
all assume that xk=X(k*IX+1), (k=0,1,...,n), n=2m=N/2 define a half-period of the function to be analysed and that the other half period is generated by even continuation. The results returned are the cosine terms
yj=y2n-j=Y(j*IY+1), &quad;(j=0,1,...,n).
Note that the full period has 2n=N points.
MODE = 3/6: Analysis/Synthesis of a real odd function.
For an odd function the transform is also identical to its inverse. All assume that xk=X(k*IX+1), (k=1,2,...,n) ;
CALL RFT(M,X,IX,Y,IY,3) or
CALL RFT(M,X,IX,Y,IY,6) or
CALL RSA(M,X,IX,Y,IY)
n=2m=N/2 define a half-period of the function to be analysed and that the other half period is generated by odd continuation. The results returned are the sine terms
yj=-y2n-j=Y(j*IY+1), &quad;(j=1,2,...,n) .
Note that y0=yn=0 and that the values returned are Y(1)=X(1) and Y(n*IY+1)=X(n*IX+1) . Again the full period has 2n=N points.

Restrictions:

These subroutines work for any input such that the full period has at least four points, i.e., m≥2 for general functions, or m≥1 for odd or even functions. If the number of data points exceeds 129 (m≤7 ), the calling program must provide sufficient working storage by using the statement

    COMMON /FWORK/ W(nnn)
where nnn=5*2m .

References:

  1. C. Iselin, An approach to fast Fourier transform, CERN 71-19.
A copy of Ref. 1 is available.

D703



next up previous index
Next: D703 Real Fast Up: CERNLIB Previous: D601 Solution of


Janne Saarela
Mon Apr 3 15:06:23 METDST 1995