Real-valued discrete Fourier transform

We consider the discrete Fourier transform of a sequence composed of \(N\) real numbers \(x_n\):

\[X_k = \dft{k}{N}{x_0}{x_1}{x_{N - 1}} = \sum_{n = 0}^{N - 1} x_n \twiddle{- 2 \pi}{n k}{N},\]

and its inverse transform:

\[x_n = \idft{n}{N}{X_0}{X_1}{X_{N - 1}} = \frac{1}{N} \sum_{k = 0}^{N - 1} X_k \twiddle{2 \pi}{n k}{N}.\]

Note that these relations satisfy the Hermitian symmetry \(X_{N - k} = X_k^*\).

Derivation
\[ \begin{align}\begin{aligned}X_{N - k} & = \sum_{n = 0}^{N - 1} x_n \twiddle{- 2 \pi}{n \left( N - k \right)}{N}\\& = \sum_{n = 0}^{N - 1} x_n \exp \left( - 2 \pi n I \right) \twiddle{2 \pi}{n k}{N}\\& = \sum_{n = 0}^{N - 1} x_n \twiddle{2 \pi}{n k}{N}\\& = \left[ \sum_{n = 0}^{N - 1} x_n \twiddle{- 2 \pi}{n k}{N} \right]^* \,\, \left( \because x_n^* = x_n \right)\\& = X_k^*.\end{aligned}\end{align} \]

This symmetry can facilitate the reduction of the computational load, along with the fact that \(x_n \in \mathbb{R}\).

Order of signals

We store \(X_k\): the output of the forward transform and the input of the backward transform in the following manner:

\[ \begin{align}\begin{aligned}x \left[ 0 \right] & \leftarrow \real{X_0} = X_0,\\x \left[ 1 \right] & \leftarrow \real{X_1},\\& \vdots\\x \left[ N / 2 - 1 \right] & \leftarrow \real{X_{N / 2 - 1}},\\x \left[ N / 2 \right] & \leftarrow \real{X_{N / 2}} = X_{N / 2},\\x \left[ N / 2 + 1 \right] & \leftarrow \imag{X_{N / 2 + 1}},\\x \left[ N / 2 + 2 \right] & \leftarrow \imag{X_{N / 2 + 2}},\\& \vdots\\x \left[ N - 2 \right] & \leftarrow \imag{X_{N - 2}},\\x \left[ N - 1 \right] & \leftarrow \imag{X_{N - 1}}.\end{aligned}\end{align} \]

Namely the real part \(\real{X_k}\) is stored first up to \(k = N / 2\), which is followed by the imaginary part. Recall the Hermitian symmetry: \(\imag{X_k} = - \imag{X_{N - k}}\), and accordingly \(\imag{X_0} = \imag{X_{N / 2}} = 0\).

Forward transform

We consider the forward transform:

\[X_k = \dft{k}{N}{x_0}{x_1}{x_{N - 1}} = \sum_{n = 0}^{N - 1} x_n \twiddle{- 2 \pi}{n k}{N}\]

for \(\seq{k}{0}{1}{N - 1}\). Although this can be computed as the discrete Fourier transform of a complex sequence having zero imaginary numbers, we aim at computing it with less computational effort here.

By utilizing the decimation in time (assuming \(N\) is a multiple of \(2\)), we have

\[ \begin{align}\begin{aligned}\dft{k}{N}{x_0}{x_1}{x_{N - 1}} & = \dft{k}{N / 2}{x_0}{x_2}{x_{N - 2}} + \twiddle{- 2 \pi}{k}{N} \dft{k}{N / 2}{x_1}{x_3}{x_{N - 1}},\\\dft{k + N / 2}{N}{x_0}{x_1}{x_{N - 1}} & = \dft{k}{N / 2}{x_0}{x_2}{x_{N - 2}} - \twiddle{- 2 \pi}{k}{N} \dft{k}{N / 2}{x_1}{x_3}{x_{N - 1}},\end{aligned}\end{align} \]

for \(\seq{k}{0}{1}{N / 2 - 1}\). Instead of taking this recurrence relation as it is, we introduce \(z_n\) defined as

\[z_n = x_{2 n} + I x_{2 n + 1}\]

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

