Fast Fourier Transform

Overview

This document introduces the techniques of fast Fourier transform.

Introduction to Fast Fourier Transform

The discrete Fourier transform is computationally intensive requiring N² complex multiplications for a set of N elements. This problem is exacerbated when working with 2-dimensional data like images. An image of size M x M will require (M²)² or M⁴ complex multiplications.

Fortunately, in 1942, it was discovered that the discrete Fourier transform of length N could be rewritten as the sum of two Fourier transforms of length N/2. This concept can be recursively applied to the data set until it is reduced to transforms of only two points. Due partially to the lack of computing power, it wasn't until the mid 1960s that this discovery was put into practical application. In 1965, J.W. Cooley and J.W. Tukey applied this finding at Bell Labs to filter noisy signals.

This divide and conquer technique is known as the fast Fourier transform. It reduces the number of complex multiplications from N² to the order of N log₂ N. Table 7.1 shows the computations and time required to perform the DFT directly and via the FFT. It is assumed that each complex multiply takes 1 microsecond.
This savings is substantial especially when image processing. The FPT is separable, which makes Fourier transforms even easier to do. Because of the separability, we can reduce the FFT operation from a 2-dimensional operation to two 1-dimensional operations. First we compute the FFT of the rows of an image and then follow up with the FFT of the columns. For an image of size M x N, this requires N + M FFTs to be computed. The order of NM log₂ NM computations are required to transform our image. Table 7.2 shows the computations and time required to perform the DFT directly and via the FFT.



There are some considerations to keep in mind when transforming data to the frequency domain via the FPT. First, since the FFT algorithm recursively divides the data down, the dimensions of the image must be powers of 2 (N = 2^j and M = 2^k where j and k can be any number). Chances are pretty good that your image dimensions are not a power of 2. Your image data set can be expanded to the next legal size by surrounding the image with zeros. This is called zero-padding. You could also scale the image up to the next legal size or cut the image down at the next valid size. Listing 7.1 shows how to implement the 2-dimensional FFT in software. This code calls the 1-dimensional FFT routine for each row" and column. The actual FFT function can be broken down into two main functions.

Listing 7.1 - Implementing 2-Dimensional FFT

/**********************************************************************
Func: twoD_FFT
                                                                                     
Desc: performs an FFT on image data

Params: complex_data - array of complex image data
rows - number of rows in source image
cols- number of cols in source image
dir - direction of FFT (l=forward -l=reverse)
**********************************************************************/

void twoD_FFT(complex_ptr complex_data, int rows, int cols, int dir)
{
unsigned long x,y;       /* x and y indices to source image */
unsigned long index;     /* index to output line buffer */
COMPLEX *col_data;       /* storage for the columns */
COMPLEX *row_data;       /* storage for the rows */  
int M, N;                /* power of 2 for each dimension */
int num;                 /* temporary variable */

/* compute power of 2s */
num = cols;
M = 0;
while(num >= 2)
{
num »= 1;
M++;          
}
num = rows;
N = 0;                 
while(num >= 2)
{                                              
num »= 1;
N++;            
}
               
/* allocate memory for storage */
col_data = (COMPLEX *) malloc(sizeof(COMPLEX)*rows);
row_data = (COMPLEX *) malloc(sizeof(COMPLEX)*cols);
if((col_data == NULL) || (row_data == NULL))
exit(l);
for(y=0; y<rows; y++)     /* FFT on all rows */
{      
printf("Processing row %lu of %d\n",y+l, rows);    
index = y * cols;                              
for(x=0; x<cols; x++)     /* copy row data into array */
{
row_data[x] .re = complex_data [ index].re;
row_data[x].im = complex_data[index++].im;
}
fft(row_data, ,M, cols, dir);/* compute forward FFT */
index = y * cols;
for(x=0; x<cols; x++)        /* copy row back */
{
complex_data[index].re = row_data[x].re;
complex_data[index++].im = row_data[x].im;
   }
}

for(x=0; x<cols; x++)
{              
printf("Processing column %lu of %d \n", x+1, cols);
index = x;
for(y=0; y<rdws; y++)        /* copy cols into array */
{
col_data[y].re = complex_data[index].re;
col_data[y].im = complex_data[index].im;
index += cols;
}  

fft(col_data, N, rows, dir); /* perform forward FFT */ index = x;
index = x;
for(y=0; y<rows; y++)
{
complex_data[index].re = col_data[y].re;
complex_data [ index] . im = col_data[y].im;
index += cols;
}
}
free (col_data);
}
               
/**********************************************************************
Func: fft

Desc: performs forward and reverse Fast Fourier Transform

Params: f- complex array of data
 logN- power of 2 (numpoints = 2^logN)
 numpoints- number of elements in the data array
 dir- direction of fft (l=forward, -l=reverse)  
**********************************************************************/

void fft(COMPLEX *f, int logN, int numpoints, int dir) '  
{
scramble(numpoints,,f) ;
butterflies(numpoints, logM, dir, f);
}  

/**********************************************************************
Func:: scramble
Desc: Scrambles the input data by swapping data with data in
element with bit reversed index  
Params: numpoints- number of values in the complex array
f - complex array of data
**********************************************************************/
   
void scramble(int numpoints, COMPLEX *f)
{
int i, j, m;   /* loop variables */
double temp;   /* temporary storage during swapping */

J = 0;    
for (i=0;i<numpoints;i++)                
{

if (i > j)
{               /* swap f[j] and f[i] */
temp = f[j] .re;    /* swap .real */
f[j] -re = f[i] .re;  
f[i].re = temp;
temp = f[j].im;     /* swap imaginary */
f[j].im = f[i] .im;
f[i].im = temp;
}
m = numpoints»l;  
while ((j >= m) & (m >= 2))
{
J -= m;
m = m » 1;
}
J += m;
}
 }

