File FFT.h

↰ Parent directory (ADS/numeric)

The fast Fourier transforms, beyond the book (4.3.3 and 4.6.4 of Knuth): the discrete transform :math:` hat a_k = sum_{j<n} a_j omega_n^{j k} , :math: omega_n = e^{-2 pi i / n} , of n complex numbers in double or single precision, for every n, in :math: O( n log n ) ` operations instead of the :math:` n^2 ` products of the definition, dft(); its inverse :math:` a_j = n^{-1} sum_k hat a_k omega_n^{-j k} `; the transforms of real data, of the cosines and of the sines, many transforms at once and those in more dimensions, and the product of polynomials.

Definition (ADS/numeric/FFT.h)

Detailed Description

The codelets. A transform of N points, N a power of 2 up to 64 or :math:` 2^a 3^b 5^c ` up to 25 (only 3 and 5 under ADS_FFT_LEAN), is one codelet, unrolled when compiled with its items in registers and its powers of :math:` omega_N ` the constants of a table: by decimation in time in radix 4 for the powers of 2, :math:`[ hat a_{k + u N/4} = sum_{t<4} ( -i )^{t u} omega_N^{t k}

hat a^{(t)}_k , qquad k < N / 4 , ` :math:` hat a^{(t)} ` the transform of the items :math:` t , t + 4 , ldots

; a size N = A B with A and B coprime (A the power of the smallest prime of N) by the prime-factor algorithm of Good and Thomas, with no power of :math: omega ` between its two stages: the input :math:` ( a B + b A ) bmod N ` and the output :math:` ( u B ( B^{-1} bmod A ) + v A ( A^{-1} bmod B ) ) bmod N ` of indices :math:` a , u < A ` and :math:` b , v < B ` make the transform one in 2 dimensions, :math:`[ hat a_{k( u , v )} = sum_{b<B} omega_B^{b v} sum_{a<A}

omega_A^{a u} , a_{( a B + b A ) bmod N} , ` the transforms of length A of B columns and of length B of A rows; the other sizes in the smallest prime radix p of N. A prime p is done by pairs, with :math:` s_j = a_j + a_{p-j} , :math: d_j = a_j - a_{p-j} ` and :math:` h = ( p - 1 ) / 2 : :math:`[ hat a_{k} , hat a_{p-k} = a_0 + sum_{j le h} s_j cos frac{2 pi j k}{p} mp i sum_{j le h} d_j sin frac{2 pi j k}{p} , ` :math: 2 h^2 ` real products by component instead of :math:` ( p - 1 )^2 `; the primes 7 , 11 , 13 have such codelets, and every other prime up to 101 the same sums in loops.

The four steps. Above the codelets, n = r m points are the matrix of m rows and r columns, :math:` a_{r i + j} `, and, by the factorization of Cooley and Tukey, :math:`[ hat a_{k + m q} = sum_{j<r} omega_r^{j q} , omega_n^{j k}

sum_{i<m} a_{r i + j} omega_m^{i k} , qquad k < m , ; q < r : ` the transforms of length m of the columns, the products by :math:`

omega_n^{j k} , and those of length r of the rows, whose outputs are those of the transform in their order, with no reversal of the digits (the four steps of Bailey, the last one, the transposition, done by the order in which the rows are written). Each step works on W columns, or rows, at a time, one in each lane of a vector of W numbers (4 doubles or 8 floats on AVX, half on SSE2 and NEON), the real and the imaginary parts in separate vectors: a codelet on such items does W transforms with the instructions of one, and the blocks of W by W between the two steps are transposed in registers, as in the short vector codes of Franchetti and Pueschel; the numbers alternated in memory and apart in the registers, as the mixed layout of Popovici, Franchetti and Low. The vectors are those of the extensions of GCC and Clang, the registers of SSE2 and AVX by the intrinsics of x86 with the compilers that have no such extensions (MSVC), or arrays of W numbers elsewhere or with ADS_FFT_SIMD defined as 0. The columns and the rows longer than the codelets go by the rounds of Stockham, each a codelet of radix :math: R_t ` with its powers of :math:` omega `, :math:`[ y_{q + s ( R p + j )} = omega_{m R}^{p j} sum_{i<R} x_{q + s ( p +

i m )} , omega_R^{i j} , ` as few rounds as the radices of the codelets up to 32 allow, the radices as equal as can be; 4 of them at a time, gathered from lines of the cache brought in some rows ahead. The rows have about 32 points up to :math:` 2^{15} ` points, 64 at :math:` 2^{12} ` and :math:` 2^{16} ` and about :math:` sqrt n ` above: the divisors r of n by their distance from those lengths, those that make r and m multiples of W (otherwise the last block of the lanes is shorter) and whose rows and columns have a codelet first (the sizes :math:` 2^a 3^b 5^c ` up to 60 there); up to :math:` 2^{17} ` points the first 4 of them are timed when the plan is made and the fastest kept, a choice among a few, as FFTW_MEASURE makes among many (ADS_FFT_TUNE sets the limit, 0 none). Below :math:` 2 W^2 ` points the lanes are half as wide. Up to 1024 points a power of 2 (2048 in single precision), and at 48 and 96, the four steps are unrolled into one function whose columns stay in an array of n / W items, the transform of 16 points in 2 codelets of 4; so are the transforms in 2 dimensions of the powers of 2 up to 32 x 32. In double precision on the vectors of AVX the many transforms of 4 to 32 points (those of fftMany(), the rows and the columns of fft2() and fftn()) go two at a time, one in each half of a vector of the numbers alternated, with no transposition (a product by i costs one swap of the parts), and a power of 2 n of 8 or 16 points by the two halves :math:` E ` and :math:` O ` of the even and the odd points so made, then :math:` hat a_k = E_k + omega_n^k

O_k , :math: hat a_{k + n/2} = E_k - omega_n^k O_k . The powers of :math: omega ` are each computed alone, the angle taken to :math:` [ 0 , pi / 4 ] ` by the symmetries of the circle, in integers, so that each has the error of a sine or a cosine of a small angle; from :math:` 2^{16} ` points on, those between the steps are the products of two tables of about :math:` n / W ` and r items, rather than one of n.

The large primes. A prime n above 13 goes by the algorithm of Rader, the transform of the nonzero indices as a cyclic convolution of n - 1 points: for a generator g of the integers mod n, :math:`[ hat a_{g^{-q}} = a_0 + sum_{p < n-1} a_{g^p} , omega_n^{g^{p -

q}} , qquad hat a_0 = sum_j a_j , ` by two transforms of n - 1 points (the transform of :math:` omega_n^{g^{-q}}

` made once); a size whose prime above 101 is not n, or a prime whose n - 1 costs more, by that of Bluestein, from :math:` j k = ( j^2 + k^2 - ( k - j )^2 ) / 2 `, :math:`[ hat a_k = w_k sum_j ( a_j w_j ) overline{w_{k-j}} , qquad w_j =

omega_{2n}^{j^2} , ` a convolution of :math:` M ge 2 n - 1 ` points of the form :math:` 2^a 3^b

5^c , by two transforms of M (that of :math: overline w ` made once). The choice among the rounds, Rader and Bluestein is by a cost of about :math:` n sum_t log_2 R_t ` for the rounds with codelets and :math:` n p / 2 ` for a prime radix p in loops. The inverse transforms by Rader and Bluestein are the conjugates of the forward ones of the conjugates.