\[ \begin{align}\begin{aligned}x_{2 n } &= + \frac{1}{2} z_n + \frac{1}{2} z_n^*,\\x_{2 n + 1} &= - \frac{I}{2} z_n + \frac{I}{2} z_n^*.\end{aligned}\end{align} \]

Their discrete Fourier transforms, due to the linearity, lead to

\[ \begin{align}\begin{aligned}\dft{k}{N / 2}{x_0}{x_2}{x_{N - 2}} & = + \frac{1}{2} \dft{k}{N / 2}{z_0}{z_1}{z_{N / 2 - 1}} + \frac{1}{2} \dft{k}{N / 2}{z_0^*}{z_1^*}{z_{N / 2 - 1}^*},\\\dft{k}{N / 2}{x_1}{x_3}{x_{N - 1}} & = - \frac{I}{2} \dft{k}{N / 2}{z_0}{z_1}{z_{N / 2 - 1}} + \frac{I}{2} \dft{k}{N / 2}{z_0^*}{z_1^*}{z_{N / 2 - 1}^*},\end{aligned}\end{align} \]

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

Now two discrete Fourier transforms are involved, and we focus on the latter:

\[\dft{k}{N / 2}{z_0^*}{z_1^*}{z_{N / 2 - 1}^*} = \sum_{n = 0}^{N / 2 - 1} z_n^* \twiddle{- 2 \pi}{n k}{N / 2}.\]

Since this is equal to

\[\left( \dft{N / 2 - k}{N / 2}{z_0}{z_1}{z_{N / 2 - 1}} \right)^*,\]
Derivation
\[ \begin{align}\begin{aligned}\dft{k}{N / 2}{z_0^*}{z_1^*}{z_{N / 2 - 1}^*} = & \sum_{n = 0}^{N / 2 - 1} z_n^* \twiddle{- 2 \pi}{n k}{N / 2}\\= & \sum_{n = 0}^{N / 2 - 1} z_n^* \left[ \twiddle{- 2 \pi}{n \left( -k \right)}{N / 2} \right]^*\\= & \sum_{n = 0}^{N / 2 - 1} z_n^* \left[ \twiddle{- 2 \pi}{n \left( -k \right)}{N / 2} \right]^* \left[ \twiddle{- 2 \pi}{n N / 2}{N / 2} \right]^* \,\, \left( \because \left[ \twiddle{- 2 \pi}{n N / 2}{N / 2} \right]^* \equiv 1 \right)\\= & \sum_{n = 0}^{N / 2 - 1} z_n^* \left[ \twiddle{- 2 \pi}{n \left( -k \right)}{N / 2} \twiddle{- 2 \pi}{n N / 2}{N / 2} \right]^*\\= & \sum_{n = 0}^{N / 2 - 1} z_n^* \left[ \twiddle{- 2 \pi}{n \left( N / 2 - k \right)}{N / 2} \right]^*\\= & \left[ \sum_{n = 0}^{N / 2 - 1} z_n \twiddle{- 2 \pi}{n \left( N / 2 - k \right)}{N / 2} \right]^* \,\, \left( \because \left( a b \right)^* \equiv a^* b^* \right)\\= & \left( \dft{N / 2 - k}{N / 2}{z_0}{z_1}{z_{N / 2 - 1}} \right)^*.\end{aligned}\end{align} \]

we notice that the latter relation can be obtained from the former; namely we need to evaluate only one discrete Fourier transform

\[\dft{N / 2 - k}{N / 2}{z_0}{z_1}{z_{N / 2 - 1}}\]

to find

\[ \begin{align}\begin{aligned}\dft{k}{N / 2}{x_0}{x_2}{x_{N - 2}} & = + \frac{1}{2} \dft{k}{N / 2}{z_0}{z_1}{z_{N / 2 - 1}} + \frac{1}{2} \left( \dft{N / 2 - k}{N / 2}{z_0}{z_1}{z_{N / 2 - 1}} \right)^*,\\\dft{k}{N / 2}{x_1}{x_3}{x_{N - 1}} & = - \frac{I}{2} \dft{k}{N / 2}{z_0}{z_1}{z_{N / 2 - 1}} + \frac{I}{2} \left( \dft{N / 2 - k}{N / 2}{z_0}{z_1}{z_{N / 2 - 1}} \right)^*,\end{aligned}\end{align} \]

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

For the corner case \(k = 0\), we see

