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 & 1 \\ 1 & - 1 \\ \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & \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)\).

Similarly, when \(N\) is a multiple of \(3\), we have

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

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

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    // size of sub problem
 61    const size_t nitems_sub = nitems / 2;
 62    // divide-and-conquer
 63    dft(nitems_sub, stride * 2, trigs, xs + 0 * stride, ys + 0 * nitems_sub);
 64    dft(nitems_sub, stride * 2, trigs, xs + 1 * stride, ys + 1 * nitems_sub);
 65    // unify sub problem results
 66    for (size_t k = 0; k < nitems_sub; k++) {
 67      const trig_t * const trig_y1 = trigs + stride * 1 * k;
 68      const double complex c_y1 = trig_y1->cos - I * trig_y1->sin;
 69      double complex * const y0_ptr = ys + k + 0 * nitems_sub;
 70      double complex * const y1_ptr = ys + k + 1 * nitems_sub;
 71      const double complex y0 = *y0_ptr;
 72      const double complex y1 = *y1_ptr * c_y1;
 73      *y0_ptr = 1. * y0 + 1. * y1;
 74      *y1_ptr = 1. * y0 - 1. * y1;
 75    }
 76  } else if (0 == nitems % 3) {
 77    // size of sub problem
 78    const size_t nitems_sub = nitems / 3;
 79    // divide-and-conquer
 80    dft(nitems_sub, stride * 3, trigs, xs + 0 * stride, ys + 0 * nitems_sub);
 81    dft(nitems_sub, stride * 3, trigs, xs + 1 * stride, ys + 1 * nitems_sub);
 82    dft(nitems_sub, stride * 3, trigs, xs + 2 * stride, ys + 2 * nitems_sub);
 83    // unify sub problem results
 84    for (size_t k = 0; k < nitems_sub; k++) {
 85      const trig_t * const trig_y1 = trigs + stride * 1 * k;
 86      const trig_t * const trig_y2 = trigs + stride * 2 * k;
 87      const double complex c_y1 = trig_y1->cos - I * trig_y1->sin;
 88      const double complex c_y2 = trig_y2->cos - I * trig_y2->sin;
 89      double complex * const y0_ptr = ys + k + 0 * nitems_sub;
 90      double complex * const y1_ptr = ys + k + 1 * nitems_sub;
 91      double complex * const y2_ptr = ys + k + 2 * nitems_sub;
 92      const double complex y0 = *y0_ptr;
 93      const double complex y1 = *y1_ptr * c_y1;
 94      const double complex y2 = *y2_ptr * c_y2;
 95      const double complex m = - 0.5 - 0.8660254037844386 * I;
 96      const double complex p = - 0.5 + 0.8660254037844386 * I;
 97      *y0_ptr = y0 +     y1 +     y2;
 98      *y1_ptr = y0 + m * y1 + p * y2;
 99      *y2_ptr = y0 + p * y1 + m * y2;
100    }
101  } else {
102    // naive O(N^2) DFT
103    for (size_t k = 0; k < nitems; k++) {
104      double complex * const y = ys + k;
105      *y = 0. + I * 0.;
106      for (size_t n = 0; n < nitems; n++) {
107        const size_t index = n * k % nitems;
108        *y += xs[stride * n] * cexp(- 2. * pi * index * I / nitems);
109      }
110    }
111  }
112  return 0;
113}

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} 1 & 1 \\ 1 & - 1 \\ \end{pmatrix} \begin{pmatrix} \frac{1}{2} & 0 \\ 0 & \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} \]

Similarly, when \(N\) is a multiple of \(3\), we have

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

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

src/dft.c
125static int idft (
126    const size_t nitems,
127    const size_t stride,
128    const trig_t * const trigs,
129    const double complex * const xs,
130    double complex * const ys
131) {
132  if (1 == nitems) {
133    ys[0] = xs[0];
134  } else if (0 == nitems % 2) {
135    const size_t nitems_sub = nitems / 2;
136    // divide-and-conquer
137    idft(nitems_sub, stride * 2, trigs, xs + 0 * stride, ys + 0 * nitems_sub);
138    idft(nitems_sub, stride * 2, trigs, xs + 1 * stride, ys + 1 * nitems_sub);
139    // unify sub problem results
140    for (size_t k = 0; k < nitems_sub; k++) {
141      const trig_t * const trig_y1 = trigs + stride * 1 * k;
142      const double complex c_y1 = trig_y1->cos + I * trig_y1->sin;
143      double complex * const y0_ptr = ys + k + 0 * nitems_sub;
144      double complex * const y1_ptr = ys + k + 1 * nitems_sub;
145      const double complex y0 = *y0_ptr        / 2.;
146      const double complex y1 = *y1_ptr * c_y1 / 2.;
147      *y0_ptr = y0 + y1;
148      *y1_ptr = y0 - y1;
149    }
150  } else if (0 == nitems % 3) {
151    const size_t nitems_sub = nitems / 3;
152    // divide-and-conquer
153    idft(nitems_sub, stride * 3, trigs, xs + 0 * stride, ys + 0 * nitems_sub);
154    idft(nitems_sub, stride * 3, trigs, xs + 1 * stride, ys + 1 * nitems_sub);
155    idft(nitems_sub, stride * 3, trigs, xs + 2 * stride, ys + 2 * nitems_sub);
156    // unify sub problem results
157    for (size_t k = 0; k < nitems_sub; k++) {
158      const trig_t * const trig_y1 = trigs + stride * 1 * k;
159      const trig_t * const trig_y2 = trigs + stride * 2 * k;
160      const double complex c_y1 = trig_y1->cos + I * trig_y1->sin;
161      const double complex c_y2 = trig_y2->cos + I * trig_y2->sin;
162      double complex * const y0_ptr = ys + k + 0 * nitems_sub;
163      double complex * const y1_ptr = ys + k + 1 * nitems_sub;
164      double complex * const y2_ptr = ys + k + 2 * nitems_sub;
165      const double complex y0 = *y0_ptr        / 3.;
166      const double complex y1 = *y1_ptr * c_y1 / 3.;
167      const double complex y2 = *y2_ptr * c_y2 / 3.;
168      const double complex m = - 0.5 - 0.8660254037844386 * I;
169      const double complex p = - 0.5 + 0.8660254037844386 * I;
170      *y0_ptr = y0 +     y1 +     y2;
171      *y1_ptr = y0 + p * y1 + m * y2;
172      *y2_ptr = y0 + m * y1 + p * y2;
173    }
174  } else {
175    // naive O(N^2) iDFT
176    for (size_t n = 0; n < nitems; n++) {
177      double complex * const y = ys + n;
178      *y = 0. + I * 0.;
179      for (size_t k = 0; k < nitems; k++) {
180        const size_t index = n * k % nitems;
181        *y += xs[stride * k] * cexp(+ 2. * pi * index * I / nitems);
182      }
183      *y /= nitems;
184    }
185  }
186  return 0;
187}