Complex DFT and FFT Algorithm Implementation Using C++ Programming Language

All four types of Fourier Transform family can be carried out with either real number or complex number. In my previous post, I shared how to implement real DFT algorithm using C++. In this post, I will implement the complex number version of DFT algorithm using C++. After that, I will also implement the Fast Fourier Transform (FFT) algorithm. FFT is another method for calculating the DFT. While it produces the same result as the DFT algorithm, it is incredibly more efficient, often reducing the computation time by hundreds. The algorithm that I use in 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 complex DFT transforms two N point time domain signals into two N point frequency domain signals. The two time domain signals are called the real part and the imaginary part, just as are the frequency domain signals. You can calculate the real DFT using complex DFT by move the N point signal into the real part of the complex DFT's time domain, then set all of the samples in the imaginary part to zero. Samples from 0 to N/2 of complex DFT correspond to the real DFT's spectrum.

This is the equation for calculating the complex DFT:

The compact version of this equation is in the exponential form. Many textbook use this form. From this form, we can expand the equation using Euler's formula, then with complex number multiplication we get final equation for calculate the real and imaginary part. This final equation is needed for implementing forward complex DFT code. The are several notation for the equation:

With equation 1, we can implement the forward complex DFT code:

// Implement forward complex DFT equation (Eq. 1)
void cdft()
{
    float sr, si;

    // Zero REX[] and IMX[] so they can be used as accumulators
    for (int k = 0; k < N; k++)
    {
        REX[k] = 0;
        IMX[k] = 0;
    }

    // Loop for each value in frequency domain
    for (int k = 0; k < N; k++)
    {
        // Loop for each value in time domain
        for (int n = 0; n < N; n++)
        {
            // Correlate with the complex sinusoid, sr and si
            sr = cos(2*PI*k*n/N);
            si = -sin(2*PI*k*n/N);
            REX[k] += xr[n]*sr - xi[n]*si;
            IMX[k] += xr[n]*si + xi[n]*sr;
        }
    }
}

As in real DFT, to calculate the inverse complex DFT we need to scale the frequency domain signals to get the sinusoids amplitudes by this equation:

After that, we can calculate the inverse complex DFT using this equation:

The equation for calculating the inverse complex DFT is similar to the forward complex DFT, but no minus sign on the exponential form. From this form we can expand the equation to get the final equation for calculating real and imaginary part of time domain signal. By using the equation 2 and final equation 3, we can implement the inverse complex DFT code:

// Implement inverse complex DFT equation (Eq. 3)
void icdft()
{
    float sr, si;

    // Find the cosine and sine wave amplitudes (Eq. 2)
    for (int k = 0; k < N; k++)
    {
        REX[k] /= N;
        IMX[k] /= N;
    }

    // Zero xr[] and xi[] so it can be used as an accumulators
    for (int n = 0; n < N; n++)
    {
        xr[n] = 0;
        xi[n] = 0;
    }

    // Loop for each value in time domain
    for (int n = 0; n < N; n++)
    {
        // Loop for each value in frequency domain
        for (int k = 0; k < N; k++)
        {
            // Correlate with the complex sinusoid, sr and si
            sr = cos(2*PI*k*n/N);
            si = -sin(2*PI*k*n/N);
            xr[n] += REX[k]*sr + IMX[k]*si;
            xi[n] += -REX[k]*si + IMX[k]*sr;
        }
    }
}

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 complex DFT. The original signal and inverse complex DFT signal should be the same. Actually there are not identical because of the round off error.

Original Signal:

 

Inverse Complex DFT Signal:

 

The FFT algorithm is another method for calculating the DFT. This is the code for forward and inverse FFT. If you want to learn more detail about this algorithm, you can read from The Scientist and Engineer's Guide to Digital Signal Processing book.

// Upon entry, REX[] and IMX[] contain the real and imaginary parts
// of the input
// Upon return, REX[] and IMX[] contain the FFT output
void fft()
{
    int nm1 = N - 1;
    int nd2 = N / 2;
    int m = log10(N) / log10(2);
    int j = nd2;
    int k;
    int le, le2;
    float ur, ui, sr, si;
    int jm1;
    int ip;
    float tr, ti;

    // Bit reversal sorting
    for (int i = 1; i <= N-2; i++)
    {
        if (i >= j) goto a;
        tr = REX[j];
        ti = IMX[j];
        REX[j] = REX[i];
        IMX[j] = IMX[i];
        REX[i] = tr;
        IMX[i] = ti;
        a:
        k = nd2;
        b:
        if (k > j) goto c;
        j -= k;
        k /= 2;
        goto b;
        c:
        j += k;
    }

    // Loop each stage
    for (int l = 1; l <= m; l++)
    {
        le = pow(2, l);
        le2 = le / 2;
        ur = 1;
        ui = 0;
        // Calculate sine and cosine values
        sr = cos(PI/le2);
        si = -sin(PI/le2);
        // Loop for each sub DFT
        for (int j = 1; j <= le2; j++)
        {
            jm1 = j - 1;
            // Loop for each butterfly
            for (int i = jm1; i <= nm1; i += le)
            {
                ip = i + le2;
                tr = REX[ip]*ur - IMX[ip]*ui;
                ti = REX[ip]*ui + IMX[ip]*ur;
                REX[ip] = REX[i] - tr;
                IMX[ip] = IMX[i] - ti;
                REX[i] = REX[i] + tr;
                IMX[i] = IMX[i] + ti;
            }
            tr = ur;
            ur = tr*sr - ui*si;
            ui = tr*si + ui*sr;
        }
    }
}

// Upon entry, REX[] and IMX[] contain the real and imaginary parts
// of the complex frequency domain
// Upon return, REX[] and IMX[] contain the complex time domain
void ifft()
{
    // Change the sign of IMX[]
    for (int k = 0; k <= N-1; k++)
    {
        IMX[k] = -IMX[k];
    }

    // Calculate forward FFT
    fft();

    // Divide the time domain by N and change the sign of IMX[]
    for (int n = 0; n <= N-1; n++)
    {
        REX[n] /= N;
        IMX[n] /= -N;
    }
}

You can get the project file for this tutorial from this link:

Complex DFT, FFT.