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

*x[n]*is an

*nth*element of input sampled signal in time domain,

*X[m]*is the

*mth*element of output sequence in the frequency domain,

*N*is a total number of sample points taken. There are mainly two methods for calculating the

*X[m]*, one is a naive method and another is FFT method. I will explain each method separately.

**The simplest naive approach:**

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

**Fast Fourier Transform approach:**

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

what does rev_bits function do??

link is broken for the source, can you please update the link. Thank you

reverse bit

I’m trying to run your code but I get an error in fft1 line 22.