Implementation of Fast Fourier Transform

Implementation of Fast Fourier Transform

Before reading the implementation details of FFT, please read Fast Fourier Transform document first.

Let us remember the FFT formula:

Fourier Formula

with k = 0,1,2,cdots,N/2-1 and $W=e^{-i 2pi/N}$

Now, let us choose N=8=2^3 as our specific example to understand easily. It can be generalized for larger afterwards. That is,
X sequence

We will examine the calculations step by step.

For a moment, just concentrate on the even and odd partitions and ignore the full formula. In the first stage:

First Stage

And in second stage:

Second Stage

Finally, in the third stage:

Third Stage

Another representation of this partition is shown as a tree:

Now, the question is Is there a pattern to this shuffling? Lets look at the binary representations which will give us a clue:

decimal	binary	operators
0	000	eee
4	100	oee
2	010	eoe
6	110	ooe
1	001	eeo
5	101	oeo
3	011	eoo
7	111	ooo

Since the sequences will include different number of elements in each step of the transform, we can modify the formula and replace by :

Modified Fourier Formula

with k = 0,1,2,cdots,n/2-1 and $W=e^{-i 2pi/N}$

There are two approaches to implement the algorithm:

Non-recursive Programming
- Without a look-up table
- With a look-up table
Recursive Programming

Nonrecursive Programming

Without a look-up table

Lets start by analysing the algorithm step by step when N=8 . The input is given as

and after bit reversal, we got

which is also eight 1-element DFTs.

We introduce a butterfly representation here which will make things clear.

Butterfly Representation

You see that the leftmost part of the butterfly is the deepest levels of the tree representation. We are doing the calculations on the original input (i.e: in place programming), there’s no memory requirement! The solid lines represent the added terms and the dashed lines reprenet the subtracted terms. Notice that we started with the bit reversed data.

stands for the stage number, n_1 is the distance between solid and dotted line pairs (i.e: distance between plus and minus terms of the same ), n_2 is the distance between same colors (i.e: distance between groups for that stage or having the same ).

The algorithm will work as follows:

Bit reverse data
for each stage ( $number of stages = log_{2}(N)$ )
- Calculate and
- for each color
  - Calculate
  - for each solid line of the selected color
    - Evaluate , then the new elements will be and .

All things considered, let us write the C code of the second part of this algorithm (the bit reversal is assumed to be done, i.e. the input is bit reversed):

n1 = 0;
n2 = 1;

We initialized the system.

for (i = 0; i < m; i++) {

for each stage

   n1 = n2;
   n2 = n2 + n2;

n_1 = n_2 of previous iteration
n_2 =n_2 + n_2 of previous iteration

  e = - twopi / n2;
  a = 0.0;

is initially 0.

  for(j = 0; j < n1; j++) {

for each color

    c = cos(a);
    s = sin(a);

Now, see that W =cos(a) + i sin(a) = c + i s by Euler expansion.

    for(k = j; k < N; k = k+n2)

for the selected color

      t1 = c*x[k+n1] - s*y[k+n1];
      t2 = s*x[k+n1] + c*y[k+n1];
 
      x[k+n1] = x[k] - t1;
      y[k+n1] = y[k] - t2;
 
      x[k] = x[k] + t1;
      y[k] = y[k] + t2;

Here, we used the multiplication of two complex numbers. Try to derive the temporary variables. Observe that the subtracted part is the dashed line and the added part is the solid line.

    }
    a = a + e;

order is as follows: $e^{i*0}, e^{i*e}, e^{i*2e}, cdots$ that is, cos(0) + i sin(0), cos(e) + i sin(e), cos(2e) + i sin(2e), cdots

}
}

With a look-up table

Observe that we are calculating many times in the FFT algorithm we introduced. But is there a way to eliminante these re-calculations?

Look at the butterfly representation above and see that we only use ‘s with k = 0,1, cdots, N/2-1 . Then why don’t we create a look-up table / array including these values initially. And then we can use these values directly.

Lets create the look-up table in C.

double *lookUpReal, *lookUpImag;