\[ \begin{align}\begin{aligned}\dft{k = N / 2}{N / 2}{z_0}{z_1}{z_{N / 2 - 1}} & = \sum_{n = 0}^{N / 2 - 1} z_n \twiddle{- 2 \pi}{n N / 2}{N / 2}\\& = \sum_{n = 0}^{N / 2 - 1} z_n\\& = \dft{k = 0}{N / 2}{z_0}{z_1}{z_{N / 2 - 1}}.\end{aligned}\end{align} \]

In summary, the discrete Fourier transform of a real sequence is computed as follows.

First, we create a signal \(z_n\) composed of \(N / 2\) complex numbers:

\[z_n = x_{2 n} + I x_{2 n + 1},\]

where \(\seq{n}{0}{1}{N / 2 - 1}\). Practically, this can be achieved merely by casting the user input as a series of complex numbers; namely no data manipulations are needed.

Next, the discrete Fourier transform of \(z_n\) is evaluated to find

\[Z_k \equiv \dft{k}{N / 2}{z_0}{z_1}{z_{N / 2 - 1}}\]

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

src/rdft.c
171dft(nh, 1, trigs, xs, zs);
src/rdft.c
193const double complex z0 =
194  + 0.5 * zs[2 * k + 0]
195  + 0.5 * zs[2 * k + 1] * I;
src/rdft.c
197const double complex z1 =
198  + 0.5 * zs[2 * (nh - k) + 0]
199  + 0.5 * zs[2 * (nh - k) + 1] * I;

The result is used to compute

\[ \begin{align}\begin{aligned}X_k^e & \equiv \dft{k}{N / 2}{x_0}{x_2}{x_{N - 2}} = + \frac{1}{2} Z_k + \frac{1}{2} Z_{N / 2 - k}^*,\\X_k^o & \equiv \dft{k}{N / 2}{x_1}{x_3}{x_{N - 1}} = - \frac{I}{2} Z_k + \frac{I}{2} Z_{N / 2 - k}^*,\end{aligned}\end{align} \]

where \(\seq{k}{0}{1}{N / 2 - 1}\):

src/rdft.c
201const double complex xe = +     z0 +     conj(z1);
src/rdft.c
203const double complex xo = - I * z0 + I * conj(z1);

Finally they are adopted to calculate

\[ \begin{align}\begin{aligned}\dft{k}{N}{x_0}{x_1}{x_{N - 1}} & = X_k^e + \twiddle{- 2 \pi}{k}{N} X_k^o,\\\dft{k + N / 2}{N}{x_0}{x_1}{x_{N - 1}} & = X_k^e - \twiddle{- 2 \pi}{k}{N} X_k^o,\end{aligned}\end{align} \]

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

Due to the order of signals, it is worthwhile to decompose the real and the imaginary parts:

\[ \begin{align}\begin{aligned}\real{\dft{k}{N}{x_0}{x_1}{x_{N - 1}}} & = \real{X_k^e} + \left\{ + \real{X_k^o} \ctwiddle{- 2 \pi}{k}{N} - \imag{X_k^o} \stwiddle{- 2 \pi}{k}{N} \right\},\\\imag{\dft{k + N / 2}{N}{x_0}{x_1}{x_{N - 1}}} & = \imag{X_k^e} + \left\{ - \real{X_k^o} \stwiddle{- 2 \pi}{k}{N} - \imag{X_k^o} \ctwiddle{- 2 \pi}{k}{N} \right\}:\end{aligned}\end{align} \]
src/rdft.c
205// X_k^e + exp(- 2 pi k / nitems) X_k^o
206// X_k^e - exp(- 2 pi k / nitems) X_k^o
207const double complex za = xe + xo * twiddle;
208const double complex zb = xe - xo * twiddle;
209xs[k     ] = creal(za);
210xs[k + nh] = cimag(zb);

Note that

\[Z_{N / 2} = Z_0\]

is used to get

\[ \begin{align}\begin{aligned}X_0^e & = \real{Z_0},\\X_0^o & = \imag{Z_0},\end{aligned}\end{align} \]

and as a consequence:

