Discrete Fourier transform

Forward transform

For a sequence of \(N\) complex numbers \(z_n\) with \(\seq{n}{0}{1}{N - 1}\), the discrete Fourier transform is defined as

\[Z_k \equiv \sum_{n = 0}^{N - 1} z_n \twiddle{- 2 \pi}{n k}{N} \equiv \dft{k}{N}{z_0}{z_1}{z_{N - 1}},\]

where \(\seq{k}{0}{1}{N - 1}\).

Assuming that \(N\) is a multiple of \(2\), we obtain the following relation which is known as the decimation in time:

\[\begin{split}\begin{pmatrix} Z_{k } \\ Z_{k + N / 2} \\ \end{pmatrix} = \begin{pmatrix} 1 & \twiddle{- 2 \pi}{k}{N} \\ 1 & - \twiddle{- 2 \pi}{k}{N} \\ \end{pmatrix} \begin{pmatrix} \dft{k}{N / 2}{z_0}{z_2}{z_{N - 2}} \\ \dft{k}{N / 2}{z_1}{z_3}{z_{N - 1}} \\ \end{pmatrix}\end{split}\]

for \(\seq{k}{0}{1}{N / 2 - 1}\).

Derivation

We decompose the sum into the even and odd components:

\[ \begin{align}\begin{aligned}Z_k = & \sum_{n = 0}^{N / 2 - 1} z_{2 n} \twiddle{- 2 \pi}{2 n k}{N} + \sum_{n = 0}^{N / 2 - 1} z_{2 n + 1} \twiddle{- 2 \pi}{\left( 2 n + 1 \right) k}{N}\\= & \sum_{n = 0}^{N / 2 - 1} z_{2 n} \twiddle{- 2 \pi}{n k}{N / 2} + \twiddle{- 2 \pi}{k}{N} \sum_{n = 0}^{N / 2 - 1} z_{2 n + 1} \twiddle{- 2 \pi}{n k}{N / 2}\\= & \dft{k}{N / 2}{z_0}{z_2}{z_{N - 2}} + \twiddle{- 2 \pi}{k}{N} \dft{k}{N / 2}{z_1}{z_3}{z_{N - 1}},\end{aligned}\end{align} \]

where \(\seq{k}{0}{1}{N - 1}\). Assigning \(k \leftarrow k + N / 2\) to this relation yields

\[ \begin{align}\begin{aligned}Z_{k + N / 2} = & \sum_{n = 0}^{N / 2 - 1} z_{2 n} \twiddle{- 2 \pi}{n \left( k + N / 2 \right)}{N / 2} + \twiddle{- 2 \pi}{\left( k + N / 2 \right)}{N} \sum_{n = 0}^{N / 2 - 1} z_{2 n + 1} \twiddle{- 2 \pi}{n \left( k + N / 2 \right)}{N / 2}\\= & \sum_{n = 0}^{N / 2 - 1} z_{2 n} \twiddle{- 2 \pi}{n k}{N / 2} - \twiddle{- 2 \pi}{k}{N} \sum_{n = 0}^{N / 2 - 1} z_{2 n + 1} \twiddle{- 2 \pi}{n k}{N / 2}\\= & \dft{k}{N / 2}{z_0}{z_2}{z_{N - 2}} - \twiddle{- 2 \pi}{k}{N} \dft{k}{N / 2}{z_1}{z_3}{z_{N - 1}},\end{aligned}\end{align} \]

indicating that considering \(\seq{k}{0}{1}{N / 2 - 1}\) is sufficient to obtain \(Z_k\) for \(\left( \seq{k}{0}{1}{N - 1} \right)\).

src/dft.c
50static int dft (
51    const size_t nitems,
52    const size_t stride,
53    const trig_t * const trigs,
54    const double complex * const xs,
55    double complex * const ys
56) {
57  if (1 == nitems) {
58    ys[0] = xs[0];
59  } else if (0 == nitems % 2) {
60    const size_t nh = nitems / 2;
61    // divide-and-conquer
62    dft(nh, stride * 2, trigs, xs         , ys     );
63    dft(nh, stride * 2, trigs, xs + stride, ys + nh);
64    // unify two sub problem results
65    for (size_t i = 0; i < nh; i++) {
66      const size_t j = i + nh;
67      const trig_t * const trig = trigs + stride * i;
68      const double complex twiddle = trig->cos - I * trig->sin;
69      double complex * const yi = ys + i;
70      double complex * const yj = ys + j;
71      const double complex ye = *yi;
72      const double complex yo = *yj * twiddle;
73      *yi = ye + yo;
74      *yj = ye - yo;
75    }
76  } else {
77    // naive O(N^2) DFT
78    for (size_t k = 0; k < nitems; k++) {
79      double complex * const y = ys + k;
80      *y = 0. + I * 0.;
81      for (size_t n = 0; n < nitems; n++) {
82        *y += xs[stride * n] * cexp(- 2. * pi * n * k * I / nitems);
83      }
84    }
85  }
86  return 0;
87}

Backward transform

The inverse transform is defined as an identical way with the opposite sign of the twiddle factor (and the pre-factor \(1 / N\)):

\[z_n \equiv \frac{1}{N} \sum_{k = 0}^{N - 1} Z_k \twiddle{2 \pi}{n k}{N} \equiv \idft{n}{N}{Z_0}{Z_1}{Z_{N - 1}}\]

for \(\seq{n}{0}{1}{N - 1}\), and similarly we can deduce