Real data, cosines and sines. The transform of n real numbers, n even, is one complex transform of h = n / 2 points, :math:` z_j = x_{2j} + i x_{2j+1} , and a pass: with :math: E_k = ( Z_k + overline{Z_{h-k}} ) / 2 ` and :math:` O_k = ( Z_k - overline{Z_{h-k}} ) / 2 i ` (the transforms of the even and of the odd items), :math:` hat x_k = E_k + omega_n^k O_k ` and :math:` hat x_{h-k} = overline{E_k - omega_n^k O_k} ; the inverse runs the same steps backwards. Up to 16 points the complex transform is one codelet whose outputs stay in registers for the pass; the inverse of 4098 to 131072 points on the four steps does its pass in their first step, the points :math: Z_K ` made from :math:` hat x_K ` and :math:` hat x_{h-K} ` as the columns read them (and convolve() there the products of the two transforms); elsewhere the pass runs on its own on the numbers alternated, VL pairs to a vector with 3 shuffles. The transforms of the cosines and of the sines are those of FFTW (REDFT00, REDFT10, REDFT01, REDFT11 and RODFT00 to RODFT11), with no normalization; the type 2 by the algorithm of Makhoul, the real transform V of :math:` v_j = x_{2j} , v_{n-1-j} = x_{2j+1} ` and :math:` y_k = 2 Re( omega_{4n}^k V_k ) , y_{n-k} = -2 Im( omega_{4n}^k V_k ) ; the type 3 by the same steps backwards; the type 4, for n even, by the complex transform of n / 2 points of :math: ( x_{2j} + i x_{n-1-2j} ) omega_{8n}^{4j+1} ` times :math:` omega_{2n}^k `, the outputs 2 k and n - 1 - 2 k its real part and minus its imaginary part, times 2; the type 1 of both by the real transform of the even, or odd, extension of 2 ( n - 1 ), or 2 ( n + 1 ), points; the sines of types 2 to 4 by the cosines, the items with the signs alternated or in the reverse order. The type 3 is the inverse of the type 2, and the type 4 of itself, times 2 n.

The cost. A transform of :math:` n = 2^K ` points does about :math:` 4.25 n K ` real operations, the radix 4 saving 15% of those of radix 2; its data cross the memory in 2 passes (the columns, the rows) plus one per round of Stockham of the columns or rows longer than a codelet, the four steps keeping a buffer of n + O( sqrt n ) items per thread. The plans, the tables and the buffers are kept per thread for each size met, made at the first transform of a size, in :math:` O( n ) ` evaluations of sines and cosines ( :math:` O( n log n ) ` for Rader and Bluestein, whose plans also hold a transform) and, up to :math:` 2^{17} ` points, the 12 transforms that time the factorizations.