Since we are using C, we defined two pointers to the real and imaginary parts of the look-up table.

void createLookupTable(int size)
{
  double minusTwoPiOverN;
  double exponent;
  int j;
 
  lookUpReal = (double *) malloc (size * sizeof (double));
  lookUpImag = (double *) malloc (size * sizeof (double));
 
  minusTwoPiOverN = -6.283185307179586/(size+size);
  exponent = 0.0;
  for (j=0; j < size; j++) {
    lookUpReal[j] = cos(exponent);
    lookUpImag[j] = sin(exponent);
    exponent = exponent + minusTwoPiOverN;
  }
}

I didn’t put any explanation here, please refer to the previous section for the explanations of calculating s.

Then in the FFT code, we need to make some modifications, let us insert an initialization function call in the beginning:

createLookupTable(N/2);

Notice that, the size input will be half of the number of sequences.

n1 = 0;
n2 = 1;
                                            
for (i=0; i < L; i++) {
  n1 = n2;
  n2 = n2 + n2;
 
  for (j=0; j < n1; j++) {
    lookupIndex = j*(N/n2);
 
    WRe = lookUpReal[lookupIndex];  /* get the real part of the W factor */
    WIm = lookUpImag[lookupIndex];  /* get the imaginary part of the W factor */
 
    for (k=j; k < N; k=k+n2) {
      tempRe = WRe*xRe[k+n1] - WIm*xIm[k+n1];
      tempIm = WIm*xRe[k+n1] + WRe*xIm[k+n1];
      xRe[k+n1] = xRe[k] - tempRe;
      xIm[k+n1] = xIm[k] - tempIm;
      xRe[k]    = xRe[k] + tempRe;
      xIm[k]    = xIm[k] + tempIm;
    }
  }
}

All the calculations are similar to the previous one. Notice that lookupIndex values point to the index of the current in the look-up table.

Recursive Programming

The second approach is the recursive programming. Although it is very easy to code a recursive program, it is not preferably used technique for FFT, because it needs extra memory space and time for passing variables from one recursion node to the other. In-place programming in nonrecursive approach is very efficient in this manner.

Lets look at the FFT tree described in the previos section (above) again. The leftmost node is the root of the tree and the 2 nodes nearest to it are its children. Then the next 4 ones are the children of these 2 nodes. It continues like this till the end of the tree. At the end there are leaves. Recall that the DFT of one element is again itself.

After observing this, we see that it can easily be coded using recursion. The general recursion function will be:

fft(X) = fft(evens(X)) + W*fft(odds(X))

where X denotes a list of elements.

And the finishing condition for this recursion will be that the length of x is equal to 1. When this equality is provided, than it should be

fft(X) = X

By this way the recursion is supposed to end if we start with a dyadic number N = 2^k

Now let us write this using C. Since we will need the even and odd operators, look at them first.

double* D_even(double x[], int N) {
	double *y;
	int i;
	y = (double *) malloc(N/2 * sizeof(double));
	for(i=0;i<N/2;i++)
		y[i] = x[2*i];
	return y;
}
 
double* D_odd(double x[], int N) {
	double *y;
	int i;
	y = (double *) malloc(N/2 * sizeof(double));
	for(i=0;i<N/2;i++)
		y[i] = x[2*i+1];
	return y;
}

These routines take an array and return even/odd elements of it.

