This is the top of the page.

Please note that JavaScript and style sheet are used in this website,
Due to unadaptability of the style sheet with the browser used in your computer, pages may not look as original.
Even in such a case, however, the contents can be used safely.

MathKeisan®

1. Home
2. User's Guide
3. User's Guide: FFT & PARFFT

# User's Guide: FFT & PARFFT

## Introduction

FFT is a library which computes Fast Fourier Transforms. The FFT is an algorithm for calculating Discrete Fourier Transforms (DFT). The DFT of the complex numbers is the complex numbers , given by

$$y_k = \sum_{j=0}^{n-1} x_j e^{\frac{-2 \pi i}{n} jk} \; \; \; \; \; \; k = 0, 1, \ldots, n-1$$

In this document, . The reverse transform is given by

$$x_j = \frac{1}{n} \sum_{k=0}^{n-1} y_k e^{\frac{ 2 \pi i }{n} jk} \; \; \; \; \; j = 0, 1, \ldots, n-1$$

In the above equations, there is no scaling on the forward transform, and the reverse transform is scaled by . An alternative convention is to scale the forward transform by and have no scaling on the reverse transform, or to scale both transforms by .

## User Interface

Software packages for FFTs differ in functionality, and there is no standard user interface. The SX MathKeisan FFTs contain the same functionality as Cray® LibSci™ 3.1 and HP MLIB VECLIB [3]. The interfaces are the same except that the size of work arrays in MathKeisan LibSci™ differ from those in Cray's LibSci™. The FFT Man Pages have information on work array sizes.

MathKeisan FFT supports these types of transforms:

• One-dimensional (1-d) FFT
• Multiple 1-d FFT
• Two-dimensional (2-d) FFT
• Three-dimensional (3-d) FFT

They are all available for these data type mappings:

• Complex-to-complex
• Real-to-complex
• Complex-to-real

User interface information is available from several sources:

• The FFT Routine List contains one-line summaries of LibSci™ and VECLIB subroutines.
• The FFT Man Pages have user interface information for individual subroutines, including LibSci™ array sizes. To view information on subroutine ZFFTMLT from the command-line, type man zfftmlt. If you have difficulty viewing man pages, please see the Man Pages Chapter.
• There is user interface information for VECLIB FFTs in the HP MLIB User's Guide [3].
• The Linking Chapter tells how to link with the FFT and PARFFT libraries.

## PARFFT

PARFFT is an OpenMP parallel version of FFT. It has the same user interface as FFT, so any code that is linked to FFT can alternatively be linked to PARFFT. If the environment variable OMP_NUM_THREADS is set to np, then PARFFT will run on np threads. If OMP_NUM_THREADS is unset, PARFFT will run on mp threads, where mp is the maximum number of processes in the resource group or node.

## A Brief Sketch of the FFT Algorithm

If we let

$$\omega^k_n = e^{\frac{-2 \pi i}{n}k} = \cos \frac{2 \pi }{n}k - i \sin\frac{2 \pi }{n}k$$

then the DFT is given by the matrix-vector multiplication:

$$\left[ \begin{array}{lllcl} \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0 \\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1} \\ \omega_n^0 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)} \\ \vdots & \vdots & \vdots & & \vdots \\ \omega_n^0 & \omega_n^{(n-1)} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \\ \end{array} \right] \left[ \begin{array}{c} x_0 \\ x_1 \\ x_3 \\ \vdots \\ x_{n-1} \\ \end{array}\right] = \left[ \begin{array}{c} y_0 \\ y_1 \\ y_3 \\ \vdots \\ y_{n-1} \\ \end{array}\right]$$

The floating-point operation count for the DFT is if the calculation of the powers of is ignored.

In 1965 Cooley and Tukey [1] developed the FFT algorithm to calculate the DFT. When is a power of 2, the FFT algorithm has the greatly reduced operation count of . For general , the operational count is . Unknown to Cooley and Tukey at the time, the FFT algorithm was proposed earlier by Gauss in 1805, and the work of Gauss was published posthumously in 1866 [2].

A derivation of the Cooley-Tukey algorithm is shown in Van Loan [4]. To get an idea of the algorithm, consider the case where is a power of 2. Start by substituting into the DFT equation to get:

$$y_k = \sum^{n-1}_{j=0} x_j \omega_n^{jk} \; \; \; \; \; k = 0, 1, \ldots, n-1$$

Now set . For , split into even components and odd components . Substituting into the equation above gives:

$$y_k = \sum^{m-1}_{j=0} \underline{x}_j \omega_n^{2jk} + \sum^{m-1}_{j=0} \overline{x}_j \omega_n^{(2j+1)k} \; \; \; \; \; \; \; k = 0, 1, \ldots, n-1$$

