pm_fftpack Module Reference

This module contains procedures and generic interfaces for computing the Discrete Fourier Transform of a real or complex sequence using a mixed-radix decimation-in-frequency Fast-Fourier Transform.
More...

Data Types
interface	getFactorFFT
	Generate and return the factorization vector factor of the specified input sequence length and the corresponding vector of trigonometric coefficients coef. More...
interface	getFFTF
	Generate and return the Forward Fourier Transform (a.k.a. Fourier Analysis) of a periodic sequence of type complex or real of arbitrary kind type parameter. More...
interface	getFFTI
	Generate and return the Inverse (normalized) Fourier Transform of a periodic sequence of type complex or real of arbitrary kind type parameter. More...
interface	getFFTR
	Generate and return the Reverse (unnormalized) Fourier Transform of a periodic sequence of type complex or real of arbitrary kind type parameter. More...
interface	setFFTF
	Return the Forward Fourier Transform (or equivalently, the Fourier coefficients) of a periodic sequence of type complex or real of arbitrary kind type parameter. More...
interface	setFFTI
	Return the Inverse (normalized) Fourier Transform of a periodic sequence of type complex or real of arbitrary kind type parameter. More...
interface	setFFTR
	Return the Reverse (unnormalized) Fourier Transform of a periodic sequence of type complex or real of arbitrary kind type parameter. More...

Variables
character(*, SK), parameter	MODULE_NAME = "@pm_fftpack"

Detailed Description

This module contains procedures and generic interfaces for computing the Discrete Fourier Transform of a real or complex sequence using a mixed-radix decimation-in-frequency Fast-Fourier Transform.

The discrete Fourier transform (DFT) converts a finite sequence of equally-spaced samples of a function into a same-length sequence of equally-spaced samples of the discrete-time Fourier transform (DTFT), which is a complex-valued function of frequency.

Since DFT deals with a finite amount of data, it can be implemented in computers by numerical algorithms or even dedicated hardware.
These implementations usually employ efficient fast Fourier transform (FFT) algorithms, so much so that the terms FFT and DFT are used interchangeably.

The DFT transforms a sequence of \(N\) complex numbers \(\{ x_j \} := x_0, x_1, \ldots, x_{N − 1}\) into another sequence of complex numbers, \(\{z_k\} := z_0, z_1, \ldots, z_{N - 1}\) which is defined by,

\begin{eqnarray} \large x_k &=& \sum_{j=0}^{N-1} z_j \cdot e^{-\frac {i 2\pi}{N}jk} ~, \\ &=& \sum_{j=0}^{N-1} z_j \cdot \left[\cos\left(\frac{2\pi}{N}jk\right) - i \cdot \sin\left(\frac{2 \pi}{N}jk\right)\right], \end{eqnarray}

where \(i\) represents the imaginary unit.
The naive evaluation of the DFT is a matrix-vector multiplication that costs \(\mathcal{O}(n^2)\) operations for \(N\) data-points.
Fast Fourier Transform (FFT) algorithms use a divide-and-conquer strategy to factorize the matrix into smaller sub-matrices, corresponding to the integer factors of the length \(N\).
If \(N\) can be factorized into a product of integers \(f_1 f_2 \ldots f_m\), then the DFT can be computed in \(\mathcal{O}(n \sum f_i)\) operations.
For a radix-2 FFT this gives an operation count of \(\mathcal{O}(n \log_2 n)\).

For every Forward DFT (setFFTF), which expresses the data sequence in terms of frequency coefficients), there is a corresponding Reverse (Backward) DFT.
The Reverse DFT (setFFTR) takes the output of the Forward DFT and returns the original data sequence given to the Forward DFT, multiplied by the length of the sequence,

\begin{eqnarray} \large N z_j &=& \sum_{k=0}^{N-1} x_k \cdot e^{\frac {i 2\pi}{N}jk} ~, \\ &=& \sum_{k=0}^{N-1} x_k \cdot \left[\cos\left(\frac{2\pi}{N}jk\right) + i \cdot \sin\left(\frac{2 \pi}{N}jk\right)\right], \end{eqnarray}

A third definition, the Inverse DFT (setFFTI), further normalizes the Reverse DFT and returns,