double** fftRecursive(double xRe[], double xIm[], int N) {
	double **y;

Since we calculate the real and imaginary parts, we need to return 2 arrays for each of them. This can be done by returning a double pointer that points to these arrays. y will be this returned variable.

	y = (double **) malloc(2 * sizeof(double *));

y is a double pointer composed of 2 pointers

	y[0] = (double *) malloc(N * sizeof(double));   /* Real values */
	y[1] = (double *) malloc(N * sizeof(double));   /* Imaginary values */

Each of its pointer points to N elements

	if(N == 1) {
		y[0][0] =xRe[0];
		y[1][0] =xIm[0];
	}

If there is one element in the signal, i.e. we are on the leaves of our recursion tree, FT of an element is again itself

	else {
		double **y_e;   /* Double pointer for even indexed transformed y's */
		double **y_o;   /* Double pointer for odd indexed transformed y's */
		double minusTwoPiOverN;
		double exponent;
		int k;
		minusTwoPiOverN = -6.283185307179586/(N);
		exponent = 0.0;
 
		y_e = fftRecursive(D_even(xRe,N), D_even(xIm,N), N/2);
		y_o = fftRecursive(D_odd(xRe,N), D_odd(xIm,N), N/2);
 
		for(k=0;k<N/2;k++) {
			double c=cos(exponent); /* cosine part of W factor */
			double s=sin(exponent); /* sine part of W factor */
            
                      /* In the equations below, finding out the real and imaginary parts are found by
                          using the equations 
                          y(k) = y_e(k) + (c+s*i)*y_o(k)
                          y(k+N/2) = y_e(k) - (c+s*i)*y_o(k) 
                          Below is just the same representation if the real and imaginary parts 
                          are considered seperately */
            
                      /* Any elements in the upper half (0..n/2-1) of the sequence is composed
                          of addition of even and scaled odd element. */
	      		y[0][k] = y_e[0][k] + c*y_o[0][k] - s*y_o[1][k];
			y[1][k] = y_e[1][k] + c*y_o[1][k] + s*y_o[0][k];
            
                      /* Similarly, any element on the other half is composed of subtraction */
			y[0][k+N/2] = y_e[0][k] - c*y_o[0][k] + s*y_o[1][k];
			y[1][k+N/2] = y_e[1][k] - c*y_o[1][k] - s*y_o[0][k];
            
                      /* Update the exponent */
			exponent = exponent + minusTwoPiOverN;
		}
	}

If there are more than 1 elements, then divide the array into odds and evens. Take the FFTs for these parts, and use the addition of evens FFT with odds FFT multiplied by to merge them. The explanation for calculation of can be found in non-recursive programming described above.

	return y;
}

Return the found real and imaginary parts.

It can also be written in MATLAB. Now, let us assume that input is a dyadic number (remember to apply zero padding if it is not).

function y = myFFTRecursive( x )
 
% If there is one element in the signal, i.e. we are on the leaves of our recursion tree
if n==1
    y = x;  % DFT of an element is again itself
else
% if there are more than 1 elements in the array
 
    % We need the FFT of the even and odd indexed cells
    y_e = myFFTRecursive(D_even(x));    
    y_o = myFFTRecursive(D_odd(x)); 
    
    for k=1:n/2
        W_nk = exp(-i*2*pi*(k-1)/n);      % The surplus factor we got in odd sequence.
        y(k)     = y_e(k) + W_nk*y_o(k);    % any elements in the upper half (1..n/2)
                                            % of the sequence is composed
                                            % of addition of even and
                                            % scaled odd element.
        y(k+n/2) = y_e(k) - W_nk*y_o(k);    % similarly, any element on the
                                            % other half is composed of
                                            % subtraction
    end         
end

And below are the even and odd operators we use in fft.

function y_e = D_even( y )
for k=1:length(y)/2
    y_e(k) = y(2*k-1);
end
 
function y_o = D_odd( y )
for k=1:length(y)/2
    y_o(k) = y(2*k);
end

Source Codes

nonrecursive programming - C : fft routines and a main function to test them written in C.
nonrecursive programming - MATLAB Mex-file : the same C code supplied within a mex function in order to use in MATLAB.
recursive programming - MATLAB Mex-file and MATLAB M-file : fft routine in mex function (in C) and in m-file.
recursive programming - MATLAB M-file with info : the m-file for fft routine with counters and displays.

References

A. Sevgen, Phys494 Applied Fourier Analysis Course Notes, Fall 2006, Bogaziçi University.
J.W.Cooley and J.W.Tukey. An algorithm for the machine calculation of complex fourier series. Mathematics of Computation, 19(2):297–301, April 1965.
Douglas Jones. "decimation-in time (dit) radix-2 fft". World Wide Web, http://cnx.rice.edu/content/m12016/latest, 2004. The Connexions Project module m12016.