Page: (fft)Top, Next:
(matrix), Prev:
(complex), Up:
(cpl)Library Index
Fast Fourier Transform
The fft.cpl library handles fast Fourier transforms. All of the following routines accept one or two parameters, being implied that when only one array is specified this acts as both input and output.
SUBROUTINE IFT(COMPLEX Rin^(∗), Rout^(∗)) computes
Rout(i) = SUM Rin(j)∗EXP(2∗PI∗I∗i∗j/Rin.LENGTH) FOR ALL j
SUBROUTINE FFT(COMPLEX Rin^(∗), Rout^(∗)) computes
Rout(j) = [SUM Rin(i)∗EXP(-2∗PI∗I∗i∗j/Rin.LENGTH) FOR ALL i]/Rin.LENGTH
These sums are actually computed by a fast-Fourier-transform algorithm which only works when Rin.LENGTH is the product of a power of two and possibly a single factor of three. The normalizing 1/Rin.LENGTH factor sits in the direct transform, whereas the inverse transform is a simple Fourier sum (the commonest convention for a Fourier series). When the discrete transform is used to approximate a continuous Fourier transform, the normalizing factor will have to be corrected explicitly. The
BOOLEAN FUNCTION FFTfit(INTEGER VARIABLE N)
can be used to test whether N is a suitable dimension for an array to be passed to the FFT routines.
The following variants produce discrete Fourier transforms of real arrays:
SUBROUTINE RFT(COMPLEX Rin^(∗), Rout^(∗))
SUBROUTINE HFT(COMPLEX Rin^(∗), Rout^(∗))
Of these RFT is defined to give the real part of the corresponding IFT, but with double the resolution. HFT is the inverse of RFT.
The output of RFT and the input of HFT is formally a COMPLEX ARRAY (so that either of these transforms may be executed in place) that actually contains a double number of REAL elements. These may be extracted by the ad hoc function
REALIFIED(POINTER TO ARRAY(∗) OF COMPLEX x)
For instance, the REAL ARRAY hidden in the COMPLEX ARRAY Rout may be accessed as REALIFIED(Rout), or equivalently as Rout.REALIFIED.
An additional complication arises from the fact that the HFT of a general array of 2N real numbers (encoded as a COMPLEX ARRAY of size 0..N-1), in addition to having an always real 0th Fourier component, contains an extra, also real, component of frequency N (or equivalently -N) which does not fit in the original array of size 0..N-1. For this reason RFT's and HFT's input/output array can either be of (even) size 0..N-1 or of (odd) size 0..N. When the even size is recognized, the component at frequency N is discarded by HFT and assumed zero by RFT. When the odd size is provided, the Nth array element is ignored in input by HFT and used to store the Nth frequency component on output, with the converse being enacted by RFT. It must be said that, while the default behaviour for even size 0..N-1 makes HFT and RFT not strictly reciprocal, to zero the Nth component (which has a neither positive nor negative frequency) is an appropriate action when the discrete Fourier transform is used to represent the truncation of a virtually infinite range of frequencies. An odd-sized array must be used in applications where the Nth component needs to be preserved.
The variants IFTU, RFTU, FFTU and HFTU of the above routines achieve a slightly shorter computation time by leaving the output array of the first two (and accepting the input array of the last two) unordered. They are, nonetheless, mutually inverse transforms and can be usefully applied to perform convolutions.
Subroutines IFT, FFT, RFT, HFT are also defined for 2D and 3D arrays, in their simplest implementation as loops of the corresponding 1D transforms.
