Using DSPLIB FFT Implementation for Real Input and Without Data Scaling

Using DSPLIB FFT Implementation for Real Input and Without Data Scaling

Trademarks

All other trademarks are the property of their respective owners.

1 Real Input 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:

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

For more information on the derivation of this method, see Implementing Fast Fourier Transform Algorithms of Real-Valued Sequences With the TMS320 DSP Platform.

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

1.1 Prerequisites

• Install the C674x DSPLIB or C66x DSPLIB based on the DSP being used.
  • Familiarize yourself with the FFT and IFFT kernels
  • Look at demo code (*_d.c) for usage examples

1.2 Computing a Length N/2 Complex FFT From a Length N Real Input Sequence

Let g(n) be an N-point real sequence (N must be even). Compute the N-point complex FFT, but only use an N/2-point FFT computation. This can be accomplished by 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).
Gr(k) = Xr(k)Ar(k) – Xi(k)Ai(k) + Xr(N/2–k)Br(k) + Xi(N/2–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) + Xrr(N/2–k)Bi(k) – Xi(N/2–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, you 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, where even indices contains the real part and odd indices contain the imaginary part.

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)); }

1.3 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). You want to get back g(n) = IFFT{G(k)}. However, you 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/2–k)IBr(k) + Gi(N/2–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/2–k)IBi(k) – Gi(N/2–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) = x1(n) + jx2(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), again, even indices contains the real part and odd indices contain the imaginary part.

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).