http://www.developer.com/

## Fun with Java, Understanding the Fast Fourier Transform (FFT) AlgorithmJanuary 11, 2005 Java Programming, Notes # 1486 ## Preface
You may find it useful to open another copy of this lesson in a separate browser window. That will make it easier for you to scroll back and forth among the different figures and listings while you are reading about them.
I recommend that you also study the other lessons in my extensive collection of online Java tutorials. You will find those lessons published at Gamelan.com. However, as of the date of this writing, Gamelan doesn't maintain a consolidated index of my Java tutorial lessons, and sometimes they are difficult to locate there. You will find a consolidated index at www.DickBaldwin.com. ## General DiscussionThe purpose of this lesson is to help you to understand how the Fast Fourier Transform (FFT) algorithm works. In order to understand the FFT, you must first understand the Discrete Fourier Transform (DFT). I explained how the DFT works in an earlier lesson entitled Fun with Java, How and Why Spectral Analysis Works. There are several different FFT algorithms in common use. In addition, there are many sites on the web where you can find explanations of the mechanics of FFT algorithm. I won't replicate those explanations. Rather, I will explain the underlying concepts that make the FFT possible and illustrate those concepts using a simple program. Hopefully, once you understand the underlying concepts, one or more of the explanations of the mechanics that you find on other sites will make sense to you.
The Fourier transform is most commonly associated with its use in transforming time-domain data into frequency-domain data. However, it is important to understand that there is nothing inherent in the Fourier transform regarding either the time domain or the frequency domain. Rather, the Fourier transform is a general-purpose transform that is used to transform a set of complex data in one domain into a different set of complex data in another domain. It is purely happenstance that it happens to be so valuable in describing the relationship between the time domain and the frequency domain.
For example, my first job after earning a BSEE degree in 1962 was in the Seismic Research Department of Texas Instruments. That is where I had my first encounter with Digital Signal Processing. In that job, I did a lot of work with Fourier transforms involving the time domain and the frequency domain. I also did a lot of work with Fourier transforms involving the space domain and the wave-number domain. Wave number is the name given to the reciprocal of wavelength for compression and shear waves propagating through a medium such as an iron bar, earth, water, or air, and also for electromagnetic waves such as radio and radar propagating through space.
For example, one of the things that we did was to compute two-dimensional Fourier transforms on diagrams representing weighted points in two-dimensional space. We would transform the weighted points in the space domain into points in the wave-number domain. The weighted points in the space domain represented the locations and amplifications of seismometers in a two-dimensional array on the surface of the earth. Each seismometer was amplified by a different gain factor and polarity. The amplified outputs of the seismometers were added together in various and complex ways intended to enhance signals and suppress noise.
In this case, the wave number was the reciprocal of the wave length of seismic waves propagating across the array. By plotting the results of the transformation in the wave-number domain, we could estimate which seismic waves would be enhanced and which seismic waves would be suppressed by the processing being applied to the seismometer outputs. We could also perform experiments on the computer where we caused the weights to vary with frequency, thus, allowing us to design and place digital filters on the seismometers to optimize the response of the array to earthquake signals while suppressing seismic noise associated with nearby cities and other sources of seismic noise.
I mention all of this simply to illustrate the general nature of the Fourier transform. Once again, the Fourier transform is simply a mathematical process that can be used to transform a set of complex values in one domain into a set of complex values in a different domain. Before getting into the details of this discussion, I want to refer you to a couple of excellent references on the FFT. Of course, you can find many more by performing a Google search for the keyword FFT.
Many of the images that you will see in this lesson were produced using an applet that I downloaded from http://sepwww.stanford.edu/oldsep/hale/FftLab.html. I changed the name of the applet class to cause it to fit into my file-naming scheme. I also made a couple of minor modifications to force its output to fit into this narrow publication format. Otherwise, I used the applet in its original form. This applet is extremely useful in performing FFT experiments very quickly and easily. I strongly recommend that you become familiar with it.
I also want to refer you to Chapter 12 of the excellent online book entitled
As mentioned earlier, the FFT algorithm is very complicated. I won't discuss the mechanics of the algorithm in this lesson. Rather, I will explain the underlying concepts that make the FFT algorithm possible. Hopefully after reading my explanation of the basic concepts, you will be able to understand the explanation of the mechanics of the algorithm provided by Smith and others.
The FFT algorithm is an algorithm that takes advantage of several reasonably well-know facts along with some less well-known facts. One of those facts is that the Fourier transform is a linear transform. By this, I mean that the transform of the sum of two or more input series is equal to the sum of the transforms of the individual input series. I will attempt to illustrate this in Figure 1, Figure 2, and Figure 3.
The images in these figures were produced using the FFT applet mentioned earlier. Figure 1 Transform of pulse with negative slope. An examination of Figure 1 shows that the display produced by the applet contains two sections. One section is labeled f(x) and the other section is labeled F(k). This is an interactive applet with the ability to transform the complex samples represented by f(x) into complex samples represented by F(k). Alternatively, the applet can be used to transform complex samples represented by F(k) into complex samples represented by f(x).
Each section contains two boxes, one labeled With one exception, each sample is represented by a black circle. In each box, one of the samples is represented by an empty circle. The empty circle represents an index value of zero. Samples to the right of the sample with the empty circle are samples at positive indices, and samples to the left of the sample with the empty circle are samples at negative indices.
A pair of values, one taken from the Real box and one taken from the Imaginary box, represents a complex sample. Any of the circles can be interactively moved up or down with the mouse. The value of each sample is represented by the distance of the corresponding circle from the horizontal line. When a change is made to the value of any sample belonging to either f(x) for F(k), the transformation is recomputed and the display of the other function is modified accordingly. If you modify the value of a sample in f(x), the values in F(k) are automatically modified to show the Fourier transform of f(x). If you modify the value of a sample in F(k), the values in f(x) are automatically modified to show the inverse Fourier transform of F(k). This is an extremely powerful interactive tool.
Most FFT algorithms require the input series to contain a number of complex samples that is a power of two such as 2, 4, 8, 16, 32, etc. Most FFT algorithms also produce the same number of complex samples in the output as are provided in the input. The FFT algorithm used in this applet is no exception to those rules. A pull-down list at the bottom of the applet lets the user specify 16, 32, or 64 complex samples for both the input and the output. All of the examples in this lesson use 16 complex samples for input and output.
The applet also provides a check box that allows the user to cause the origin
The other pull-down list and the button at the bottom of the applet provide other control features that don't need to be discussed here. I strongly urge you to download this applet and experiment with it. The results can be very enlightening.
Having discussed the features of the interactive FFT tool that I used to produce many of the images in this lesson, it is time to get back to the discussion of the Fourier transform as a linear transform. The fact that the Fourier transform is a linear transform is illustrated in Figure 1, Figure 2, and Figure 3. In these three figures, the input series is shown in the real area in the upper left. For simplification, the values of the imaginary part of the input series shown in the upper right are all zero. Also, for simplification, the zero origin is shown in the center by the value with the empty circle. The real and imaginary parts of the transform output are shown in the bottom of each figure. Figure 1 shows an input series consisting of a pulse that starts with a high value at the origin and extends down and to the right for five samples, ending in a large negative value. This input series produces a rather complicated transform output series, as can be seen in the bottom two boxes in Figure 1. I will come back to a discussion of the transform output later.
Figure 2 shown an input series consisting of a pulse that begins with a large negative value four samples to the left of the origin and extends up and to the right ending with a large positive value at the origin. The input series in Figure 2 is the mirror image of the input series in Figure 1 relative to the origin. Figure 2 Transform of pulse with positive slope.
Once again, the output from the transform of the input series is shown in the bottom two boxes of Figure 2. A comparison of the real part of each of the transforms for Figure 1 and Figure 2 shows that the real parts are the same, at least insofar as I was able to control the input by interactively adjusting the locations of the circles using the mouse. A comparison of the imaginary part of each of the transforms shows that the imaginary parts are the same except for the algebraic sign of each of the values in the imaginary part. The algebraic sign of each of the values in Figure 2 is the reverse of the algebraic sign of each of the values in Figure 1.
To demonstrate that the Fourier transform is a linear transform, I will create a new input series that is the sum of the input series from Figure 1 and Figure 2. I will show that the transform of the sum is the sum of the transforms. This is illustrated in Figure 3. Figure 3 Transform of the sum of two pulses.
Figure 3 shows an input series that is the sum of the individual input series from Figure 1 and Figure 2. This produces a pulse that is symmetric around the origin indicated by the value with the empty circle.
Note that the display of the transform values produced by this applet is normalized so as to keep them in a reasonable range for plotting. As a result, absolute values don't have much meaning. Only relative values have meaning. The real part of the transform of the input series in Figure 3 has the same shape as the real parts of the transforms of the input series in Figure 1 and Figure 2. This is what would be produced by adding the real parts of the transforms of the pulses in Figure 1 and Figure 2, and then normalizing the result.
The imaginary part of the transform of the input series in Figure 3 is zero at all sample values. This is what would be produced by adding the imaginary parts of the transforms of the input series in Figure 1 and Figure 2.
Thus, Figure 1, Figure 2, and Figure 3 demonstrate that the transform of the sum of two or more input series is equal to the sum of the transforms of the individual input series. The Fourier transform is a linear transform.
The real part of the transform of a single real sample with a shift relative to the origin has the shape of a cosine curve with a period that is proportional to the reciprocal of the shift. Negative sample values produce cosine curves with negative amplitudes. The imaginary part of the transform of a single real sample with a shift relative to the origin has the shape of a sine curve with a period that is proportional to the reciprocal of the shift. Negative sample values produce sine curves with negative amplitudes. The magnitude of the transform is the square root of the sum of the squares of the real and imaginary parts at each output sample point. For the case of a single input sample with a shift, that magnitude is constant for all output sample points and is proportional to the absolute value of the sample. The above facts are illustrated in Figure 4, Figure 5, Figure 6, and Figure 7. Figure 4 Transform of a real single sample with no shift.
Figure 4 shows the transform of a single real pulse with a shift of zero relative to the origin.
Although it isn't obvious, the real part of the transform in Figure 4 has the shape of a cosine curve with a period that is the reciprocal of the shift. Because the shift is zero, the period of the cosine curve is infinite, producing real values that are constant at all output sample values. Similarly, the imaginary part of the transform in Figure 4 has a shape that is a sine curve with an infinite period. Thus, it is zero at all output sample values.
Figure 5 shows the transform of a single real sample with a negative value and a shift of one sample interval relative to the origin. Figure 5 Transform of a real single sample with a shift equal to one sample interval and a negative value.
The shape of the real part of the transform output is an upside down cosine curve. It is upside down because it has a negative amplitude. This is caused by the fact that the input sample has a negative value. The shape of the imaginary part of the transform is an upside down sine curve.
This transform program computes real and imaginary values from zero to an output index that is one output sample interval less than the sampling frequency. The number of output values is equal to the number of samples in the input series. This is very typical of FFT algorithms. In this case, I set the applet up to accept sixteen input samples and to produce sixteen output samples.
For the moment, lets think in terms of time and frequency. Assume that the input series f(x) is a time series and the output series F(k) is a frequency spectrum. To make the arithmetic easy, let's assume that the sampling interval for the input time series in the upper left box of Figure 5 is one second. This gives a sampling frequency of one sample per second, and a total elapsed time of sixteen seconds. The sine and cosine curves in Figure 5 each go through one complete period between a frequency of zero and the sampling frequency, which one sample per second. Thus, the period of the sine and cosine curves along the frequency axis is one sample per second. This is the reciprocal of the time shift of one sample interval at a sampling frequency of one sample per second. Stated differently, the number of periods of the sine and cosine curves in the real and imaginary parts of the transform between a frequency of zero and a frequency equal to the sampling frequency is equal to the shift in sample intervals. A shift of one sample interval produces sine and cosine curves having one period in the frequency range from zero to the sampling frequency. A shift of two sample intervals produces sine and cosine curves having two periods in the frequency range from zero to the sampling frequency, etc. This is illustrated by Figure 6.
Figure 6 shows the transform of a real single sample with a shift equal to two sample intervals and a positive value. Figure 6 Transform of a real single sample with a shift equal to two sample intervals and a positive value. The real part of the transform has the shape of a cosine curve with two complete periods between zero and an output index equal to the sampling frequency. The imaginary part of the transform has the shape of a sine curve with two complete periods within the same output interval. This agrees with the conclusions stated in the previous section.
Finally, Figure 7 shows the transform of a real single sample with a shift equal to four sample intervals. Figure 7 Transform of a real single sample with a shift equal to four sample intervals and a positive value. The cosine and sine curves that represent the real and imaginary parts of the transform each have four complete periods between zero and an output index equal to the sampling frequency.
The main point is that if you know the value of a single real sample and you know its position in the series relative to the origin, you can write equations that describe the real and imaginary parts of the transform of that single sample without any requirement to actually perform a Fourier transform. Those equations are simple sine and cosine equations as a function of the units of the output domain. This is an important concept that contributes greatly to the implementation of the FFT algorithm.
The FFT algorithm is an algorithm that transforms a series of complex values in one domain into a series of complex values in another domain. The images in the figures discussed so far indicate a transformation of a complex function given by f(x) into another complex function given by F(k). There is nothing in these images to indicate anything about time and frequency. If the complex part of the input series f(x) is not zero, things get somewhat more complicated. For example, the real and imaginary parts of the transform of a single delayed sample having both real and imaginary parts are not necessarily cosine and sine curves. This is illustrated in Figure 8. Figure 8 Transform of a complex single sample with a shift equal to two sample intervals. Figure 8 shows the results of transforming a single sample having both real and imaginary parts and a shift of two sample intervals. Although both the real and imaginary parts of the transformed result have the shape of a sinusoid, neither is a cosine curve and neither is a sine curve. Both of the curves are sinusoidal curves that have been shifted along the horizontal output axis moving their peaks and zero crossings away from the origin.
Because the Fourier transform is a linear transform, you can transform the real and imaginary parts of the input separately and add the two resulting transforms. The sum of the two transforms represents the transform of the entire input series including both real and imaginary parts. The program that I will discuss later takes advantage of this fact. Even for a complex input series, if you know the values of the real and imaginary parts of a sample and you know the value of the shift associated with that sample, you can write equations that describe the real part and the imaginary part of the transform results.
That brings us to the crux of the matter. Given an input series consisting of a set of sequential samples taken at uniform sampling intervals, we know how to write equations for the real and imaginary parts that would be produced by performing a Fourier transform on each of those samples individually.
We know that we can consider the input series to consist of the sum of the individual samples, each having a specified value and a different shift. We know that the Fourier transform is a linear transform. Therefore, the Fourier transform of an input series is the sum of the transforms of the individual samples. If we are clever enough, we can use these facts to develop a computational algorithm that can compute the Fourier transform of a time series much faster than can be obtained using a brute force DFT algorithm. Fortunately, some very clever people have already developed that algorithm. It goes by the name of the Fast Fourier Transform, or FFT algorithm.
In truth, there are several different forms of the FFT algorithm, and the mechanics of each may be slightly different. At least one, and probably many of the algorithms operate by performing the following steps: - Decompose an N-point complex series into N individual complex series,
each consisting of a single complex sample. The order of the decomposition
in an FFT algorithm is rather complicated. It is this order of
decomposition, and the order of the subsequent recombination of transform
results that causes the FFT algorithm to be so fast. It is also that order
that makes the algorithm somewhat difficult to understand. Note that the program that
I will discuss later
implement that special order of decomposition and recombination.**does not** - Calculate the transform of each of the N complex series, each consisting of a single complex sample. This treats each complex sample as if it is located at the origin of a complex series. This step is trivial. The real part of the transform of a single complex sample located at the origin of the series is a complex constant whose values are proportional to the real and imaginary values that make up the complex sample. Since the complex input series consists of only one complex sample, there is only one complex value in the complex transform.
- Correct each of the N transform results to reflect the original position of the complex sample in the input series. This involves the application of sine and cosine curves to the real and imaginary parts of the transform. This step is usually combined with the recombination step that follows.
- Recombine the N transform results into a single transform result that represents the transform of the original complex series. This is a very complicated operation in a real FFT algorithm. It must reverse the order of decomposition in the first step described earlier. As mentioned earlier, it is the order of the decomposition and subsequent recombination that minimizes the arithmetic operations required and gives the FFT its tremendous speed. The program that I will discuss later does not implement the special order of decomposition and recombination used in an actual FFT algorithm.
## A Sample ProgramI want to emphasize at the outset that this program DOES NOT implement an FFT algorithm. Rather, this program illustrates the underlying signal processing concepts that make the FFT possible in a form that is more easily understood than is normally the case with an actual FFT algorithm.
In summary, a typical FFT algorithm performs the following processes: - Decompose an N-point complex series into N individual complex series, each consisting of a single complex sample.
- Recognize that the complex transform of a single complex sample is equal to the value of the complex sample.
- Correct the transform for each complex sample to reflect the original position of the complex sample in the input series.
- Recombine the N transform results into a single transform result that represents the transform of the original complex series.
This program performs each of the processes listed above. However, it does not perform those processes in the special order used by an FFT algorithm that causes the FFT algorithm to be able to perform those processes at very high speed.
The decomposition process in this program takes the complex samples in the
order that they appear in the input complex series.
No attempt is made to manage the decomposition and the subsequent recombination in the manner of a true FFT algorithm. Therefore, this program is designed to illustrate the processes involved, and is not designed to provide the speed of a true FFT algorithm. This program was tested using SDK 1.4.2 under WinXP.
As is my usual approach, I will discuss and explain this program in fragments. A complete listing of the program is provided in Listing 9 near the end of the lesson. The program begins in Listing 1, which shows the beginning of the controlling
class named
The first statement in the The purpose of an object of the I will put the
Listing 2 presents the beginning of the class named
The The third parameter is a scale factor that is applied to the transform output in an attempt to keep the values in a range suitable for plotting if desired. The last two parameters are references to array objects of type
The body of the
Each complex value in the incoming arrays represents both a complex sample and the transform of that complex sample under the assumption that the complex sample appears at the origin of the input series.
The method named Then the method named
The
This method accepts an incoming complex sample value and the position in the series associated with that sample. The method corrects the real and imaginary transform values for that complex sample to reflect the specified position in the input series. After correcting the transform values for the sample on the basis of position, the method updates the corresponding real and imaginary values contained in array objects that are used to accumulate the real and imaginary values for all of the samples. References to the array objects are received as input parameters. Outgoing results are scaled by an incoming parameter in an attempt to cause the output values to fall within a reasonable range in case someone wants to plot them. The incoming parameter named Hopefully this explanation will make it possible for you to understand the code in Listing 4. Note in particular the use of the Note the use of the Also note how the correction is made separately on the real and imaginary parts of the input. This produces results similar to those shown in Figure 7 after those results are added in the accumulators.
Returning now to the
Note that for Case A, the input complex series contains non-zero values only in the real part. Also, most of the values in the real part are zero.
Case A is shown in graphic form in Figure 9. As you can see, the input series consists of two non-zero values in the real part. All the values in the imaginary part are zero. Figure 9 Case A. Transform of a real sample with two non-zero values. The real part of the transform of the complex input series looks like one cycle of a cosine curve. All of the values in the imaginary part of the transform result are zero.
As you saw in Listing 5, the code in the
If you plot the real and imaginary values in Figure 10, you will see that they match the transform output shown in graphic form in Figure 9.
The code from the
Case B is shown in graphical form in Figure 11. Figure 11 Case B. Transform of a simple complex series.
The output from the code in Listing 6 is shown in Figure 12.
If you plot the values for the real and imaginary parts from Figure 12, you will see that they match the real and imaginary output shown in Figure 11.
The code extracted from the
The complex input series for Case C is a little more complicated than that for either of the previous two cases. Note in particular that the input complex series contains non-zero values in both the real and imaginary parts. In addition, very few of the values in the complex series have a value of zero.
Case C is shown in graphic form in Figure 13. Figure 13 Case C. Transform of a more complicated complex series.
One of the interesting things to note about Figure 13 is the similarity of Figure 13 and Figure 5. These two figures illustrate the reversible nature of the Fourier transform. If I had used a positive input real value instead of a negative input real value in Figure 5, the input of Figure 5 would look exactly like the output in Figure 13, and the output of Figure 5 would look exactly like the input of Figure 13. With that as a hint, you should now be able to figure out how I used a mouse and drew the perfect sine and cosine curves in Figure 13. In fact, I didn't draw them at all. Rather, I used my mouse and drew the output, and the applet gave me the corresponding input automatically.
The output produced by the code in Listing 7 is shown in Figure 14.
If you plot the real and imaginary input values from Listing 7, you will see that they match the input values in Figure 13. If you plot the real and imaginary output values in Figure 14, you will see that they match the output values shown in Figure 13. Listing 7 signals the end of the
Listing 8 shows the code for a simple method named
Listing 8 also signals the end of the controlling class named ## Run the ProgramI encourage you to copy and compile the program that you will find in Listing 9. Experiment with different complex input series. I also encourage you to download the applet from http://sepwww.stanford.edu/oldsep/hale/FftLab.html and experiment with it as well. Compare the numeric output produced by this program with the graphic output produced by the applet. I also encourage you to read what Steven Smith has to say about the FFT algorithm at http://www.dspguide.com/ch12.htm. Pay particular attention to his explanation of the order in which the input series is decomposed and the order in which the individual transform outputs are recombined. Finally, I encourage you to examine the source code for the applet. You will find that source code at http://sepwww.stanford.edu/oldsep/hale/FftLab.java. Concentrate on that portion of the source code that performs the FFT. Hopefully, what you have learned in this lesson in addition to what you learn from Steven Smith's book will make it easier for you to understand the source code for the FFT. ## SummaryNow that you understand those concepts, you should be able to better understand explanations of the mechanics of the FFT algorithm that appear on various websites. ## Complete Program ListingA complete listing of the program is provided in Listing 9 below.
Copyright 2004, Richard G. Baldwin. Reproduction in whole or in part in any form or medium without express written permission from Richard Baldwin is prohibited. ## About the authorRichard Baldwin
is a college professor (at Austin Community College in Austin, TX) and
private consultant whose primary focus is a combination of Java, C#, and
XML. In addition to the many platform and/or language independent benefits
of Java and C# applications, he believes that a combination of Java, C#,
and XML will become the primary driving force in the delivery of structured
information on the Web.
-end- |