Fast Fourier transform

The fast Fourier transform or FFT, is the efficient algorithm for the numerical evaluation of the Discrete Fouler transform defined with

$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\displaystyle (1)~~~B_{k}=\sum _{m=0}^{N-1}A_{m}\exp {\Big (}-\mathrm {i} {\frac {2\pi }{N}}~k~m{\Big )}$

where $N$ is natural number; usually, $N=2^{q}$ for some natural $q$;
$A$ is the input array of length $N$, numbered beginning from 0;
$B$ is the output array of length $N$, numbered beginning from 0.
$k$ is also integer number, $0\leq k<N$.

There are many ways of evaluation of the expression (1). The trivial one would requires of order of $N^{2}$ operations. The efficient way requires of order on $N~\log _{2}N$ operations and called FFT. At large $N$, the time of evaluation by the direct implementation is many orders of magnitude larger than that with FFT; therefore, FFT is often used for the harmonic analysis, expansion by modes, calculation of the Fourier coefficients or approximation of the Fourier operator. There are several efficient algorithms to perform the FFT.

One of the simplest algorithms realized in the implementation of the C++ routine fft(*complex(double), int, int) is stored in file fafo.cin. This routine is used in the implementation fafo(*complex(double), int, int) of the Fourier operator.

More advanced algorithms allow the efficient evaluation even if $N$ is not an integer power of 2, although for the fast evaluation $N$ should not have large prime factors.

Common confusion

The often confusion about FFT is the believe that FFT somehow implements, or approximates the Fourier operator.

Unfortunately, this believe is just wrong. Actually, FFT does not approximate the Fourier operator at all.

However FFT can be used to make the numerical implementation of the Fourier operator. One of such implementations is loaded as fafo.cin.

The difference between FFT and the numerical implementation of the Fourier operator is discussed in the next section.

FFT and fafo

The routine "fafo" is used in calculations of field in lasers, described in articles and for plotting of some pictures in TORI. Before the similar method is used by E.A.Sziklas and A.E.Siegman and other researchers. Application of Fourier methods in optics is discussed also in the book by Okan Ersov .

Since year 2011, fafo.cin is available with examples. This routine takes less than one screen of code, and it seems to be compatible with various operational systems, the only C++ is required. The reason of use of the fafo.cin is not the performance (perhaps some factor of order of 0.5 can be achieved at the optimization of the CPU time), but its compactness, readability and portability; no any "installing" is necessary after the downloading.

Comparison of the self–Fourier Gaussian exponential (dotted curve) to its discrete representation $A$ by (2), blue line (and it coincides with its Fourier transform), and its FFT by (1), red line

The operator defined with (1) is often confused with the implementation of the Fourier operator; both replace the input array to the output one. But the result is not the same. The simplest way to see the difference is to perform the transform of the array of values of the Gaussian exponential, let

$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\displaystyle (2)~~~A_{k}=\exp(-x_{k}^{2}/2~)$

$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\displaystyle (3)~~~x_{k}=(k-N/2)d$

where $d={\sqrt {2\pi /N}}$ is step of the equidistant grid. It is assumed that $k$ is integer; $0\!\leq \!k\!<\!N$. The expression (2) defines the discrete analogy of the simplest self-Fourier function. With such an input array, The fafo returns almost the same array as that at the input (at least, beginning with $N=8$); the FFT returns some set of values that are difficult to interpret in terms of the Fourier operator.

The case of the Gaussian self-Fourier function by (2) is shown at the figure at right for $N\!=\!16$.
The dotted curve shows $y=\exp(-x^{2}/2)$; this is initial exponential.
The blue segmented line shows its discrete representation by (2). This representation practically coincide with the numerical estimate for its Fourier transform by routine fafo.
The Red curve shows the Discrete Fourier Transform, that is equivalent to equation (1).

If for the Gaussian exponential, the routine fft returns the array similar to that shown at the figure with red line, then this indicates that, perhaps, the routine works correctly and returns the approximation for the sum by equation (1).

The fft is used for the implementation of fafo; and it can be easy extracted from there. The algorithm of expression of fafo through the FFT is realized through the displacement of the elements of array. This can be done by various ways; the most straightforward (as simplest to understand) is realized.

Matlab and other languages

For the translation of fafo to other languages (for example, Matlab), where FFT is already implemented, there is no need to translate FFT routine. In particular, the conventional Matlabs's fft can be used "as is".

The numeric implementation of the fafo through the matlab built-in function fft can be even simplified using the built-in functions "fftshift" and "ifftshift"; For the array X, the correct call may have form

ifftshift(fftshift(X))

The way of realization of these functions may be dependent of the matlab version; one can play with these functions and experimentally find the way of calling, that approximates the Fourier operator

$\!\!\!\!\!\!\!\!\!\!\!\!\!(4)~~~\displaystyle \mathrm {Fourier} (A)(x)={\frac {1}{\sqrt {2\pi }}}\int _{-\infty }^{\infty }\exp(-\mathrm {i} xy)A(y)~\mathrm {d} y$

Keywords

Fourier operator, Fourier transform

References

The content of this article is adopted from http://tori.ils.uec.ac.jp/TORI/index.php/FFT