FFT Algorithm

Explanations on the complex FFT algorithm

We explain the FFT algorithm on the cplx layout. The reim layout follows the same mathematical principle, but separates real and imaginary parts.

FFT is usually described by separating the polynomial into odd and even powers of X (e.g. $P(X)=P_\textrm{even}(X^2)+X.P_\textrm{odd}(X^2)$). If we have the $m/2$ evaluations of $P_{\textrm{even}}$ and $P_{\textrm{odd}}$, we deduce the $m$ evaluations of $P$, grouping each root $P(\omega)$ with its opposite $P(-\omega)$ by computing twiddle factors between $a=P_{\textrm{even}}(\omega^2)$, $b=P_{\textrm{even}}(\omega^2)$, and $\omega$. Namely, $P(\omega)=a+\omega.b$ and $P(-\omega)=a-\omega.b$.

The above presentation has $O(m.\log_2(m))$ complexity, but requires either to move data around (permute the coefficients of P to isolate the even and odd parts), or to input the coefficient spaces in reversed bit order.

For practical applications, we are in general interested in keeping the coefficients in natural order, but are on the other hand more relaxed on the order in which the evaluations are presented in the FFT vector, as long as FFT and iFFT are consistent. We therefore use a dual description of the same algorithm.

Description of the FFT algorithm

Let $m$ be a power of two and $$P(X)=c_0+c_1.X+\dots+c_{m-1}.X^{m-1} \mod X^m-\exp(2i\pi.\omega_0)$$

For $\ell\in[0,\log_2(m)]$, $k\in[0,2^\ell-1]$ and $j\in[0,m/2^\ell-1]$, we call $P_{\ell,j}(X)$ the sub-polynomial of degree $2^\ell$ that contains all coefficients of index $j$ modulo $m/2^\ell$: $$ P_{\ell,j}(X) = \sum_{n=0}^{2^\ell-1} c_{j+nm/2^{\ell}}.X^j,$$ and we call

$$\omega_{\ell,k}=\exp\left(2i\pi.\frac{\omega_0}{2^\ell}+\textrm{frb}_k\right)$$

The main idea is to compute all the $u_{\ell,j,k}=P_{\ell,j}(\omega_{\ell,k})$, by increasing value of $\ell$, using the layout indexing of the following figure, from left to right in the case of FFT, and right to left for iFFT:

There are many interesting properties we can derive from this figure:

• On the first layer $(\ell=0,k=0)$, $(u_{0,0,0},\dots,u_{0,m-1,0})$ happen to be the $m$ coefficients of $P$ in natural order, which are the inputs of the FFT.
• On the last layer $(\ell=\log_2(m),j=0)$, $(u_{\log_2(m),0,0},\dots,u_{\log_2(m),0,m-1})$ happen to be the $m$ evaluations of $P$ in fracrevbit order, which are the outputs of the FFT.
• There is a nice recurrence formula to compute layer $\ell+1$ out of the previous one: $(u_{\ell+1,2k,j},u_{\ell+1,2k+1,j}) = \textrm{twiddle}(u_{\ell,k,j},u_{\ell,k,j+m/ 2^{\ell+1}},\omega_{\ell+1,2k})$ This means that the total number of twiddles is $m\log_2(m)/2$ and that only the current layer needs to reside in memory.
• The twiddle formula allow to update the layer in place, and in a mostly sequential manner in memory.