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 3
In 1965 two mathematicians, James Cooley and John Tukey, published a paper which proposed a method that made the calculation of the DFT much more efficient. In our example of a 16 point DFT (A DFT with 16 samples in it), we saw that the key to calculating it efficiently was to use the Divide and Conquer method:
- 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.
In the last post, we looked at how to compute a 2-point DFT for each of the sample pairs in the “Calculate” stage mentioned above and found that it was quite easy to do. However, before we move onto the “Combine” stage of the method, I would first like to take a look at how we would normally calculate a 4-point DFT without Cooley and Tukey’s help.
The 4-Point DFT
In a 4-point DFT, we divide our signal into 4 samples over its duration. We multiply the 4 samples in the signal by 4 samples of a Cosine wave and 4 samples of a Sine wave at frequency index 0. We then add all the results together to give ourselves 2 results, one for the Cosine multiplication and one for the Sine multiplication. This gives us the Cosine and Sine contributions of our signal for frequency 0. We then repeat this for each of the 4 different frequency indexes available in a 4-point DFT.
Let’s look at the source waves for each of the 4 frequencies, starting with frequency zero which I’m going to call b(0):
In the left hand graph above, we can see our signal (shown in green and which I’m going to call “x”) divided into 4 points at samples 0, 4, 8 and 12. The four points on our signal, I’m going to call therefore x(0), x(4), x(8) and x(12).
This is multiplied by a Cosine wave with a frequency of zero. (A Cosine wave with a frequency of zero is simply a straight line with an amplitude of 1). Any number multiplied by 1 remains unchanged. So when we add together all the multiplied points for our Cosine contribution at a frequency of zero, we’re simply going to get:
In the right hand graph above, we can see our signal (again shown in green) multiplied by a Sine wave with a frequency of zero. A Sine wave with a frequency of zero is simply a straight line with an amplitude of zero. Any number multiplied by zero equals zero, so there is no Sine contribution for this first frequency.
Next we’ll look at frequency 1 which I’m going to call b(1):
Looking at the left hand graph first, we see that at sample index 0, the Cosine wave is 1 so our signal at x(0) will simply be multiplied by 1. At sample index 4, the Cosine wave is zero so sample x(4) will be multiplied by zero. This is true also of sample x(12). Sample 8 is multiplied by -1. So for the Cosine contribution, we’re going to get:
–If you’ve been following the previous posts in this blog, the above equation might begin to stir certain memories–
Looking at the right hand graph, we see that at sample index 0, the Sine wave is zero, so x(0) will be multiplied by zero. This is true also for sample x(8). At sample index 4, the Sine wave has an amplitude of 1, so x(4) will be multiplied by 1. At sample index 12, the Sine wave has an amplitude of -1, so x(12) will be multiplied by -1. So for the Sine contribution, we’re going to get:
–any more memories stirring yet? —
Now I know that my old math professors would be screaming at me at this point for sloppy, ambiguous writing. We cannot have b(1) being equal to 2 completely different things as I have written it in the two equations above. So let’s go back to our complex notation and represent the Cosine contribution as a real number and the Sine contribution as an imaginary number (by multiplying it by “i”), so to keep everyone happy, I’ll rewrite my sloppy notation as follows:
If this last equation confuses you like it would have done me not so long ago, ignore it and just remember that the “i” is only there to try and keep the Cosine and Sine terms separate as I have done in the 2 equations preceding it. It’s not even all that important. What I’m trying to drive with all these graphs –and those cryptic comments in green– will hopefully begin to become clearer as we continue.
So now for frequency 2 which I’m going to call b(2):
Looking at the left hand graph, you’ll excuse me if I cut all the description and just go straight to the formula. Hopefully, you’ll be getting the gist of it by now. For the Cosine contribution for b(2):
–This equation looks similar to the equation for b(0) save for a couple of minus signs! —
Looking at the right hand graph… oh look, all the Sine samples are zero again so there is no Sine component for this frequency.
— bit like b(0) in that respect! —
Let’s move on to frequency 3 (the final frequency in a 4-point DFT) which I’m going to call b(3).
Looking at the left hand graph for the Cosine contribution:
— Now I’ve definitely seen that somewhere before! —
Looking at the right hand graph for the Sine contribution:
— OK so the minuses and pluses are switched round, but I’ve seen this equation before too! —
Again, to keep my old math lecturer happy, let’s combine these equations properly:
Alright! So enough of the green comments and all the nods and winks to previous blog posts! Let’s spell it out in words of one syllable and chuck in a nice, big red font for good measure.
The Calculations are repeating themselves!
We might need to change a plus sign to a minus sign here and there, but all we are doing is adding or subtracting our 4 samples x(0), x(4), x(8) and x(12) in different orders.
There is a definite similarity between the even frequencies b(0) and b(2) and the odd frequencies b(1) and b(3), and what is more, we’ve done some of these calculations before when we were calculating the 2-point DFTs in the previous blog post.
It seems quite plain that calculations involving x(0) and x(8) are repeating themselves with regularity, at least in the Cosine calculation, but x(4) and x(12) seem to come and go as they please. In order to help us ground what we have instinctively understood, we need to look now to the “Combine” stage of Cooley and Tukey’s algorithm.
3. Combine
The beauty of the Cooley Tukey method is that it allows us to treat each stage of the calculation as if we were only calculating a 2-point DFT. The reason that we can do it is that, at each stage, we have already done half of the work in the previous stage and we use the results we have already calculated to help us calculate the next stage.
As we saw last time, to calculate a 2-point DFT, we simply add together the 2 samples for the first frequency term and subtract them for the second. Let’s look at how we can still us this same method, combining what we calculated for the 2-point DFTs to calculate a 4-point DFT.
We have already calculated a 2-point DFT for samples 0 and 8, so we can reuse that result. We can look at this on a butterfly diagram as shown below:
Starting at the top left of the diagram, x(0) and x(8) are used to calculate the two frequency terms a(0) and a(1) for a 2-point DFT. These terms, a(0) and a(1) are then fed into the top 2 inputs of the 4-point DFT. However, as this is a 4-point DFT, we need 2 more samples to make up the numbers. These are samples x(4) and x(12) which produced the frequency components a(2) and a(3) from our second 2-point DFT in the bottom left hand corner of the diagram.
However, if we now look at a(2) and a(3) going into the second butterfly (the 4-point DFT) in the middle of the diagram, we see that something happens to these samples as they are fed in. They both get multiplied by a weird term which I’ve called “W”. This W term is known as a “Twiddle Factor” (also called a phase factor).
We spoke about “twiddle factors” at the end of the last post, It is now that the effects of these twiddle factors come in to play. In order to remind ourselves of what is happening, let’s look at the 2-point DFT for samples x(4) and x(12). Remember that in the 2-point DFT, there is no Sine contribution as both frequencies of Sine waves at the points we need to measure a 2-point DFT are zero. Therefore, we’re going to look at what happens to the Cosine wave.
In order to perform the 2-point DFT on samples x(4) and x(12), we had to shift our Cosine wave along the x-axis by π/2 radians as shown in the animation below.
This is because samples x(4) and X(12) did not lie at the correct points on the x-axis for a 2-point DFT to be performed on them. However, now we are dealing with a 4-point DFT, X(4) and x(12) do lie at the correct point, so we have to “twiddle” the result we obtained from the 2-point DFT by the amount we shifted the cosine wave. Watch what happens to the 2 samples on the Cosine wave at sample index 4 and sample index 2 as the wave is phase shifted back to where it should be. I’ve highlighted the samples with yellow circles.
You can see how they slide down the Cosine wave to the values they should be when the Cosine wave is in position for a 4-point DFT.
As we saw in the 2-point DFT:
- the frequency term a(2) is made up of the sum of the samples x(4) + x(12).
- the frequency term a(3) are simply the difference between the samples x(4) – x(12).
If both of these samples x(4) and x(12) are modified by the shifting back of the Cosine wave, then the corresponding frequency terms a(2) and a(3) will also be similarly affected.
Therefore before a(2) and a(3) can be fed into the 4-point DFT, they undergo an extra multiplication.
If we look at the term a(2) on the butterfly diagram above, as it gets fed into the bottom half of the 4-point DFT, it gets multiplied by a twiddle factor “W” with a 0 written just above it and a 4 written just below it.
As we see from the above equation, this twiddle factor is a complex number. The “0” above the W refers to the frequency index of a(2) in the 2-point DFT (remember a(2) was the first of 2 frequency terms, therefore its index is “0”). The “4” below the W refers to the number of points in the DFT we are currently calculating.
However, here’s the interesting thing.
In the equation above, if we disregard the “0” for a second, notice that 2π/4 is equal to π/2 which is exactly the number of radians we had to shift the cosine wave through when we were trying to perform the 2-point DFT on the 2 samples x(4) and x(12). So now we are calculating a 4-point DFT, in order to prepare the term a(2), we have shifted the Cosine wave back (by π/2 radians) to where it should be. It just so happens that this particular twiddle factor has no effect on the input a(2) as Cos(0 x π/2) = 1. If we look at what we are doing on a graph for this first frequency index:
As we are talking about a frequency of zero, it doesn’t matter by how much we phase shift the Cosine wave, it will always have an amplitude of 1. Also it doesn’t matter how much we phase shift the sine wave, it will always have an amplitude of zero.
So there is no imaginary part to the result as the sine term is zero. Sine(0 x π/2) = 0.
So putting this all together for our first frequency term of the 4-point DFT which I called b(0), the twiddle factor for a(2) was simply equal to 1 as we just saw so:
This is exactly the b(0) calculation we calculated when we did it the long way, only we already did half the work in the last stage as:
and:
so:
Cooley and Tukey’s algorithm just doesn’t stop making life easier! By combining the results from the previous stage 2-point DFT stage, we can now treat our 4-point DFT as if it was made up of interleaved 2-point DFTs as a(0) and a(2) have already been calculated. This means that no sooner have we calculated our term for b(0), when we get b(2) thrown in for free. We can reuse the a(0) and a(2) terms we used to calculate b(0), but just change the sign of a(2) to calculate b(2) just as we did with samples x(0) and x(8) in the 2-point DFT.
So putting this all together for our third frequency term of the 4-point DFT which I called b(2), the twiddle factor for a(2) is the same as before (=1), so:
This is exactly the b(2) calculation we calculated when we did it the long way, only we already did half the work in the last stage as:
and:
so:
We do exactly the same process for the other two frequency terms in the 4-point DFT, b(1) and b(3). The only difference is the twiddle factor.
Looking first at b(1), the frequency of our Cosine wave is 1, as we saw in the animation (here it is again):
At a frequency of 1 the twiddle factor does have an effect on samples x(4) and x(12), so a(3), which is made up of x(4) and x(12) will be similarly affected.
We already saw in the animation that the Cosine term at sample indexes 4 and 12 at this frequency is zero so for b(1) there is no Cosine component. However, look what happens to the Sine component as we shift the Sine wave back to where it should be for the 4-point DFT.
There is a Sine component as Sin (2π/4), which is the same as Sin (π/2), is equal to 1. This is why the “i” remains in the twiddle factor equation above as it signifies a Sine component with an amplitude of 1.
Using the butterfly diagram we can plug what we know back into the calculation for b(1).
For the Cosine component, a(3) is zero so:
and from before we know that:
For the first time, we have a Sine component in our signal, a(3) = 1, so for the Sine component:
and from before we know that:
So putting this all together:
This is exactly the b(1) calculation we calculated when we did it the long way as we can rewrite this as:
This being the Cooley Tukey algorithm, having calculated our frequency b(1), we get frequency b(3) for free simply by changing the sign of a(3).
Just as with the 2-point DFTs, we repeated the process to cover all the samples in the signal, we do exactly the same for the all the 4-point DFTs. This combines all the eight 2-point DFTs into four 4-point DFTs.
This is shown on a butterfly diagram as follows:
Consecutive groups of 2 samples from the 2-point DFTs combine together to form the 4 inputs of the 4 point DFTs. In each of the 4-point DFTs, the second pair of inputs are each multiplied by twiddle factors then their twiddled values interleaved with the first pair of inputs to calculate the result.
We have now covered all the relevant stages for the FFT. Next time we’ll apply the methods we have learned in this post to combine the four 4-point DFTs into two 8-point DFTs then into to one 16-point DFT which will give us our final result. We’ll then follow this up with a numeric example and actually calculate a 16-point FFT for a real life signal.
The Fast Fourier Transform Part 4 >>
