Robin's Blog

Stuff that doesn't fit elsewhere

February 13, 2013 / by Robin Scheibler / 4 comments

Real FFT Algorithms

Practical information on basic algorithms might be sometimes challenging to find. In this article, I break down two fundamental algorithms to compute the discrete Fourier transform (DFT, inverse transform is iDFT) of real-valued data using fast Fourier transform algorithm (FFT/iFFT).

When dealing with real data, these algorithms allow to divide the complexity of whatever FFT you are using by a factor of two. This is sometimes not negligible.

Nevertheless, don’t try to implement these algorithms, except for educational purpose. If you need to use these algorithms in practice, very efficient implementations, such as FFTW, already exists. The exception is if you are dealing with an architecture where FFTW is not available, or when the open source license of FFTW is not suitable.

Notation

We use lower case and the index \(n\) for time-domain sequences whereas the corresponding upper case letter, along with the index \(k\) is used for the corresponding frequency-domain sequence. As an example

\( X[k] = \operatorname{DFT}_k^N\{x\} = \sum\limits_{n=0}^{N-1} x_n e^{-j2\pi\frac{nk}{N}}, \qquad k=0,\ldots,N-1, \)

is the definition of the length-\(N\) DFT of the \(N\)-point sequence \(x[n]\), \(n=0,\ldots,N-1\). Its inverse DFT is defined as

\( x[n] = \operatorname{iDFT}_n^N\{X\} = \sum\limits_{k=0}^{N-1} X_n e^{j2\pi\frac{nk}{N}}, \qquad n=0,\ldots,N-1. \)

The complex conjugate is denoted by a star \(^*\) so that the famous conjugate symmetry of the DFT of real sequences is written as

\(X[k] = X^*[N-k],\) for a length-\(N\) real DFT.

Two for the price of one

How to compute two real DFTs of length \(N\) with a single complex DFT of the same size.

\( x[n], y[n] \in \mathbb{R},\qquad n=0,\ldots,N-1 \)

FFT

  1. Create complex sequence:

    \( z[n] = x[n] + jy[n],\qquad n=0,\ldots,N-1 \) Since \(x[n]\) and \(y[n]\) are real sequences, we have \( x[n] = \frac{z[n] + z^*[n]}{2}, \)
    \( y[n] = -j\frac{z[n]-z^*[n]}{2}. \)

  2. Using the linearity of the DFT we get:

    \( X[k] = \frac{Z[k] + Z_c[k]}{2} \)
    \( Y[k] = -j\frac{Z[k]-Z_c[k]}{2} \) where \(k=0,\ldots,N-1\) and \(Z_c[k] = \mathrm{DFT}_{k}^{N}\{z^*\}\). Since \(Z_c[k] = Z^*[N-k]\) (see proposition in appendix), we finally have: \( X[k] = \frac{Z[k] + Z^*[N-k]}{2} \)
    \( Y[k] = -j\frac{Z[k]-Z^*[N-k]}{2} \)

Algo: Two-for-One FFT algorithm

Input \( x[n],y[n]\in\mathbb{R}\), \(n=0,\ldots,N-1\)
Output \( X[k] = \mathrm{DFT}_{k}^{N}\{x\}\), \(Y[k]=\mathrm{DFT}_{k}^{N}\{y\} \)

  1. \( z[n] = x[n] + jy[n] \)
  2. \( Z[k] = \mathrm{DFT}_{k}^{N}\{z\} \)
  3. \( X[k] = \frac{Z[k] + Z^*[N-k]}{2} \)
  4. \( Y[k] = -j \frac{Z[k] - Z^*[N-k]}{2} \)

iFFT

By the linearity of the FFT we also have \(Z[k] = X[k] + jY[k]\). Thus the same algorithm is applied to iFFT.

Algo: Two-for-One iFFT algorithm