\begin{eqnarray} \large z_j &=& \frac{1}{N} \sum_{k=0}^{N-1} x_k \cdot e^{\frac {i 2\pi}{N}jk} ~, \\ &=& \frac{1}{N} \sum_{k=0}^{N-1} x_k \cdot \left[\cos\left(\frac{2\pi}{N}jk\right) + i \cdot \sin\left(\frac{2 \pi}{N}jk\right)\right], \end{eqnarray}

such that the resulting normalized Reverse (Backward) DFT or Inverse DFT becomes a true inverse of the Forward DFT.
The Reverse DFT becomes relevant when the overall scale of the result is unimportant.
In such cases, it is convenient to use the Reverse (Backward) DFT instead of the Inverse DFT to avoid the redundant divisions in the normalization step.
A 64-bit real division is typically on the order of \(\sim20\) CPU cycles on the contemporary hardware.
Avoiding the redundant normalizations can therefore become a significant saving if the DFT is to computed repeatedly for an array of millions of elements.

The Cooley–Tukey algorithm

The Cooley–Tukey algorithm, named after J. W. Cooley and John Tukey, is the most common fast Fourier transform (FFT) algorithm.
It re-expresses the discrete Fourier transform (DFT) of an arbitrary composite size \(N = N_{1} N_{2}\) in terms of \(N_1\) smaller DFTs of sizes \(N_2\), recursively, to reduce the computation time to \(\mathcal{O}(N log N)\) for highly composite \(N\) (smooth numbers).
Because of the algorithm importance, specific variants and implementation styles have become known by their own names, as described below.
Because the Cooley–Tukey algorithm breaks the DFT into smaller DFTs, it can be combined arbitrarily with any other algorithm for the DFT.
For example, The Rader or Bluestein algorithm can be used to handle large prime factors that cannot be decomposed by Cooley–Tukey, or the prime-factor algorithm can be exploited for greater efficiency in separating out relatively prime factors.

The algorithm, along with its recursive application, was invented by Carl Friedrich Gauss.
Cooley and Tukey independently rediscovered and popularized it 160 years later.

Intuition

The Cooley–Tukey algorithms recursively re-express a DFT of a composite size \(N = N_1 N_2\) as:

1. Perform \(N_1\) DFTs of size \(N_2\).
2. Multiply by complex roots of unity (often called the twiddle factors).
3. Perform \(N_2\) DFTs of size \(N_1\).

Typically, either \(N_1\) or \(N_2\) is a small factor (not necessarily prime), called the radix (which can differ between stages of the recursion).
If \(N_1\) is the radix, it is called a decimation in time (DIT) algorithm, whereas if \(N_2\) is the radix, it is decimation in frequency (DIF) algorithm, also called the Sande–Tukey algorithm.

Note

The routines of this module are generic mixed-radix implementation, meaning that the can be used to compute the forward and reverse FFT of arbitrary-length real and complex data sequences.
By contrast, radix-2 FFT routines (such as those of the Numerical Recipes), require a potentially-padded sequence (with trailing zeros) such that the length of the input sequence is always a power of \(2\).
The mixed-radix algorithms are generally faster than the radix-2 algorithms at the cost of requiring an extra memory storage of the same length as the input sequence.
The mixed-radix algorithm uses optimized small length sub-transform FFTs which are combined to create larger FFT of the sequence.
FFTPACK implements efficient sub-transform for factors of 2, 3, 4, and 5.
The computations for other factors fall back to a general length- \(N\) algorithm which uses Singleton method for efficiently computing a DFT.
This module is \(\mathcal{O}(N^2)\), slower than an explicitly dedicated factor implementation. However, it works for arbitrary length sub-transforms.
Arbitrary-length sub-transforms are factorized as much as possible.
For example, a sub-sequence length of \(143\) will be factorized into \(11\times13\).
Large prime factors (e.g., \(n = 2\times 3\times 3,351,773\)) are the worst case scenario, because they cannot be further factorized.
The large prime factors in sequence should be avoided as much as possible because their \(\mathcal{O}(N^2)\) computational cost scaling will dominate the run-time.
The mixed-radix initialization procedures under the generic interface getFactorFFT return the list of factors chosen for a given input sequence length \(N\).
The output factors from these procedures can be checked to ensure an efficient factorization and to estimate the run-time.
The run-time of the FFT algorithms of this module roughly scale as \(N\sum f_i\), where the \(f_i\) are the factors of the sequence length \(N\).
If specific data lengths appear in the FFT problem that cannot be efficiently factorized using the existing small-prime factor implementations, the performance can be improved by explicitly implementing the algorithms for the specific factors of interest.