/**********************************************************************
Func: butterflies
Desc: Calculates the butterflies for data in fft
if  a reverse fft, points are divided by numpoints
Params: numpoints- number of points in the complex array
logN- power of 2 (numpoints = 2^logN)    
dir- direction of fft (l=forward, -l=reverse) k        
f- array of data
**********************************************************************/

void butterflies(int numpoints, int logN, int dir, COMPLEX *f)
{      
double angle;              /* polar angle */
COMPLEX w, wp, temp;       /* intermediat complex numbers */
int i, j, k, offset;       /* loop variables */
 int N, half_M;            /* storage for powers of 2 */
double wtemp;              /* temporary storage for w.re */

N = 1;
for(k=0; k<logN; k++)
{
half_N = M;
N «= 1;                   /* multiply N by 2 */
angle = -2.0 * PI / N * dir;
wtemp = sin(0.5 * angle);
wp.re = -2.0 * wtemp * wtemp;

wp.im = sin(angle);
w.re = 1.0;
w.im = 0 . 0 ;
for(offset=0; offset<half_N; offset++)
{
for(i=offset; i<numpoints; i+=N)
{
j = i + half_N;            
temp.re = (w.re * f[j].re) - (w.im * f[j].im);
temp.im = (w.im * f[j].re) + (w.re * f[j].im);
f[j].re = f[i].re - temp.re;
f[i].re += temp.re;
f[j].im = f[i].im - temp.im;
f[i].im += temp.im;
}
wtemp = w.re;
w.re = wtemp * wp.re - w.im * wp.im + w.re;
w.im = w.im * wp.re + wtemp * wp. im + w. im ;
}
}
if (dir == -1)    ,         /*if. inverse fft, divide by numpoints */
for.(i=0;i<numpoints;i++)
{
f[i].re /= numpoints;
f[i].im /= numpoints;
}
}

Explanation of Listing 7.1

The first is the scrambling routine. Proper reordering of the data can take advantage of the periodicity and symmetry of recursive DFT computation. The scrambling routine is very simple. A bit reversed index is computed for each element in the data array. The data is then swapped with the data pointed to by the bit-reversed index. For example, suppose you are computing the FFT for an 8 element array. The data element at address 1 (001) will be swapped with the data at address 4 (100). Not all data is swapped since some indices are bit-reversals of themselves (000, 010, 101, and 111) (Figure 7.10).

The second part of the FFT function is the butterflies function. The butterflies function divides the set of data points down and performs a series of two point discrete Fourier transforms. The function is named after the flow graph that represents the basic operation of each stage: one multiplication and two additions (Figure 7.11).

Listing 7.1 shows the code to perform a forward and inverse FFT. The dir variable determines whether to do a forward or inverse transform. The twoD_FFT function performs a 1-dimensional FFT on all the rows and then the columns. The FFT function consists of only the bit scramble function and the butterflies function. The data is stored in an array of COMPLEX structures. COMPLEX is defined in ip.h and consists of a real and imaginary component.


Remember that the FFT is not a different transform than the DFT, but a family of more efficient algorithms to accomplish the data transform. Usually when one speeds up an algorithm, this speed up comes at a cost. With the FFT, the cost is complexity. There is complexity in the bookkeeping and algorithm execution. The computational savings, however, do not come at the expense of accuracy.

Now that you can generate image frequency data, it's time to display it. There are some difficulties to overcome when displaying the frequency spectrum of an image. The first arises because of the wide dynamic range of the data resulting from the discrete Fourier transform. Each data point is represented as a floating point number and is no longer limited to values from 0 to 255. This data must be scaled back down to put in a displayable format. A simple linear quantization does not always yield the best results, as many times the low amplitude data points get lost. The zero frequency term is usually the largest single component. It is also the least interesting point when inspecting the image spectrum. A common solution to this problem is to display the logarithm of the spectrum rather than the spectrum itself. The display function is

where c is a scaling constant and | H(u,v) | is the magnitude of the frequency data to display. The addition of 1 insures that the pixel value 0 does not get passed to the logarithm function.

Sometimes the logarithm function alone is not enough to display the range of interest. If there is high contrast in the output spectrum using only the logarithm function, you can clamp the extreme values. The rest of the data can be scaled appropriately using the logarithm function above (Figure 7.12).

Since scientists and engineers were brought up using the Cartesian coordinate system, they like image spectra displayed that way. An unaltered image spectrum will have the zero component displayed in the upper left hand corner of the image corresponding to pixel zero. The conventional way of displaying image spectra is by shifting the image both horizontally and vertically by half the image width and height. Figure 7.13 shows the image spectrum before and after this shifting. All spectra shown thus far have been displayed in this conventional way. This format is referred to as ordered (as opposed to unordered).

Now that we can view the image frequency data, how do we interpret it? Each pixel in the spectrum represents a change in the spatial frequency of one cycle per image width. The origin (at the center of the ordered image) is the constant term, sometimes referred to as the DC term (from electrical engineering's direct current). If every pixel in the image were gray, there would only be one value in the frequency spectrum. It would be at the origin. The next pixel to the right of the origin represents 1 cycle per image width. The next pixel to the right represents 2 cycles per image width and so forth. The further from the origin a pixel value is, the higher the spatial frequency it represents. You will notice that typically the higher values cluster around the origin. The high values that are not clustered about the origin are usually close to the u or v axis.

Purchase A Simplified Approach to Image Processing from Prentice Hall.

Related Links:
Introduction to the Frequency Domain
Discrete Fourier Transform

Reader Comments | Submit a comment »