Input \(X[n],Y[n]\in\mathbb{C}\), \(n=0,\ldots,N-1\), DFTs of real sequences. Output \(x[n] = \mathrm{iDFT}_{n}^{N}\{X\}\), \(y[n]=\mathrm{iDFT}_{n}^{N}\{Y\}\) with \(x[n], y[n]\in\mathbb{R}\)

  1. \(Z[k] = X[k] + jY[k]\)
  2. \(z[n] = \mathrm{iDFT}_{n}^{N}\{Z\}\)
  3. \(x[n] = \mathrm{Re}\{z[n]\}\)
  4. \(y[n] = \mathrm{Im}\{z[n]\}\)

A single real FFT

Using a one-stage decimation in time (DIT) and the method used in the previous sequence, we can compute the DFT of a single real sequence \(x[n]\in\mathbb{R}\), \(n=0,\ldots,N-1\), \(N\) even, using a complex DFT of half the length.

FFT

  1. Use DIT:

    \( X[k] = \sum\limits_{n=0}^{N-1} x[n] e^{-j2\pi\frac{kn}{N}} \)
    \( \ = \sum\limits_{n_1=0}^{N/2-1} x[2n_1] e^{-j2\pi\frac{k(2n_1)}{N}} + \sum\limits_{n_2=0}^{N/2-1} x[2n_2+1] e^{-j2\pi\frac{k(2n_2+1)}{N}} \)
    \( \ = \sum\limits_{n_1=0}^{N/2-1} x[2n_1] e^{-j2\pi\frac{kn_1}{N/2}} + e^{-j2\pi\frac{k}{N}}\sum\limits_{n_2=0}^{N/2-1} x[2n_2+1] e^{-j2\pi\frac{kn_2}{N/2}} \)
    \( \ = X_e[k] + X_o[k]\,e^{-j2\pi\frac{k}{N}} \) where \(X_e[k] = \mathrm{DFT}{k}^{N/2}\{x_e\}\), \(X_o[k] = \mathrm{DFT}{k}^{N/2}\{x_o\}\) with \(x_e[n] = x[2n]\) and \(x_o[n] = x[2n+1]\), respectively, and \(k=0,\ldots,N/2-1\).

  2. Now use the Two-for-One method to compute \(X_e[k]\) and \(X_o[k]\) using a \(N/2\)-points complex FFT.
  3. \(z[n] = x_e[n] + jx_o[n]\)
  4. \(Z[n] = \mathrm{DFT}_{k}^{N/2}\{z\}\)
  5. \(X_e[n] = \frac{Z[k] + Z^*[N/2-k]}{2}\)

  6. \(X_o[n] = -j \frac{Z[k] - Z^*[N/2-k]}{2}\)

  7. We can now compute \(X[k]\), \(k=0,\ldots,N-1\), using \(X[k] = X_e[k\bmod{N/2}] + X_o[k\bmod{N/2}] e^{-j2\pi\frac{k}{N}}\)

Algo: Single Real FFT

Input \(x[n]\in\mathbb{R}\), \(n=0,\ldots,N-1\)
Output \(X[k] = \mathrm{DFT}_{k}^{N}\{x\}\)

  1. \(x_e[n] = x[2n]\), \(n=0,\ldots,N/2-1\)
  2. \(x_o[n] = x[2n+1]\), \(n=0,\ldots,N/2-1\)
  3. \(z[n] = x_e[n] + jx_o[n]\)
  4. \(Z[k] = \mathrm{DFT}_{k}^{N/2}\{z\}\)
  5. \(X_e[k] = \frac{Z[k] + Z^*[N/2-k]}{2}\)
  6. \(X_o[k] = -j \frac{Z[k] - Z^*[N/2-k]}{2}\)
  7. \(X[k] = X_e[k] + X_o[k] e^{-j2\pi\frac{k}{N}}\), for \(k=0,\ldots,N/2-1\)
  8. \(X[k] = X_e[k] - X_o[k]\), for \(k=N/2\)
  9. \(X[k] = X^*[N-k]\), for \(k=N/2+1,\ldots,N-1\)

