Calculation of Discrete Fourier Transform(DFT) in C/C++ using Naive and Fast Fourier Transform (FFT) method

void DFT::idft1()
{
    double pi2 = 2.0 * M_PI;
    double angleTerm,cosineA,sineA;
    double invs = 1.0 / size;
    for(unsigned int y = 0;y < size;y++) {
        output_seq[y] = 0;
        for(unsigned int x = 0;x < size;x++) {
            angleTerm = pi2 * y * x * invs;
            cosineA = cos(angleTerm);
            sineA = sin(angleTerm);
            output_seq[y].real += input_seq[x].real * cosineA - input_seq[x].imag * sineA;
            output_seq[y].imag += input_seq[x].real * sineA + input_seq[x].imag * cosineA;
        }
        output_seq[y] *= invs;
    }
}

This is the fastest method of calculating DFT. Many algorithms are developed for calculating the DFT efficiently. Using FFT to calculate DFT reduces the complexity from O(N^2) to O(NlogN) which is great achievement and reduces complexity in greater amount for the large value of N. The code snippet below implements the FFT of  1-dimensional sequence

Complex* FFT::fft1(Complex *data, unsigned int size, int log2n,  bool inverse)
{
        double angle, wtmp, wpr, wpi, wr, wi;
        int n = 1, n2;
        double pi2 = M_PI * 2.0;
        double scale = 1.0/size;
        Complex tc; 
        for (int k = 0; k < log2n; ++k)
        {
            n2 = n;
            n <<= 1;
            angle = (inverse)?pi2/n:-pi2/n;
            wtmp=sin(0.5*angle);
            wpr = -2.0*wtmp*wtmp;
            wpi = sin(angle);
            wr = 1.0;
            wi = 0.0; 
 
            for (int m=0; m < n2; ++m) {
                for (unsigned int i=m; i < size; i+=n) {
                    int j = i+n2;
                    tc.real = wr * data[j].real - wi * data[j].imag;
                    tc.imag = wr * data[j].imag + wi * data[j].real;
                    data[j] = data[i] - tc;
                    data[i] += tc;
                }
                wr=(wtmp=wr)*wpr-wi*wpi+wr;
                wi=wi*wpr+wtmp*wpi+wi;
            }
        }
        if(inverse) {
            for(int i = 0;i < n;i++) {
                data[i] *= scale;
            }
        }
        return data;
}

The Complex is a user-defined data type, you can use C++ built-in complex data type under too. Here size denotes the number of points of sequence, if the input sequence is not equal to the power of 2 of some number then it is padded with 0.
The above example is only for the 1-dimensional sequence, but most of the practical systems are 2 dimensional. For example, an Image is a two-dimensional function f(x, y). So to calculate the Fourier transform of an image, we need to calculate 2 dimensional FFT. Due to the separability property of DFT, we can compute the FFT along one direction and then other direction separately. For example first performing along the row and then al the ng column. The code below illustrates the use of separablility.

Complex** FFT::fft2(Complex **data, int r, int c, bool inverse)
{
    Complex *row = new Complex[r];
    Complex *column = new Complex[c];
    int log2w = r >> 1;
    int log2h = c >> 1; 
    // Perform FFT on each row
    for (int y = 0; y < c; ++y) {
      for (int x = 0; x < r; ++x) {
        int rb = rev_bits(x, r);
        row[rb] = data[x][y];
      }
      row = fft1(row, r, log2w, inverse);
      for (int x = 0; x < r; ++x) {
        data[x][y] = row[x];
      }
    } 
 
    // Perform FFT on each column
    for (int x = 0;x < r; ++x) {
      for (int y = 0; y < c; ++y) {
            int rb = rev_bits(y, c);
            column[rb] = data[x][y];
      }
      column = fft1(column, c, log2h, inverse);
      for (int y = 0; y < c; ++y) {
        data[x][y] = column[y];
      }
    }
    return data;
}

Full source code is available on GitHub.

The sample output of above program for 2D sequence is given below

Note: The output is verified using Matlab’s function fft2(X)