\[\begin{split}\begin{pmatrix} z_n \\ z_{n + N / 2} \\ \end{pmatrix} = \begin{pmatrix} \frac{1}{2} & \frac{1}{2} \twiddle{2 \pi}{n}{N} \\ \frac{1}{2} & - \frac{1}{2} \twiddle{2 \pi}{n}{N} \\ \end{pmatrix} \begin{pmatrix} \idft{n}{N / 2}{Z_0}{Z_2}{Z_{N - 2}} \\ \idft{n}{N / 2}{Z_1}{Z_3}{Z_{N - 1}} \\ \end{pmatrix},\end{split}\]

for \(\seq{n}{0}{1}{N / 2 - 1}\).

Derivation

The even-odd decomposition yields

\[ \begin{align}\begin{aligned}\idft{n}{N}{Z_0}{Z_1}{Z_{N - 1}} & = \frac{1}{N} \sum_{k = 0}^{N - 1} Z_k \twiddle{2 \pi}{n k}{N}\\& = \frac{1}{N} \sum_{k = 0}^{N / 2 - 1} Z_{2 k} \twiddle{2 \pi}{n \left( 2 k \right)}{N} + \frac{1}{N} \sum_{k = 0}^{N / 2 - 1} Z_{2 k + 1} \twiddle{2 \pi}{n \left( 2 k + 1 \right)}{N}\\& = \frac{1}{2} \frac{1}{N / 2} \sum_{k = 0}^{N / 2 - 1} Z_{2 k} \twiddle{2 \pi}{n k}{N / 2} + \frac{1}{2} \twiddle{2 \pi}{n}{N} \frac{1}{N / 2} \sum_{k = 0}^{N / 2 - 1} Z_{2 k + 1} \twiddle{2 \pi}{n k}{N / 2}\\& = \frac{1}{2} \idft{n}{N / 2}{Z_0}{Z_2}{Z_{N - 2}} + \frac{1}{2} \twiddle{2 \pi}{n}{N} \idft{n}{N / 2}{Z_1}{Z_3}{Z_{N - 1}}.\end{aligned}\end{align} \]

Assigning \(n \leftarrow n + N / 2\) yields

\[ \begin{align}\begin{aligned}\idft{n + N / 2}{N}{Z_0}{Z_1}{Z_{N - 1}} & = \frac{1}{2} \frac{1}{N / 2} \sum_{k = 0}^{N / 2 - 1} Z_{2 k} \twiddle{2 \pi}{\left( n + N / 2 \right) k}{N / 2} + \frac{1}{2} \twiddle{2 \pi}{n + N / 2}{N} \frac{1}{N / 2} \sum_{k = 0}^{N / 2 - 1} Z_{2 k + 1} \twiddle{2 \pi}{\left( n + N / 2 \right) k}{N / 2}\\& = \frac{1}{2} \frac{1}{N / 2} \sum_{k = 0}^{N / 2 - 1} Z_{2 k} \twiddle{2 \pi}{n k}{N / 2} \exp \left( 2 \pi k I \right) + \frac{1}{2} \twiddle{2 \pi}{n}{N} \exp \left( \pi I \right) \frac{1}{N / 2} \sum_{k = 0}^{N / 2 - 1} Z_{2 k + 1} \twiddle{2 \pi}{n k}{N / 2} \exp \left( 2 \pi k I \right)\\& = \frac{1}{2} \frac{1}{N / 2} \sum_{k = 0}^{N / 2 - 1} Z_{2 k} \twiddle{2 \pi}{n k}{N / 2} - \frac{1}{2} \twiddle{2 \pi}{n}{N} \frac{1}{N / 2} \sum_{k = 0}^{N / 2 - 1} Z_{2 k + 1} \twiddle{2 \pi}{n k}{N / 2}\\& = \frac{1}{2} \idft{n}{N / 2}{Z_0}{Z_2}{Z_{N - 2}} - \frac{1}{2} \twiddle{2 \pi}{n}{N} \idft{n}{N / 2}{Z_1}{Z_3}{Z_{N - 1}}.\end{aligned}\end{align} \]
src/dft.c
 99static int idft (
100    const size_t nitems,
101    const size_t stride,
102    const trig_t * const trigs,
103    const double complex * const xs,
104    double complex * const ys
105) {
106  if (1 == nitems) {
107    ys[0] = xs[0];
108  } else if (0 == nitems % 2) {
109    const size_t nh = nitems / 2;
110    // divide-and-conquer
111    idft(nh, stride * 2, trigs, xs         , ys     );
112    idft(nh, stride * 2, trigs, xs + stride, ys + nh);
113    // unify two sub problem results
114    for (size_t i = 0; i < nh; i++) {
115      const size_t j = i + nh;
116      const trig_t * const trig = trigs + stride * i;
117      const double complex twiddle = trig->cos + I * trig->sin;
118      double complex * const yi = ys + i;
119      double complex * const yj = ys + j;
120      const double complex ye = 0.5 * *yi;
121      const double complex yo = 0.5 * *yj * twiddle;
122      *yi = ye + yo;
123      *yj = ye - yo;
124    }
125  } else {
126    // naive O(N^2) iDFT
127    for (size_t n = 0; n < nitems; n++) {
128      double complex * const y = ys + n;
129      *y = 0. + I * 0.;
130      for (size_t k = 0; k < nitems; k++) {
131        *y += xs[stride * k] * cexp(+ 2. * pi * n * k * I / nitems);
132      }
133      *y /= nitems;
134    }
135  }
136  return 0;
137}