Fast Fourier Transform

Discrete Fourier Transform (DFT)

Let $\mathbb{F}$ be a field and $\omega \in {\mathbb{F}}^{\times }$ a primitive $n$-th root of unity (so by definition the natural number $n$ divides the order of the multiplicative group ${\mathbb{F}}^{\times }$). The DFT matrix is defined by:

$$F_{\omega }=\left[\begin{array}{ccccc}1& 1& 1& \dots & 1\\ 1& \omega & {\omega }^2& \dots & {\omega }^{n-1}\\ 1& {\omega }^2& {\omega }^4& \dots & {\omega }^{2(n-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }^{n-1}& {\omega }^{2(n-1)}& \dots & {\omega }^{(n-1{)}^2}\end{array}\right]$$

and the inverse DFT matrix is given by $F_{{\omega }^{-1}}$.

The discrete Fourier transform of length $n$ is the $\mathbb{F}$-linear map:

$$\begin{array}{rl}{\text{DFT}}_{\omega }:{\mathbb{F}}^n& \to {\mathbb{F}}^n\\ \mathbf{x}& \mapsto \mathbf{y}=F_{\omega }\mathbf{x}=(\sum _{i=0}^{n-1}{x}_i\cdot {\omega }^{ij}{)}_{j\in [[0,n-1]]}\end{array}$$

which evaluates a polynomial from $\mathbb{F}[x]$ of degree strictly less than $n$ at the distinct points $1,\omega ,\dots ,{\omega }^{n-1}$. The group elements $1,\omega ,\dots ,{\omega }^{n-1}$ being distinct for the following reason: let $x\in {\mathbb{F}}^{\times }$ a group element of order $n$; for any $0<m<n,x^m\ne 1$, and $x^n=1$. If $x^j=x^k$ for $1\le j<k\le n$, then $x^{k-j}=1$ with $0<k-j<n$, a contradiction.

For any $0\le i,j<n$:

$$\begin{array}{rl}(F_{\omega }{F}_{{\omega }^{-1}}{)}_{i,j}& =\sum _{k=0}^{n-1}{\omega }^{jk}({\omega }^{-1}{)}^{ik}\\ & =\sum _{k=0}^{n-1}{\omega }^{(j-i)k}\\ & =0\text{ if }i\ne j\text{ or }n\text{ if }i=j.\end{array}$$

As in any field, $d\in \mathbb{N}\setminus \{0\}$, for a primitive $n$-th root of unity $\omega$, $1+{\omega }^d+\cdots +{\omega }^{d(n-1)}=(1-({\omega }^n{)}^d)/(1-{\omega }^d)=0$ since ${\omega }^n=1$, and ${\omega }^{-1}$ is also a primitive $n$-th root of unity.

Thus $F_{\omega }{F}_{{\omega }^{-1}}=n\times {\text{Id}}_n$, which in turn implies $({\text{DFT}}_{\omega }{)}^{-1}=n^{-1}{\text{DFT}}_{{\omega }^{-1}}$. The inverse discrete Fourier transform is thus defined if $n$ is invertible in $\mathbb{F}$ (the characteristic of $\mathbb{F}$ is strictly greater than $n$).

When $\mathbb{F}$ is a finite field this transform is also called Number Theoretic Transform (NTT).

Note that we can consider a variant of the DFT where the vector holds points of an elliptic curve group defined over a finite field $E({\mathbb{F}}_p)$, which we call EC-DFT. The addition is then the elliptic curve point addition and the multiplication with roots of unity is the elliptic curve scalar multiplication (the roots of unity being defined over a prime field can be mapped to scalars).

Convolution product with the DFT

One can check that the $n$-DFT is an homomorphism from $({\mathbb{F}}^n,\ast )$ to $({\mathbb{F}}^n,\odot )$ where $\ast$ is the convolution product ¹ and $\odot$ is the element-wise multiplication: ${\text{DFT}}_{\omega }(A\ast B)={\text{DFT}}_{\omega }(A)\odot {\text{DFT}}_{\omega }(B)$ from which $A\ast B={\text{DFT}}_{\omega }^{-1}({\text{DFT}}_{\omega }(A)\odot {\text{DFT}}_{\omega }(B))$ follows.

Aside:

One way to see it is to remark that the DFT for field elements is a ring isomorphism (with inverse DFT^-1) given by the Chinese Remainder Theorem (CRT)

${\mathbb{F}}^n\cong \mathbb{F}[x]/⟨x^n-1⟩\cong {\Pi }_{0\le i<n}\mathbb{F}[x]/⟨x-{\omega }^i⟩.$

The first isomorphism corresponds to the fact that $\mathbb{F}[x]/⟨x^n-1⟩$ is a ring (of polynomials in $\mathbb{F}[x]$ taken modulo $x^n-1$) and an $\mathbb{F}$-vector space of dimension $n$. The second one is the CRT factorization of $\mathbb{F}[x]/⟨x^n-1⟩$ as the ideals $⟨x-{\omega }^i⟩$ are pairwise coprime so their intersection is their product $⟨x^n-1⟩$. Indeed, the DFT maps $f\in \mathbb{F}[x{]}_{<n}$ to $(f({\omega }^i)=f(x)\mathrm{mod}(x-{\omega }^i){)}_{i\in [[0,n-1]]}$. The inverse map comes down to the Lagrange interpolation (based on partial fraction decomposition) or Bézout (based on Bézout’s identity, whose coefficients can be computed with the extended Euclidean algorithm).

The Fast Fourier Transform (FFT)

The FFT from Cooley and Tukey (and Gauss) splits the computation of a DFT of size $n=n_1{n}_2$ into the computations of $n_2$ DFTs of size $n_1$ (inner sum) and $n_1$ DFTs of size $n_2$ (outer sum), for $0\le k<n_1$ and $0\le l<n_2$:

$$\begin{array}{c}P({\omega }^{k+n_{1}l})=\sum _{i=0}^{n_2-1}(\sum _{j=0}^{n_1-1}{P}_{i+n_{2}j}({\omega }^{n_2}{)}^{jk})({\omega }^{n_1}{)}^{il}{\omega }^{ik}.\end{array}$$

What follows is a reproduction of the derivation of the above FFT expression from the DFT from the section 2.3 of An Approach to Low-power, High-performance, Fast Fourier Transform Processor Design, Bevan M. Baas, 1999. The input vector $\mathbf{x}$ and output vector $\mathbf{y}$ are reshaped to form $n_1\times n_2$ matrices $X$ and $Y$ as follows: for integers $A,B,C,D$, $0\le i,k<n_1$, $0\le j,l<n_2$, $X_{i,j}=x_{Ai+Bj\mathrm{mod}n}$ and $Y_{k,l}=y_{Ck+Dl\mathrm{mod}n}=P({\omega }^{Ck+Dl})$.

Plugging this reindexing into the DFT equation we obtain (where ${\omega }_k={\omega }^{n/k}$, and is a primitive $k$-th root of unity):

$$y_j=\sum _{i=0}^{n-1}{x}_i{\omega }_n^{ij}⟹Y_{k,l}=\sum _{i=0}^{n_1-1}\sum _{j=0}^{n_2-1}{X}_{i,j}{\omega }_n^{(Ai+Bj)(Ck+Dl)}.$$

With a column-wise reshape of $\mathbf{x}$ (ie. $A=n_2,B=1$) and row-wise reshape of $\mathbf{y}$ (ie. $C=1,D=n_1$), we find the Cooley-Tuckey FFT:

$$\begin{array}{rl}{Y}_{k,l}& =\sum _{i=0}^{n_1-1}\sum _{j=0}^{n_2-1}{X}_{i,j}{\omega }_n^{(n_{2}i+j)(k+n_{1}l)}\\ & =\sum _{i=0}^{n_1-1}\sum _{j=0}^{n_2-1}{X}_{i,j}{\omega }_n^{n_{2}ik+jk+j{n}_{1}l}\\ & =\sum _{i=0}^{n_1-1}\left(\right[\sum _{j=0}^{n_2-1}{X}_{i,j}{\omega }_{n_2}^{jl}\left]{\omega }_n^{jk}\right){\omega }_{n_1}^{ik}(\ast )\\ & =\sum _{j=0}^{n_2-1}\left(\sum _{i=0}^{n_1-1}{X}_{j,i}{\omega }_{n_1}^{ik}\right){\omega }_{n_2}^{jl}{\omega }_n^{jk}.\end{array}$$

Radix-2 decimation in time vs. decimation in frequency

The terminology comes from the Digital Signal Processing community where the DFT maps a signal (DFT input) to its frequencies (DFT output).

The radices are the prime factors from the decomposition of the DFT domain size $n$. Typically $n=2^k$ so the FFT is radix-2. We consider $n=2^k$.

Decimation in time (DIT)

From $(\ast )$, choosing $n_1=n/2$, $n_2=2$: the input sequence of the DFT $\mathbf{x}$ is decomposed into the even- and odd-indexed subsequences $X_{i,0}=x_{2i}$ and $X_{i,1}=x_{2i+1}$:

$$Y_{k,l}=P({\omega }^{k+(n/2)l})=\sum _{i=0}^{n/2-1}{x}_{2i}{\omega }_{n/2}^{ik}+{\omega }_n^{(n/2)l+k}\sum _{i=0}^{n/2-1}{x}_{2i+1}{\omega }_{n/2}^{ik}$$

which amounts to the “textbook FFT formula”, with a change of variable, and for $0\le k<n$:

$$P({\omega }^k)=\sum _{i=0}^{n/2-1}{x}_{2i}{\omega }_{n/2}^{ik}+{\omega }_n^k\sum _{i=0}^{n/2-1}{x}_{2i+1}{\omega }_{n/2}^{ik}.$$

Thus we obtain a recursive definition of the DFT: ${\text{DFT}}_{{\omega }_n}(\{x_k\})={\text{DFT}}_{{\omega }_{n/2}}(\{x_{2i}\})+{\omega }_n^k{\text{DFT}}_{{\omega }_{n/2}}(\{x_{2i+1}\})$.

Decimation in frequency (DIF)

From $(\ast )$, choosing $n_1=2$, $n_2=n/2$: the output sequence of the DFT $\mathbf{y}$ is decomposed into the even- and odd-indexed subsequences $Y_{0,l}=P({\omega }^{2l})$ and $Y_{1,l}=P({\omega }^{2l+1})$:

$$\begin{array}{rl}{Y}_{k,l}=P({\omega }^{k+2l})& =\sum _{i=0}^{n/2-1}{x}_i{\omega }_n^{i(2l+k)}+\sum _{i=0}^{n/2-1}{x}_{n/2+i}{\omega }_n^{(n/2+i)(2l+k)}\\ & =\sum _{i=0}^{n/2-1}(x_i+(-1{)}^k{x}_{n/2+i}){\omega }_n^{i(2l+k)}\\ & =\sum _{i=0}^{n/2-1}(x_i+(-1{)}^k{x}_{n/2+i}){\omega }_{n/2}^{il}{\omega }_n^{ik}.\end{array}$$

Hence the recursive definition of the DFT: $\{y_{2i}\}={\text{DFT}}_{{\omega }_{n/2}}(\{x_i+x_{n/2+i}{\}}_{i\in [[0,n/2-1]]})$ and $\{y_{2i+1}\}={\text{DFT}}_{{\omega }_{n/2}}(\{(x_i-x_{n/2+i}){\omega }_n^i{\}}_{i\in [[0,n/2-1]]})$.

There is a nice symmetry: in the DIT, the input sequence decomposed in even- and odd- indexed subsequences, and the output sequence is decomposed in top and bottom subsequences; and vice-versa for the DIF.

Prime factor algorithm (FFT variant)

This algorithm allows computing FFTs on domains of size $n$ whose factorization is included in the factorization of $|{\mathbb{F}}_r^{\times }|=r-1$. Let $\omega$ be a primitive $n$-th root of unity.

When $n_1$ and $n_2$ are coprime then the Chinese Remainder Theorem (CRT) allows to re-index the DFT of size $n$ in such a way the inner DFTs don’t have to be multiplied by the $n$ extra factors ${\omega }^{iu}$.

We start from the expression of the DFT for $k=0,\dots ,n-1$: $P({\omega }^k)={\sum }_{l=0}^{n-1}{P}_l{\omega }^{lk}.$ We re-index both input and output coefficients $l$ and $k$ thanks to the ring isomorphism ${\mathbb{Z}}_n\cong {\mathbb{Z}}_{n_1}\times {\mathbb{Z}}_{n_2}$ given by the CRT. It turns out that certain combinations of input and output mappings allow getting rid of the extra factors ${\omega }^{iu}$. Let’s take as input mapping the CRT forward mapping ${\mathbb{Z}}_n\to {\mathbb{Z}}_{n_1}\times {\mathbb{Z}}_{n_2}:l\mapsto (l\mathrm{mod}{n}_1,l\mathrm{mod}{n}_2)$. Hence, by Bézout for $i=l\mathrm{mod}{n}_1$, $j=l\mathrm{mod}{n}_2$, $t_1=(1/n_1\mathrm{mod}{n}_2)$ and $t_2=(1/n_2\mathrm{mod}{n}_1)$, we obtain $P_l=P_{i{n}_2{t}_2+j{n}_1{t}_1}$. As output mapping, we take the Good’s mapping ${\mathbb{Z}}_{n_1}\times {\mathbb{Z}}_{n_2}\to {\mathbb{Z}}_n:(i,j)\mapsto j{n}_1+i{n}_2\mathrm{mod}n$. We could also have permuted the input and output mappings. Applying the two mappings, we get

$$\begin{array}{rl}P({\omega }^{k{n}_1+l{n}_2})& =\sum _{i=0}^{n_1-1}\sum _{j=0}^{n_2-1}{P}_{i{n}_2{t}_2+j{n}_1{t}_1}{\omega }^{(i{n}_2{t}_2+j{n}_1{t}_1)(l{n}_2+k{n}_1)}\\ & =\sum _{i=0}^{n_1-1}\sum _{j=0}^{n_2-1}{P}_{i{n}_2{t}_2+j{n}_1{t}_1}{\omega }^{i{n}_2{t}_{2}l{n}_2}{\omega }^{i{n}_2{t}_{2}k{n}_1}{\omega }^{j{n}_1{t}_{1}l{n}_2}{\omega }^{j{n}_1{t}_{1}k{n}_1}\\ & =\sum _{i=0}^{n_1-1}\sum _{j=0}^{n_2-1}{P}_{i{n}_2{t}_2+j{n}_1{t}_1}{\omega }^{i{n}_2{t}_{2}l{n}_2}{\omega }^{j{n}_1{t}_{1}k{n}_1}\\ & =\sum _{i=0}^{n_1-1}\sum _{j=0}^{n_2-1}{P}_{i{n}_2{t}_2+j{n}_1{t}_1}{\omega }_{n_1}^{i{t}_{2}l{n}_2}{\omega }_{n_2}^{j{t}_{1}k{n}_1}\\ & =\sum _{i=0}^{n_1-1}\left(\sum _{j=0}^{n_2-1}{P}_{i{n}_2{t}_2+j{n}_1{t}_1}{\omega }_{n_2}^{jk}\right){\omega }_{n_1}^{il}.\end{array}$$

So in the same way as implied by the equation $(\ast )$, we compute $n_1$ FFTs of length $n_2$, transpose the $n_1\times n_2$ matrix, and compute $n_2$ FFTs of length $n_1$. The cost of the transposition is somewhat negligible compared to the FFTs since it only changes the data layout.

Complexity: assuming $n_2=2^k$ and $n_1=p$ for small prime $p$: $O(2^k(p\times k+p^2))$.

Footnotes

1. $(f\times g)(x)={\sum }_k{h}_k$ where $h_k={\sum }_j{a}_{k-j}{b}_j$. ↩