See this catalog for a list of prime numbers.

Attention

There are two possible choices for the sign of the exponential in the forward and inverse FFT equations in the above.
While both conventions are commonly used, the FFTPACK library (as given here) uses the above convention (a negative exponential for the forward transform).
A prominent implementation that uses the opposite convention is provided by the Numerical Recipes.

Remarks

This reimplementation of the original FFTPACK library achieves the following goals:

1. Extension of the functionality of the library to arbitrary complex and real precision kinds (e.g., any supported by the processor (e.g., RK, RK32, RK64, or RK128), any supported by the processor (e.g., CK, CK32, CK64, or CK128)).
   
2. Usage and flexibility improvements to the syntax and interface of the library.
   
3. Performance improvements to the original library.
   
4. Safety improvements to the library, particularly,
   
   1. the removal of all assumed-shape dummy arguments,
      
   2. the removal of all implicit conversions between complex and real actual and dummy arguments,
      
   3. the removal of all implicit conversions between complex and integer actual and dummy arguments,
      
   4. the addition of abundant runtime array-bound checks to reduce misspecified input arguments.
      
5. Removal of all fixed-format FORTRAN77 syntax and their replacement with equivalent free-format syntax.
   
6. Removal of all obsolescent language features such as goto statements from the original library.
   
7. Significant reduction in the library size with the help of generics and macros.
   

Note

An acyclic time series may need to be tapered by applying split-cosine-bell tapering to the time series prior to a Fourier Analysis.
See Chapter 5 of the book Fourier Analysis of Time Series, Peter Bloomfield, Wiley-Interscience, 1976.

See also

FFTPACK, the original mixed-radix FFT library in FORTRAN77 by Paul Swarztrauber of NCAR.
FFTPACK 5.1, an unofficial mirror of FFTPACK 5.1 by Paul Swarztrauber of NCAR and Richard A. Valent (2011).
FFTPACK 5.1, an unofficial reassembly of FFTPACK 5.1 by the original authors.
Modernized FFTPACK, a modernized version of the original FFTPACK library in modern Fortran by the Fortran-lang community.
The GNU FFTPACK in C, an incomplete translation of the original FORTRAN77 FFTPACK library in C language.
Paul Swarztrauber, 1982, Vectorizing the Fast Fourier Transforms, Parallel Computations, G. Rodrigue, ed., Academic Press, New York
Paul Swarztrauber, 1984, Fast Fourier Transforms Algorithms for Vector Computers, Parallel Computing, pp.45-63
Clive Temperton, 1983, Self-sorting Mixed-radix FFTs, Journal of Computational Physics, 52, 340-350

Test:

test_pm_fftpack

Todo:

Critical Priority: The sine and cosine forward and reverse transforms must be added to this module.

Final Remarks ⛓

---

If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

1. If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
   
2. If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.
   

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Copyright

Computational Data Science Lab This module is a complete modernization and significant extension of the original netlib.org public-domain FFTPACK library in FORTRAN77 programming language.
If you use any or parts of this module, you should also cite the following reference in addition to the relevant ParaMonte library reference.

• P. N. Swarztrauber, Vectorizing the FFTs, in Parallel Computations (G. Rodrigue, ed.), Academic Press, 1982, pp. 51-83.
  

                        FFTPACK

  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

                    version 4  april 1985

       a package of fortran subprograms for the fast fourier
        transform of periodic and other symmetric sequences

                           by

                    paul n swarztrauber

    national center for atmospheric research  boulder,colorado 80307

     which is sponsored by the national science foundation

  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

Author:

Amir Shahmoradi, April 25, 2015, 2:21 PM, National Institute for Fusion Studies, The University of Texas Austin

Variable Documentation

MODULE_NAME

character(*, SK), parameter pm_fftpack::MODULE_NAME = "@pm_fftpack"

Definition at line 212 of file pm_fftpack.F90.