iFFT

  1. To recover the \(X_e[k]\) and \(X_o[k]\) sequence, we use:

    \( X[k] + X^*[N/2-k] \)
    \( \ = \left(X_e[k]+X_o[k]e^{-j2\pi\frac{k}{N}}\right) + \left(X_e^*[N/2-k]+X_o^*[N/2-k]e^{j2\pi\frac{k-N/2}{N}}\right) \)
    \( \ = 2X_e[k] + e^{-j2\pi\frac{k}{N}}\left(1 + e^{j\pi}\right)X_o[k] \)
    \( \ = 2X_e[k] \)

    and

    \( X[k] - X^*[N/2-k] \)
    \( \ = \left(X_e[k]+X_o[k]e^{-j2\pi\frac{k}{N}}\right) - \left(X_e^*[N/2-k]+X_o^*[N/2-k]e^{j2\pi\frac{k-N/2}{N}}\right) \)
    \( \ = e^{-j2\pi\frac{k}{N}}\left(1 - e^{j\pi}\right)X_o[k] \)
    \( \ = 2X_o[k] e^{-j2\pi\frac{k}{N}} \)

    for \(k=0,\ldots,N/2-1\). And thus:

    \( X_e[k] = \frac{1}{2}\left(X[k]+X^*[N/2-k]\right) \)
    \( X_o[k] = \frac{1}{2}\left(X[k]-X^*[N/2-k]\right)e^{j2\pi\frac{k}{N}} \)

    for \(k=0,\ldots,N/2-1\).

  2. Now we can use the iDFT of length \(N/2\) on the sequence \(Z[k] = X_e[k] + jX_o[k]\) and we obtain \(z[n] = x_e[n] + jx_o[n]\).

  3. We reconstitute the original sequence \(x[n]\) from the real and imaginary parts of \(z[n]\).

Algo: Single Real iFFT

Input \(X[k]\in\mathbb{C}\), \(k=0,\ldots,N/2\)
Output \(x[n] = \mathrm{iDFT}_{k}^{N}\{X\}\), with \(x[n]\in\mathbb{R}\)

  1. \(X_e[k] = \frac{1}{2}\left(X[k] + X^*[N/2-k]\right)\), \(k=0,\ldots,N/2-1\)
  2. \(X_o[k] = \frac{1}{2}\left(X[k] - X^*[N/2-k]\right)e^{j2\pi\frac{k}{N}}\), \(k=0,\ldots,N/2-1\)
  3. \(Z[k] = X_e[k] + jX_o[k]\), \(k=0,\ldots,N/2-1\)
  4. \(z[n] = \mathrm{iDFT}_{k}^{N/2}\{Z\}\)
  5. \(x[2n] = \mathrm{Re}\{z[n]\}\), \(n=0,\ldots,N/2-1\)
  6. \(x[2n+1] = \mathrm{Im}\{z[n]\}\), \(n=0,\ldots,N/2-1\)

Appendix

Proposition

Let \(z[n] \in \mathbb{C}\), \(n=0,\ldots,N-1\), be a complex sequence. Then the DFT of the complex conjugate sequence \(z^*[n]\) is

\( \mathrm{DFT}_{k}^{N}\{z^*\} = \left[\mathrm{DFT}_{N-k}^{N}\{z\}\right]^* \)

Proof

\( \mathrm{DFT}_{k}^{N}\{z^*\} = \sum\limits_{n=0}^{N-1}z^*[n] e^{-j2\pi \frac{nk}{N}} \)
\( \ = \left[ \sum\limits_{n=0}^{N-1}z[n]\, e^{j2\pi \frac{nk}{N}} \right]^* \)
\( \ = \left[ \sum\limits_{n=0}^{N-1}z[n]\, e^{-j2\pi \frac{nN}{N}}\, e^{j2\pi \frac{nk}{N}} \right]^* \)
\( \ = \left[ \sum\limits_{n=0}^{N-1}z[n]\, e^{-j2\pi \frac{n(N-k)}{N}} \right]^* \)
\( \ = \left[ \mathrm{DFT}_{N-k}^{N}\{z\} \right]^* \)

\(\square\)