The rounding. With powers of :math:` omega ` computed to a relative error :math:` mu , a transform by L rounds of radix 2 has :math: | tilde y - y |_2 / | y |_2 le L eta / ( 1 - L eta ) , :math: eta = mu + gamma_4 ( sqrt 2 + mu ) , :math: gamma_4 = 4 u / ( 1 - 4 u ) ` and u the unit roundoff (Higham, Theorem 24.2); with :math:` mu ` near u the bound is about :math:` 6.7 u log_2 n , while the error of the rounding observed on random data grows as :math: u sqrt{ log n } ` (Schatzman). Bluestein adds the error of a convolution of M points, a few times larger. The tests check every size up to 300 and the sizes of each plan against the definition in long double.

The transform turns a convolution into a product item by item; convolve() multiplies two polynomials so, by the real transforms of the smallest even size :math:` 2^a 3^b 5^c ` that holds their product, or of the power of 2 if that size is not below 3 / 4 of it, and the product of integers of many digits [see File BigInteger.h] is a convolution of their digits followed by the carries.

The contracts [see File Contract.h]: the sizes of the arguments of irfft(), the types of dct() and dst(); the inverse transform of the transform gives the vector back, up to the rounding, which the tests check. The functions are defined in FFT.tpp.

    1. Knuth, The Art of Computer Programming, vol. 2, 3rd ed., Addison-Wesley, 1998, sections 4.3.3 and 4.6.4

    1. Cooley, J. W. Tukey, An algorithm for the machine calculation of complex Fourier series, Mathematics of Computation 19(90), 297-301, 1965, doi:10.1090/S0025-5718-1965-0178586-1

    1. Good, The interaction algorithm and practical Fourier analysis, Journal of the Royal Statistical Society, Series B 20(2), 361-372, 1958, doi:10.1111/j.2517-6161.1958.tb00300.x

    1. Thomas, Using computers to solve problems in physics, in W. F. Freiberger, W. Prager (eds.), Applications of Digital Computers, Ginn, Boston, 1963 (no doi)

    1. Singleton, An algorithm for computing the mixed radix fast Fourier transform, IEEE Transactions on Audio and Electroacoustics 17(2), 93-103, 1969, doi:10.1109/TAU.1969.1162042

  1. Frigo, S. G. Johnson, The design and implementation of FFTW3, Proceedings of the IEEE 93(2), 216-231, 2005, doi:10.1109/JPROC.2004.840301

    1. Bailey, FFTs in external or hierarchical memory, The Journal of Supercomputing 4(1), 23-35, 1990, doi:10.1007/BF00162341

  1. Franchetti, M. Pueschel, Short vector code generation for the discrete Fourier transform, Proceedings of the International Parallel and Distributed Processing Symposium (IPDPS), 2003, doi:10.1109/IPDPS.2003.1213153

    1. Popovici, F. Franchetti, T. M. Low, Mixed data layout kernels for vectorized complex arithmetic, IEEE High Performance Extreme Computing Conference (HPEC), 1-7, 2017, doi:10.1109/HPEC.2017.8091024

    1. Stockham, High-speed convolution and correlation, AFIPS Spring Joint Computer Conference, 229-233, 1966, doi:10.1145/1464182.1464209

  1. Van Loan, Computational Frameworks for the Fast Fourier Transform, SIAM, 1992, doi:10.1137/1.9781611970999

    1. Rader, Discrete Fourier transforms when the number of data samples is prime, Proceedings of the IEEE 56(6), 1107-1108, 1968, doi:10.1109/PROC.1968.6477

    1. Bluestein, A linear filtering approach to the computation of discrete Fourier transform, IEEE Transactions on Audio and Electroacoustics 18(4), 451-455, 1970, doi:10.1109/TAU.1970.1162132

    1. Cooley, P. A. W. Lewis, P. D. Welch, The fast Fourier transform algorithm: programming considerations in the calculation of sine, cosine and Laplace transforms, Journal of Sound and Vibration 12(3), 315-337, 1970, doi:10.1016/0022-460X(70)90075-1

    1. Sorensen, D. L. Jones, M. T. Heideman, C. S. Burrus, Real-valued fast Fourier transform algorithms, IEEE Transactions on Acoustics, Speech, and Signal Processing 35(6), 849-863, 1987, doi:10.1109/TASSP.1987.1165220

  1. Makhoul, A fast cosine transform in one and two dimensions, IEEE Transactions on Acoustics, Speech, and Signal Processing 28(1), 27-34, 1980, doi:10.1109/TASSP.1980.1163351

  1. Wang, Fast algorithms for the discrete W transform and for the discrete Fourier transform, IEEE Transactions on Acoustics, Speech, and Signal Processing 32(4), 803-816, 1984, doi:10.1109/TASSP.1984.1164399

    1. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, chapter 24, doi:10.1137/1.9780898718027

    1. Schatzman, Accuracy of the discrete Fourier transform and the fast Fourier transform, SIAM Journal on Scientific Computing 17(5), 1150-1166, 1996, doi:10.1137/S1064827593247023

Donato Meoli

Dipartimento di Informatica Universita’ di Pisa

by Donato Meoli

Includes

Included By

Namespaces

Functions

Defines

Typedefs