Use the symmetry property to get:

$$y_k = \sum^{m-1}_{j=0} \underline{x}_j \omega_{m}^{jk} + \omega_n^{k} \sum^{m-1}_{j=0} \overline{x}_j \omega_{m}^{jk} \; \; \; \; \; k = 0, 1, \ldots, n-1$$

Next, split into the first entries and the second entries as follows:

$\begin{eqnarray} y_k & = & \sum^{m-1}_{j=0} \underline{x}_j \omega_{m}^{jk} + \omega_n^{k} \sum^{m-1}_{j=0} \overline{x}_j \omega_{m}^{jk} \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; k = 0, 1, \ldots, m-1 \nonumber \\ y_{k+m} & = & \sum^{m-1}_{j=0} \underline{x}_j \omega_{m}^{j(k+m)} + \omega_n^{k+m} \sum^{m-1}_{j=0} \overline{x}_j \omega_{m}^{j(k+m)} \; \; \; \; \; k = 0, 1, \ldots, m-1 \nonumber \end{eqnarray}$

For the equation use the symmetry properties and to get:

$\begin{eqnarray} y_k & = & \sum^{m-1}_{j=0} \underline{x}_j \omega_{m}^{jk} + \omega_n^{k} \sum^{m-1}_{j=0} \overline{x}_j \omega_{m}^{jk} \; \; \; \; \; k = 0, 1, \ldots, m-1 \nonumber \\ y_{k+m} & = & \sum^{m-1}_{j=0} \underline{x}_j \omega_{m}^{jk} - \omega_n^{k} \sum^{m-1}_{j=0} \overline{x}_j \omega_{m}^{jk} \; \; \; \; \; k = 0, 1, \ldots, m-1 \nonumber \end{eqnarray}$

We have taken the original DFT of size and an operation count and replaced it with 2 DFTs of size and operation count . This can be continued for steps to replace the original DFT of length by DFTs of length 1. This is called a radix-2 FFT. The algorithm can be extended beyond powers of 2, to general . The splitting process above corresponds to the Decimation In Time (DIT) Cooley-Tukey algorithm. There are a total of 8 splittings. The splitting used in the MathKeisan FFT library is the Decimation In Frequency (DIF) Stockham algorithm. A derivation of this algorithm is shown in Van Loan [4]. The DIF Stockham algorithm has an operation count equal to that of Cooley-Tukey, but the order in which the operations are performed is different.

## Multiple 1-d FFTs and Vector Length

The FFT calculation proceeds as a series of steps of radix , where the are generally small prime numbers from the factorization .

The vector length during the calculation of the radix step will be ; at the radix step it will be the product ; at the radix step it will be:

$$\prod_{j=1}^{p}n_j$$

If is small, then the vector lengths in the FFT calculation will be small.

Some applications allow for the calculation of a number of 1-d FFTs of the same length simultaneously. If the length of the FFTs is , and the number of FFTs is , then the vector length at the radix step of the calculation will be:

$$h \prod_{j=1}^{p}n_j$$

This larger vector length gives performance advantages for multiple 1-d over 1-d FFTs.

The multiple 1-d FFT can be written:

$$y_{kl} = \sum_{j=0}^{n-1} x_{jl} \omega_n^{jk} \; \; \; \; \; k = 0, 1, \ldots, n-1 \; \; \; \; \; l = 0, 1, \ldots, h-1$$

Here it is understood that and are 2-dimensional arrays of complex numbers.

## 2-d and 3-d FFTs

The 2-d FFT of the array of complex numbers , for and is the array of complex numbers given by

$\begin{eqnarray} y_{kp} & = & \sum_{l=0}^{h-1} \sum_{j=0}^{n-1} x_{jl} \omega_h^{lp} \omega_n^{jk} \; \; \; \; \; k = 0, 1, \ldots, n-1 \; \; \; \; \; p = 0, 1, \ldots, h-1 \nonumber \\ & = & \sum_{l=0}^{h-1} \omega_h^{lp} \sum_{j=0}^{n-1} x_{jl} \omega_n^{jk} \; \; \; \; \; k = 0, 1, \ldots, n-1 \; \; \; \; \; p = 0, 1, \ldots, h-1 \nonumber \end{eqnarray}$

Now write:

$$z_{kl} = \sum_{j=0}^{n-1} x_{jl} \omega_n^{jk} \; \; \; \; \; k = 0, 1, \ldots, n-1 \; \; \; \; \; l = 0, 1, \ldots, h-1$$

Then is given by

$$y_{kp} = \sum_{l=0}^{h-1} z_{kl} \omega_h^{lp} \; \; \; \; \; k = 0, 1, \ldots, n-1 \; \; \; \; \; p = 0, 1, \ldots, h-1$$