\[ \begin{align}\begin{aligned}\dft{0}{N}{x_0}{x_1}{x_{N - 1}} & = X_0^e + X_0^o,\\\dft{N / 2}{N}{x_0}{x_1}{x_{N - 1}} & = X_0^e - X_0^o:\end{aligned}\end{align} \]
src/rdft.c
174const double xe = zs[0];
175const double xo = zs[1];
176xs[ 0] = xe + xo;
177xs[nh] = xe - xo;

The whole process is given below for completeness:

src/rdft.c
158int rdft_exec_f (
159    rdft_plan_t * const plan,
160    double * const xs
161) {
162  if (NULL == plan) {
163    puts("uninitialized plan is passed");
164    return 1;
165  }
166  const size_t nitems = plan->nitems;
167  const size_t nh = nitems / 2;
168  const trig_t * const trigs = plan->trigs;
169  double * const zs = plan->buf;
170  // compute complex dft to find Z_k | 1
171  dft(nh, 1, trigs, xs, zs);
172  // compute FFT of the original real signal, corner cases | 6
173  {
174    const double xe = zs[0];
175    const double xo = zs[1];
176    xs[ 0] = xe + xo;
177    xs[nh] = xe - xo;
178  }
179  // compute FFT of the original real signal, common cases
180  for (size_t k = 1; k < nh; k++) {
181    // trigonometric table stores
182    //   cos(+ 2 pi k / nitems)
183    //   +
184    //   sin(+ 2 pi k / nitems) I
185    // here we need
186    //   cos(- 2 pi k / nitems)
187    //   +
188    //   sin(- 2 pi k / nitems) I
189    const double complex twiddle =
190      + 1. * trigs[k].cos
191      - 1. * trigs[k].sin * I;
192    // compute 1/2 Z_k | 3
193    const double complex z0 =
194      + 0.5 * zs[2 * k + 0]
195      + 0.5 * zs[2 * k + 1] * I;
196    // compute 1/2 Z_{N / 2 - k} | 3
197    const double complex z1 =
198      + 0.5 * zs[2 * (nh - k) + 0]
199      + 0.5 * zs[2 * (nh - k) + 1] * I;
200    // compute X_k^e = + 1/2 Z_k + 1/2 Z_{N / 2 - k}^* | 1
201    const double complex xe = +     z0 +     conj(z1);
202    // compute X_k^o = - I/2 Z_k + I/2 Z_{N / 2 - k}^* | 1
203    const double complex xo = - I * z0 + I * conj(z1);
204    // real and imaginary parts are stored separately | 6
205    // X_k^e + exp(- 2 pi k / nitems) X_k^o
206    // X_k^e - exp(- 2 pi k / nitems) X_k^o
207    const double complex za = xe + xo * twiddle;
208    const double complex zb = xe - xo * twiddle;
209    xs[k     ] = creal(za);
210    xs[k + nh] = cimag(zb);
211  }
212  return 0;
213}

Backward transform

We consider the inverse transform:

\[x_n = \idft{n}{N}{X_0}{X_1}{X_{N - 1}} = \frac{1}{N} \sum_{k = 0}^{N - 1} X_k \twiddle{2 \pi}{n k}{N}\]

for \(\seq{n}{0}{1}{N - 1}\), where \(X_k\) should satisfy the Hermitian symmetry:

\[X_{N - k}^* = X_k.\]

Also it should be ordered in a specific manner.

The inverse transform is obtained essentially by backtracing the forward procedures. To start, we calculate \(X_k^e\) and \(X_k^o\) following

\[ \begin{align}\begin{aligned}X_k^e & = \frac{1}{2} \dft{k}{N}{x_0}{x_1}{x_{N - 1}} + \frac{1}{2} \dft{k + N / 2}{N}{x_0}{x_1}{x_{N - 1}},\\X_k^o \twiddle{- 2 \pi}{k}{N} & = \frac{1}{2} \dft{k}{N}{x_0}{x_1}{x_{N - 1}} - \frac{1}{2} \dft{k + N / 2}{N}{x_0}{x_1}{x_{N - 1}},\end{aligned}\end{align} \]

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

Note that, since \(X_k\) does not store the entire signal, we need to recover the series of complex numbers:

src/rdft.c
247const double complex x0 =
248  + 0.5 * xs[         k]
249  - 0.5 * xs[nitems - k] * I;
src/rdft.c
251const double complex x1 =
252  + 0.5 * xs[nh - k]
253  + 0.5 * xs[nh + k] * I;

