Real DFT Algorithm Implementation Using C++ Programming Language

Fourier analysis is a family of mathematical techniques for decomposing signals into sinusoids. Discrete Fourier Transform (DFT) is one of the Fourier transform family. There are four types of Fourier transform:

• Fourier Transform: is used for signals that are continious and aperiodic.
• Fourier Series: is used for signals that are continious and periodic.
• Discrete Time Fourier Transform (DTFT): is used for signals that are discrete and aperiodic.
• Discrete Fourier Transform (DFT): is used for signals that are discrete and periodic.

The only type of Fourier transform that can be used in DSP is the DFT. Because digital computers can only work with information that is discrete and finite in length. Real DFT is a version of the DFT that uses real numbers to represent the input and output signals. While the DFT is uses complex numbers to represent the input and output signals. In this post I will focus on how to implement real DFT algorithm using C++. This post is based on The Scientist and Engineer's Guide to Digital Signal Processing book. This book is very good for beginner to learn DSP.

The real DFT changes N point input signals to N/2+1 point output signals. The input signals is time domain signals that will be decomposed. The output signals is in frequency domain that consist of amplitudes of the component sine and cosine waves. The process of calculating the frequency domain is called decomposition, analysis, forward DFT, or simply, DFT. The calculation of the time domain signals from frequency domain is called synthesis, or the inverse DFT. This is the equation for calculating real DFT:

The are several notation for this equation:

We can calculate forward DFT using equation 1. This is the code in C++ for calculate forward DFT:

// Implement forward real DFT equation (Eq. 1)
void rdft()
{
    // Zero REX[] and IMX[] so they can be used as accumulators
    for (int k = 0; k < N_FREQ; k++)
    {
        REX[k] = 0;
        IMX[k] = 0;
    }

    // Loop through each sample in the frequency domain
    for (int k = 0; k < N_FREQ; k++)
    {
        // Loop through each sample in the time domain
        for (int i = 0; i < N_TIME; i++)
        {
            REX[k] += x[i] * cos(2*PI*k*i/N_TIME);
            IMX[k] += -x[i] * sin(2*PI*k*i/N_TIME);
        }
    }
}

Before we calculate the inverse DFT from REX[k] and IMX[k] we need to scale the frequency domain signals to get the sinusoids amplitudes by equation 2. After that we can calculate the inverse DFT by equation 3. This is the code in C++ for calculate inverse DFT:

// Implement inverse real DFT equation (Eq. 3)
void irdft()
{
    // Find the cosine and sine wave amplitudes (Eq. 2)
    for (int k = 0; k < N_FREQ; k++)
    {
        REX[k] /= (N_TIME/2);
        IMX[k] /= (N_TIME/2);
    }
    REX[0] /= 2;
    REX[N_FREQ] /= 2;

    // Zero x[] so it can be used as an accumulator
    for (int i = 0; i < N_TIME; i++)
    {
        x[i] = 0;
    }

    // Loop through each sample in the time domain
    for (int i = 0; i < N_TIME; i++)
    {
        // Loop through each sample in the frequency domain
        for (int k = 0; k < N_FREQ; k++)
        {
            x[i] += REX[k] * cos(2*PI*k*i/N_TIME);
            x[i] += -IMX[k] * sin(2*PI*k*i/N_TIME);
        }
    }
}

To check if this code is operating properly we can use MATLAB to plot the original time domain signal and time domain signal from inverse DFT. The original signal and inverse DFT signal should be the same. Actually there are not identical because of the round off error.

Original Signal:

Inverse DFT Signal:

You can get the project file from here. I made this project using Qt Creator IDE.