Here the 2-d FFT is replaced by the two multiple 1-d FFTs of the equations for and above.

By similar methods the 3-d FFT can be replaced by three multiple 1-d FFTs. Internally in MathKeisan subroutines 2-d and 3-d FFTs are calculated by calls to multiple 1-d computational kernels. This is known as the row-column algorithm.

## Real-to-Complex and Complex-to-Real FFTs

The real-to-complex FFT has symmetry properties which give about a two times improvement in operation count and memory usage over the complex-to-complex FFT. To show the symmetry properties, write the DFT in terms of the trigonometric functions and to get:

$$y_k=\sum_{j=0}^{n-1} x_j \cos \frac{2\pi }{n} jk - i \sum_{j=0}^{n-1} x_j \sin \frac{2\pi }{n} jk \; \; \; \; \; k = 0, 1, \ldots, n-1$$

If is real, then is complex conjugate symmetric about . This allows for storing only half of the values of . The value is refered to as the DC frequency. If is even, then the value is refered to as the Nyquist frequency. By the identities and it can be seen that and are real.

For the real-to-complex FFT with real input vector of size , to obtain the speed and memory advantages, the output vector should be complex with size . For the complex to real FFTs if the input vector is complex with size , the output vector should be real with size . The reason for the size rather than is because the redundant zero parts of the DC and Nyquist frequencies are stored. For interface information on the size requirements for the real and complex input and output vectors, see the man pages.

The multiple 1-d, 2-d and 3-d FFT subroutines have data stored in two or three dimensional arrays. If the leading dimension of an array is a high powers of 2, bank conflicts will result in poor performance. It is best to set the leading dimension of the arrays to an odd numbers to avoid bank conflicts.

## FFT Routine List

Name† Prefixes Description ?FFT CC ZZ 1-d Complex-Complex FFT. ?FFT C Z 1-d Complex-Complex FFT. ?FFT2 C Z 1-d Complex-Complex FFT. ?FFT SC DZ 1-d Real-Complex FFT. ?FFT CS ZD 1-d Complex-Real FFT. ?FFT2 RC DZ 1-d Real-Complex FFT. ?FFT2 CR ZD 1-d Complex-Real FFT. ?FFT2D CC ZZ 2-d Complex-Complex FFT. ?FFT2D C Z 2-d Complex-Complex FFT. ?FFT2D SC DZ 2-d Real-Complex FFT. ?FFT2D CS ZD 2-d Complex-Real FFT. ?FFT3D CC ZZ 3-d Complex-Complex FFT. ?FFT3D C Z 3-d Complex-Complex FFT. ?FFT3D SC DZ 3-d Real-Complex FFT. ?FFT3D CS ZD 3-d Complex-Real FFT. ?FFTM CC ZZ Multiple 1-d Complex-Complex FFT. M?FFT C Z Multiple 1-d Complex-Complex FFT. ?FFTMLT C Z Multiple 1-d Complex-Complex FFT. ?FFTM SC DZ Multiple 1-d Real-Complex FFT. ?FFTM CS ZD Multiple 1-d Complex-Real FFT. ?FFTMLT R D Multiple 1-d Real-Complex, Complex-Real FFT. ?FTFAX F D Generate FFT factors and Calculate Trig tables, Real. ?FTFAX C Z Generate FFT factors and Calculate Trig tables, Complex. XERFFT FFT error handler. ?1DFFT C Z 1-d Complex-Complex FFT, input/output Complex. ?1DFFT S D 1-d Complex-Complex FFT, input/output Real. ?RC1FT C Z 1-d Real-Complex, Complex-Real FFT, input/output Complex. ?RC1FT S D 1-d Real-Complex, Complex-Real FFT, input/output Real. ?2DFFT C Z 2-d Complex-Complex FFT, input/output Complex. ?2DFFT S D 2-d Complex-Complex FFT, input/output Real. ?RC2FT C Z 2-d Real-Complex, Complex-Real FFT, input/output Complex. ?RC2FT S D 2-d Real-Complex, Complex-Real FFT, input/output Real. ?3DFFT C Z 3-d Complex-Complex FFT, input/output Complex. ?3DFFT S D 3-d Complex-Complex FFT, input/output Real. ?RC3FT C Z 3-d Real-Complex, Complex-Real FFT, input/output Complex. ?RC3FT S D 3-d Real-Complex, Complex-Real FFT, input/output Real. ?FFTS C Z Multiple 1-d Complex-Complex FFT, input/output Complex. ?FFTS S D Multiple 1-d Complex-Complex FFT, input/output Real. ?RCFTS C Z Multiple 1-d Real-Complex, Complex-Real FFT, input/output Complex. ?RCFTS S D Multiple 1-d Real-Complex, Complex-Real FFT, input/output Real.