fft.h File Reference

Discrete Fourier transform of arbitrary length, in double complex. More...

#include <cmath>
#include <complex>
#include <cstddef>
#include <vector>
#include "line/util/error.h"

Include dependency graph for fft.h:

Go to the source code of this file.

Namespaces
namespace	line

Functions
void	line::dft (std::vector< std::complex< double > > &a, bool inverse)
	In-place DFT of a.

Detailed Description

Discrete Fourier transform of arbitrary length, in double complex.

WHY IT IS HERE. MG1_CR, the Bini-Meini point-wise cyclic reduction behind every M/G/1-type G matrix, evaluates four matrix polynomials at the (nj+1)-th roots of unity, combines them through a point-wise inverse, and interpolates the result back. The transform length is nj+1 where nj+1 doubles from the degree of the block sequence, so it is NOT a power of two in general: a BMAP with three batch sizes gives a degree-4 sequence and a first transform of length 4, but a degree-5 one gives length 6. Both must work, and both must agree with MATLAB's fft/ifft to roundoff, because the reference's truncation test compares the interpolated tail against an absolute epsilon.

The implementation is radix-2 Cooley-Tukey when the length is a power of two and Bluestein's chirp-z otherwise, which reduces the arbitrary length to a power-of-two convolution. Bluestein is exact in the same sense the radix-2 transform is – it is a rearrangement, not an approximation – so the two paths differ only in rounding, and neither introduces the O(n) error a naive quadratic DFT accumulates at n = 2048 (the reference's MaxNumRoot).

CONVENTION, matching MATLAB. dft(a, false) returns sum_k a_k exp(-2 pi i j k / n); dft(a, true) returns (1/n) sum_k a_k exp(+2 pi i j k / n). The inverse therefore carries the 1/n, exactly as ifft does, so a forward followed by an inverse is the identity.

DOUBLE ONLY. The twiddle factors are cos/sin, which the exact and multiprecision number types do not provide; callers that need the transform at another arithmetic must say so rather than silently losing precision here.

Definition in file fft.h.