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 time domain. The Fourier Transform actually converts the function in time domain to frequency domain,  some processing is done in frequency domain, and finally inverse Fourier transform converts the signal back into 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 certain interval (also called period). The sampling frequency depends upon the frequency of original signal and this must satisfy Nyquist criteria. A simple example of Fourier transform is applying filters in 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 nth element of input sampled signal in time domain, X[m] is the mth element of output sequence in frequency domain, N is total number of sample points taken. There are mainly two methods for calculating the X[m], one is 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 used 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 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 row and then along 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;
}


You can download the source code from here.

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