Filming is currently underway on a special online course based on this blog which will include videos, animations and work-throughs to illustrate, in a visual way, how the Fourier Transform works, what all the math is all about and how it is applied in the real world.
Click here to reserve your free module
The module will be emailed to you the moment the course goes live.
The Fast Fourier Transform – Part 1
Those of you who have been following this blog since I started it last year will have noticed that this “final” post has been a long time in coming. The Fast Fourier Transform (FFT) is quite ingenious and it has taken me a lot of research to finally get my head round just how it works. My claim at the end of the last post that this would be my “final” post on the Fourier Transform was a little premature. The FFT is almost a whole subject in itself. Therefore this is Part 1 of another two or three posts dedicated to the FFT.
In the last post, we actually managed to get a reasonably good Fourier Transform of a real world signal. That is, we took the signal and managed to break it down into the frequency components that make it up. This is all well and good and you could program the algorithm into a computer and it would work. Computers are very good at repetitive math functions, which is what the Fourier Transform is. However, although computers can calculate things very fast, even they can be a little slow when the number of calculations gets very large.
If we look again at the Discrete Fourier Transform (DFT) equation in its polar form:
We have our sampled signal x(n) which is made up of N samples. In the example I’ve been using up till now, my signal had 40 samples in total (N=40).
In order to perform the Fourier Transform on this signal we first have to multiply every point on our sampled signal by every point on a sampled cosine wave and every point on a sampled sine wave at our first frequency index (k=0). This is known as a complex multiplication.
Secondly, we have to add together all the results of the two multiplication calculations to produce 2 summation results, one for the cosine term and one for the sine term at our first frequency index k=0. This is known as a complex addition.
We then have to repeat the above complex multiplication and complex addition calculations for every frequency index in the signal. As there are 40 samples, there will be 40 frequency indexes. This means that we are going to need to perform 40 complex multiply and add calculations for each of the 40 frequency indexes in our signal. This gives us a total of 40 x 40 = 1600 complex calculations.
So for just 40 samples, we need to do 1600 complex calculations. OK, so 1600 complex calculations isn’t much for today’s PCs, but any useful signal, like a sound signal, is going to be sampled at a much higher rate than 40Hz, and to do any meaningful sort of Fourier analysis of the signal, we are going to need a lot more than 40 samples. In the world of sound, we tend to sample signals at a rate of 44,100Hz (or higher) and send groups of 1024 samples at a time (or often much more than this) through our DFT. In such a case we are talking about over a million complex calculations for only 1024 samples. It is clear that the DFT is going to be very expensive on processing power.
We need some way of making the DFT calculation more efficient. Luckily for us, the repetitive nature of the sines and cosines within the DFT enables us to use a strategy known as “Divide and Conquer”.
Divide and Conquer
To demonstrate how the principle of “Divide and Conquer” works, let’s consider a simple problem that is not related to the Fourier Transform in any way.
Imagine we have a list of 16 numbers and we want to find which is the largest number in the list.
The most obvious way is to take a sequential approach – take each number and ask if it is greater than each other number. So starting with the first number, 8:
Question 1: Is 8 greater than 36?
Answer 1: No
Conclusion: 8 is not the largest number.
Question 2: Is 36 greater than 20?
Answer 2: Yes
Question 3: Is 36 greater and 54?
Answer 3: No
Conclusion: 36 is not the largest number.
Question 4: etc…
We would carry on with this question, answer and conclusion method until we had compared all the numbers in the list. Depending on where the largest number actually was in the list, this could take us a very long time. A far better way of doing it is to divide the problem into smaller and smaller problems until the problem becomes so small, the answer is trivial to find.
So let’s divide our list of 16 values into 2 groups of 8 values. The first 8 values will be in my first group and the last 8 values will be in my second group.
We’ll then divide the 2 groups of 8 values into 4 groups of 4 values in the same way.
…and again into 8 groups of 2 values.
Once we are down to only a pair of values in each group, we look for the largest number in each pair. This is very easy to do as it requires only one “is this number larger than this number” question per sample pair. I’ve marked the largest value in each group in yellow.
We then keep this number and discard the other. This we do for each pair. We then group together, again in pairs, the numbers we are left with. We do this again and again, finding the largest number in each pair, discarding the other number and grouping the results in new pairs, until we are left with only one number – the answer to our problem.
All in all, using “Divide and Conquer”, we have only had to ask “is this number bigger than this number?” 15 times whereas we would have had to ask the question over 200 times if we had been using the sequential method (depending on where, in the original list of 16 numbers, the largest value was located).
So “Divide and Conquer” is a much more efficient way of getting our result. How can we adapt this approach to help us calculate the Discrete Fourier Transform in a more efficient way?
Methodology
The first thing we have to do is to devise a methodology for using the “Divide and Conquer” technique. We could break down how we used the technique into the following stages:
- SPLIT – Keep splitting the samples into groups of half the number of samples until you are left with only a pair of samples in each group.
- CALCULATE – Perform your algorithm on each of the sample pairs.
- COMBINE – Use the results from the calculation you just performed to form the input of the next stage of the problem.
- REPEAT – Keep repeating 2 and 3 until you have an overall answer.
1. Split
In 1965, 2 mathematicians, James Cooley and John Tukey published a paper which revolutionized how the DFT was calculated. Interestingly, they weren’t the first to find a way of making the DFT more efficient. Carl Friedrich Gauss, the German mathematician famous for his work in statistics, had used it to interpolate the trajectories of asteroids over a century and a half earlier, however, his work was ahead of its time and needed the invention of the computer to give it a wider application.
At the heart of Cooley and Tukey’s discovery lies a very simple observation: Cosine and Sine waves repeat themselves every 2π radians. They wrote this observation as a mathematical identity:
Remember that the above equation is the same as saying:
What this identity is saying is that if you halve the number of samples in your DFT, the magnitude of the Cosine wave would be the same at N samples as it will be at N/2 samples. The same is true of a Sine wave. You can see this graphically below:
A Cosine wave completes one cycle every 2π radians. If we want our Cosine to complete 2 cycles, it is going to take 4π radians to do it. If we want 2 cycles of our cosine wave to complete in 16 samples (N=16), then we have to divide 4π by 16. Therefore the division by N (the number of samples in our DFT) is simply a way of mapping the sample index onto the x-axis of a Cosine graph as you can see in the graph above. The sample indexes are the black numbers and the radians are the red numbers. If we were now to halve the number of samples (N=8), this is the same as halving the length of the x-axis to 2π radians which would allow only 1 cycle of our Cosine wave to complete. Now you can see why the magnitude of the Cosine wave is the same. It’s because the Cosine wave repeats itself with a period of 2π radians. So if I look at the magnitude of my Cosine wave at N=16 samples, it will be the same magnitude as if I looked at it at N/2 = 8 samples as shown by the yellow circles on the graphs above. Exactly the same is true of a Sine wave.
This repeating nature of the Cosine and Sine waves means that as we increase their frequency, the amplitude of the wave at certain points is going to recur. This is demonstrated in the following animation. The blue line is the Cosine wave and the green line is our signal. Remember, the Fourier Transform multiplies the signal by the Cosine wave at each sample point. The blue dots show the amplitude of the Cosine wave at each sample point and the green dots show the amplitude of the our signal at each sample point.
The animation shows the a Cosine wave increasing in frequency from 0 Hz to 8 Hz. I have marked 4 points on the graph with yellow circles. The first pair of points is at a sample index of 0. As the frequency of the Cosine wave increases, these 2 points never change their value.
The second pair of points is at sample index 8. As the frequency increases, notice how that particular sample point on the Cosine wave falls out of the yellow circle as the amplitude of the sample changes with the increase of frequency. However although initially it falls away from the yellow circle, as the frequency increases even more, it begins to come back so that at a frequency of 2 Hz, the amplitude of the Cosine wave is exactly the same as it was at a frequency of 0 Hz. The same happens again at 4 Hz, 6 Hz and 8 Hz.
This means that the multiplication of the samples of the Cosine wave and the signal at sample index 0 and the multiplication of the samples of the Cosine wave and the signal at sample index 8 are going to be exactly the same at 0 Hz, 2 Hz, 4 Hz, 8 Hz etc. In short, every 2 Hz. If these values recur, we only need do the multiplication once then remember the result and use it again whenever the same values recur.
Now let’s do the same for another 2 pairs of sample points: sample index 4 and sample index 12.
This times the recurring values happen at 0 Hz, 4 Hz, and 8 Hz, every 4 Hz.
If we continue doing this for all of the 16 points at all the frequencies on the graph we find that each of the pairs of points we need to multiply and add together recur at differing frequency intervals. Therefore, all we need to do is the first time each pair of points occur, we take the points, multiply them together, then add the results and remember the answer for use later on when the same values recur instead of having to calculate them again.
This saving on the number of complex multiplication operations we’re going to have to do which is what makes the Fast Fourier Transform fast.
If we count up all the multiply and add operations we have to do using the Fast Fourier Transform, we’ll find that instead of needing 256 complex operations for a 16-point DFT, we now only need 64. That’s quite a saving.
So now, armed with Cooley and Tukey’s observation, let’s go back to the divide and conquer method we spoke about before. We’re going to divide the samples in such a way that the order of the sample pairs ends up with the pair of repeating sample points that recur the most are first and the pair the recur the least often are last. It just so happens that in order to do this, we have to split the samples up into even and odd samples as shown in the animation below:
On the bottom line we are left with 8 groups of 2 samples whose order from left to right reflects how often their multiplications recur. Samples 0 and 8 recur most often while samples 7 and 15 recur least often.
If a computer were to do this splitting operation, there is a very easy way to get from the single group of 16-samples down to the eight groups of 2 samples in one operation. The method used is known as “Bit Reversal”. If we rewrite the sample indexes in binary form, then reverse the order of the bits, this will give us the order of samples as they appear in the eight groups at the bottom of the above diagram.
| Original Order | Binary | Bit Reversal | New Order |
|---|---|---|---|
| 0 | 0000 | 0000 | 0 |
| 1 | 0001 | 1000 | 8 |
| 2 | 0010 | 0100 | 4 |
| 3 | 0011 | 1100 | 12 |
| 4 | 0100 | 0010 | 2 |
| 5 | 0101 | 1010 | 10 |
| 6 | 0110 | 0110 | 6 |
| 7 | 0111 | 1110 | 14 |
| 8 | 1000 | 0001 | 1 |
| 9 | 1001 | 1001 | 9 |
| 10 | 1010 | 0101 | 5 |
| 11 | 1011 | 1101 | 13 |
| 12 | 1100 | 0011 | 3 |
| 13 | 1101 | 1011 | 11 |
| 14 | 1110 | 0111 | 7 |
| 15 | 1111 | 1111 | 15 |
In the last row of the split. we only have 2 samples in each group so N=2. To work out a 2-point DFT is simplicity itself. We’ll see how this works when we look at the “Calculate” stage of our methodology in the next post.
The Fast Fourier Transform Part 2 >>