They are combined to yield \(X_k^e\) and \(X_k^o\):

src/rdft.c
255const double complex xe = (x0 + x1);
src/rdft.c
257const double complex xo = (x0 - x1) * twiddle;

Using these values, \(Z_k\) is recovered following

\[Z_k = X_k^e + I X_k^o:\]
src/rdft.c
259// real and imaginary parts are stored alternately
260const double complex z = xe + xo * I;
261zs[2 * k + 0] = creal(z);
262zs[2 * k + 1] = cimag(z);

For the corner case \(k = 0\), we have

\[ \begin{align}\begin{aligned}X_0^e & = \frac{1}{2} X_0 + \frac{1}{2} X_{N / 2}^*, = \frac{1}{2} X_0 + \frac{1}{2} X_{N / 2},\\X_0^o & = \frac{1}{2} X_0 - \frac{1}{2} X_{N / 2}^* = \frac{1}{2} X_0 - \frac{1}{2} X_{N / 2},\end{aligned}\end{align} \]

and due to \(\imag{X_0} = \imag{X_{N / 2}} = 0\), we obtain

\[Z_0 = X_0^e + I X_0^o.\]
src/rdft.c
230const double x0 = 0.5 * xs[ 0];
231const double x1 = 0.5 * xs[nh];
232const double complex xe = x0 + x1;
233const double complex xo = x0 - x1;
234const double complex z = xe + xo * I;
235zs[0] = creal(z);
236zs[1] = cimag(z);

An inverse transform with respect to \(Z_k\) is performed to yield \(z_n\):

\[z_n = \idft{n}{N / 2}{Z_0}{Z_1}{Z_{N / 2 - 1}}:\]
src/rdft.c
265idft(nh, 1, trigs, zs, xs);

From the result, we obtain

\[ \begin{align}\begin{aligned}x_{2 n } &= \real{z_n},\\x_{2 n + 1} &= \imag{z_n},\end{aligned}\end{align} \]

for \(\seq{n}{0}{1}{N / 2 - 1}\). Since the output signal is already ordered in a desired manner, no additional manipulation is necessary here.

The whole process is given below for completeness:

src/rdft.c
216int rdft_exec_b (
217    rdft_plan_t * const plan,
218    double * const xs
219) {
220  if (NULL == plan) {
221    puts("uninitialized plan is passed");
222    return 1;
223  }
224  const size_t nitems = plan->nitems;
225  const size_t nh = nitems / 2;
226  const trig_t * const trigs = plan->trigs;
227  double * const zs = plan->buf;
228  // create a complex signal, edge cases | 9
229  {
230    const double x0 = 0.5 * xs[ 0];
231    const double x1 = 0.5 * xs[nh];
232    const double complex xe = x0 + x1;
233    const double complex xo = x0 - x1;
234    const double complex z = xe + xo * I;
235    zs[0] = creal(z);
236    zs[1] = cimag(z);
237  }
238  // create a complex signal, common cases
239  for (size_t k = 1; k < nh; k++) {
240    // trigonometric table
241    //   cos(+ 2 pi k / nitems)
242    //   sin(+ 2 pi k / nitems)
243    const double complex twiddle =
244      + 1. * trigs[k].cos
245      + 1. * trigs[k].sin * I;
246    // compute 1/2 X_k | 3
247    const double complex x0 =
248      + 0.5 * xs[         k]
249      - 0.5 * xs[nitems - k] * I;
250    // compute 1/2 X_{k + N / 2} | 3
251    const double complex x1 =
252      + 0.5 * xs[nh - k]
253      + 0.5 * xs[nh + k] * I;
254    // X_k^e = (+ 1/2 X_k + 1/2 X_{k + N / 2}) | 1
255    const double complex xe = (x0 + x1);
256    // X_k^o = (+ 1/2 X_k - 1/2 X_{k + N / 2}) * exp(arg) | 1
257    const double complex xo = (x0 - x1) * twiddle;
258    // Z_k = X_k^e + X_k^o * I | 4
259    // real and imaginary parts are stored alternately
260    const double complex z = xe + xo * I;
261    zs[2 * k + 0] = creal(z);
262    zs[2 * k + 1] = cimag(z);
263  }
264  // compute complex idft to find z_n | 1
265  idft(nh, 1, trigs, zs, xs);
266  return 0;
267}

Reference