Efficient FFT Computation of Real Input

From Texas Instruments Embedded Processors Wiki

Introduction

Algorithms to perform forward and inverse Fast Fourier Transforms (FFT and IFFT) typically assume complex input and output data. However, many applications use only real-valued data in the time domain. A simple but inefficient solution to this problem is to pad N-length real input signals to N-length complex input signals with zero-valued imaginary components.

xreal = { 1, 2, 3, ... }
xcplx = { 1, 0, 2, 0, 3, 0, ... }

The complex FFT can then be applied to the double-length sequence. However, this method is obviously inefficient. This topic explains how to use the complex-valued FFT and IFFT algorithms to efficiently process real-valued sequences without padding those sequences. There are two key advantages to this approach:

1. Lower memory footprint - M bytes per sequence instead of 2M bytes
2. Fewer cycles - length N/2 FFT and IFFT computation instead of length N

We will jump straight to the algorithm and how to use it. For detailed information on the derivation of this method, please refer to the full application report on ti.com.

Example C code for this method is also available for download.

Prerequisites

• Install the C674x DSPLIB
  • Familiarize yourself with the FFT and IFFT kernels
  • Look at demo code (*_d.c) for usage examples

Computing a Length N/2 Complex FFT from a Length N Real Input Sequence

Let g(n) be an N-point real sequence. We want to compute the N-point complex FFT, but we only want to use an N/2-point FFT computation. We can accomplish this using the following steps.

1. Form the the N/2-point complex valued sequence, x(n) = x₁(n) + jx₂(n), where x₁(n) = g(2n) and x₂(n) = g(2n + 1).
2. Compute the N/2-point complex FFT on the complex valued sequence x(n) to obtain X(k) = FFT{x(n)}. Note that the FFT should be performed with bit reversal.
3. Perform an additional computation to get G(k) from X(k)

X(k)Gr(k) = Xr(k)Ar(k) – Xi(k)Ai(k) + Xr(N–k)Br(k) + Xi(N–k)Bi(k), for k = 0, 1, ..., N/2–1 and X(N/2) = X(0)

Gi(k) = Xi(k)Ar(k) + Xr(k)Ai(k) + Xr(N–k)Bi(k) – Xi(N–k)Br(k) 

Note that only N/2 points of the N-point sequence of G(k) are computed in the above equations. Because the FFT of a real-sequence has symmetric properties, we can easily compute the remaining N/2 points of G(k) with the following equations.

Gr(N/2) = Xr(0) – Xi(0) 
Gi(N/2) = 0 
Gr(N–k) = Gr(k), for k = 1, 2, ..., N/2–1 
Gi(N–k) = –Gi(k) 

As you can see, the above equations require A(k) and B(k), which are sine and cosine coefficients. These values can be precomputed with the following C code.

for (i = 0; i < N/2; i++)
{ 
  A[2 * i]     = 0.5 * (1.0 - sin (2 * PI / (double) (2 * n) * (double) i));
  A[2 * i + 1] = 0.5 * (-1.0 * cos (2 * PI / (double) (2 * n) * (double) i));
  B[2 * i]     = 0.5 * (1.0 + sin (2 * PI / (double) (2 * n) * (double) i)); 
  B[2 * i + 1] = 0.5 * (1.0 * cos (2 * PI / (double) (2 * n) * (double) i));
}

Returning to a Length N Real Sequence Using a Length N/2 Complex IFFT

Let G(k) be a N-point complex valued sequence derived from a real-valued sequence g(n). We want to get back g(n) = IFFT{G(k)}. However, we want to do this with an N/2-point IFFT. This can be accomplished using the following procedure.

Xr(k) = Gr(k)IAr(k) – Gi(k)IAi(k) + Gr(N–k)IBr(k) + Gi(N–k)IBi(k), for k = 0, 1, ..., N/2–1 and G(N/2) = G(0)
Xi(k) = Gi(k)IAr(k) + Gr(k)IAi(k) + Gr(N–k)IBi(k) – Gi(N–k)IBr(k)

1. Find X(k) using the above equations.
2. Compute the N/2-point inverse FFT of X(k) to get x(n) = x₁(n) + jx₂(n). Note that the IFFT should be performed with bit reversal.
3. Get g(n) from x(n) using the following equations.

g(2n) = x1(n), for n = 0, 1, ..., N/2–1
g(2n+1) = x2(n)

The above equations look similar to those used in our forward FFT computation. However, the pre-computed coefficients are slightly different. The following C code can be used to initialize IA(k) and IB(k).

for (i = 0; i < N/2; i++)
{ 
  IA[2 * i]     = 0.5 * (1.0 - sin (2 * PI / (double) (2 * n) * (double) i));
  IA[2 * i + 1] = 0.5 * (1.0 * cos (2 * PI / (double) (2 * n) * (double) i));
  IB[2 * i]     = 0.5 * (1.0 + sin (2 * PI / (double) (2 * n) * (double) i)); 
  IB[2 * i + 1] = 0.5 * (-1.0 * cos (2 * PI / (double) (2 * n) * (double) i));
}

Note that IA(k) is the complex conjugate of A(k) and IB(k) is the complex conjugate of B(k).

BenchMark

This table compares the performance of two methods of calculating the FFT of an N-point real sequence. Complex FFT refers to using an N-point complex FFT in the standard way, while Real FFT refers to using an N/2-point complex FFT as described in this topic. As expected, the Real FFT method yields superior performance.