Tuesday, May 14, 2013

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

Discrete Fourier Transform has great importance on Digital Signal Processing (DSP). There are many situations where analyzing the signal in frequency domain is better than that in the time domain. The Fourier Transform actually converts the function in the time domain to frequency domain,  some processing is done in the frequency domain, and finally, inverse Fourier transforms converts the signal back into the time domain. The term discrete means the signal is not continuous rather it is sampled by some sampling frequency i.e. only some samples are taken in a certain interval (also called period). The sampling frequency depends upon the frequency of the original signal and this must satisfy Nyquist criteria. A simple example of Fourier transform is applying filters in the frequency domain of digital image processing. Before looking into the implementation of DFT, I recommend you to first read in detail about the Discrete Fourier Transform in Wikipedia. If you are already familiar with it, then you can see the implementation directly.  The 1 dimensional DFT can be calculated by using the following formula
dfte 

Where 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:
In this method, we simply implement the above formula as it is without any optimization. We apply the Euler’s theorem and divide the formula to calculate real and imaginary components. This is simple to implement but has more computational complexity. As the value of N increases the computational complexity raises exponentially. Implementation in C/C++ is shown below

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

dftoutput

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




2 comments:

  1. what does rev_bits function do??

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

    ReplyDelete