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.
In this post, we’re going to develop an algorithm to implement all that we have learned in the last few posts about the FFT. In developing this algorithm, I’ve started from the smallest part of the computation I could think of and worked outwards. The basic building block of the FFT is the “Butterfly” calculation. This calculation is iterated many times over the course of the FFT.
The snippets of code that appear in this post are written in Javascript.
Before we start, let’s define some terms:
Any size of FFT will be broken down into stages. For example, I’ve shown a 16-point FFT in the diagram above. The number of stages can be calculated by the following formula:
If you plug the number 16 into the FFTSize of the formula above you’ll find that there are 4 stages required to calculate the FFT as shown in the diagram above. These stages are numbered 0-3.
Each stage of the calculation has a number of groups of butterflies in it. In stage 0, there are 8 groups, each group containing 2 samples being fed into 1 butterfly.
In stage 1, there are 4 groups, each group containing 4 samples being fed into 2 overlapped butterflies.
In stage 2, there are 2 groups, each group containing 8 samples being fed into 4 overlapped butterflies.
…and so on.
The number of samples in each group can be calculated by the following formula:
The number of groups in each stage can be calculated by the following formula:
Each butterfly takes two samples and adds them together for the first term and subtracts them for the second term. Things are complicated slightly (but not much) by the fact that the second sample needs to be multiplied by a twiddle factor. This addition and subtraction repeats itself constantly throughout the computation of the FFT which means we are going to reuse the piece of code that does it a lot so it makes sense to program it into it’s own function. The generalized butterfly for any sample pair in the FFT can be calculated as follows.
x is the input sample
F is the output frequency term
FFTSize is the total number of samples going into the FFT
n is the sample index (0, 1, 2, 3, …, FFTSize-1)
N is the number of samples in each group of the current stage of the FFT
k is the order of the twiddle factor within each group ( 0, 1, 2, 3, …, (N/2)-1 )
W is the twiddle factor
However, although the butterfly lies right at the heart of the calculation, it itself is made up of 3 even more basic calculations. These are one complex add (to add the 2 samples together), one complex subtract (to subtract the 2 samples one from each other) and one complex multiply (to multiply the second sample by the twiddle factor).
So we are first going to develop 3 functions to do each of these 3 operations. Throughout the whole algorithm, I’m going to define my own data structure to store the complex numbers used throughout the calculation. This is a structure made up of 2 floating point numbers, one for the Real term and one for the Imaginary term as follows:
var ComplexNumber = {Real: 0, Imaginary: 0};
So lets start simple and write the function which adds together 2 complex numbers. To do this we need to add the the real and imaginary terms separately.
function ComplexAdd(ComplexNumber0, ComplexNumber1)
{
var ComplexResult = {Real: 0, Imaginary: 0};
ComplexResult.Real = ComplexNumber0.Real + ComplexNumber1.Real;
ComplexResult.Imaginary = ComplexNumber0.Imaginary + ComplexNumber1.Imaginary;
return ComplexResult;
}
The function accepts 2 complex numbers (ComplexNumber0 and ComplexNumber1), adds together the real and imaginary parts separately and returns the result (ComplexResult).
Writing the subtract function is just as simple. We only need to change the addition sign to a subtraction sign:
function ComplexSubtract(ComplexNumber0, ComplexNumber1)
{
var ComplexResult = {Real: 0, Imaginary: 0};
ComplexResult.Real = ComplexNumber0.Real - ComplexNumber1.Real;
ComplexResult.Imaginary = ComplexNumber0.Imaginary - ComplexNumber1.Imaginary;
return ComplexResult;
}
The third function is the complex multiply function. Remember that multiplying two complex numbers is like multiplying 2 brackets. We have to use the FOIL method I mentioned in the previous post.
function ComplexMultiply(ComplexNumber0, ComplexNumber1)
{
var ComplexResult = {Real: 0, Imaginary: 0};
var First;
var Outside;
var Inside;
var Last;
// First - Produces real result
First = ComplexNumber0.Real * ComplexNumber1.Real;
// Outside - Produces imaginary result
Outside = ComplexNumber0.Real * ComplexNumber1.Imaginary;
// Inside - Produces imaginary result
Inside = ComplexNumber0.Imaginary * ComplexNumber1.Real;
// Last - Produces real result multiplied by i-squared (i-squared = -1)
Last = -1 * ComplexNumber0.Imaginary * ComplexNumber1.Imaginary;
ComplexResult.Real = First + Last;
ComplexResult.Imaginary = Inside + Outside;
return ComplexResult;
}
Now that we have the 3 mathematical operations, we can write the code for the butterfly. The butterfly has to add together the first and second samples to produce the first frequency term and subtract the second sample from the first to produce the second frequency term. The second sample is always multiplied by the relevant twiddle factor. Therefore the inputs to the function are the two samples and the twiddle factor. All of these are complex numbers. The function returns an array containing the two frequency terms, again as complex numbers.
function Butterfly(Sample0, Sample1, TwiddleFactor)
{
var Frequency=[];
var TwiddledSample1=ComplexMultiply(Sample1, TwiddleFactor);
Frequency[0] = ComplexAdd(Sample0, TwiddledSample1);
Frequency[1] = ComplexSubtract(Sample0, TwiddledSample1);
return Frequency;
}
Next we need to have a function which calculates the twiddle factors. The number of twiddle factors depends on which stage of the calculation we are currently doing and can be calculated by raising 2 to the power of the index of the current stage. So, for example, if we are in the second stage of the FFT the stage index will be equal to 1 (remember the index starts from zero) and there will be 2 twiddle factors (). Therefore we need to have a “StageIndex” at the input to the function to tell it which stage of the calculation we need the twiddle factors for. The functions then returns an array of complex numbers which are the twiddle factors. The twiddle factors are calculated using the following formula:
where:
N is the number of samples in each group of the current stage of the FFT
k is the order of the twiddle factor within the current group.
For the second stage of the FFT, k will range between 0 and 1.
function CalculateTwiddleFactors(StageIndex)
{
var NumberOfTwiddleFactors = Math.pow(2,StageIndex);
var TwiddleFactors=Array(NumberOfTwiddleFactors);
var NumberOfSamples = Math.pow(2,StageIndex+1);
var i;
for (i=0; i<NumberOfTwiddleFactors; i++)
{
var ComplexNumber = {Real: 0, Imaginary: 0};
ComplexNumber.Real = Math.cos(2 * Math.PI * i / NumberOfSamples);
ComplexNumber.Imaginary = -1 * Math.sin(2 * Math.PI * i / NumberOfSamples);
TwiddleFactors[i] = ComplexNumber;
}
return TwiddleFactors;
}
Now we’ve got the butterfly sorted out with all its twiddle factors, we need to prepare the samples that are going to be input into the first stage of the FFT. Remember the samples have to be placed in a special order, a bit-reversed order, not the order they occur in naturally. Therefore the following function will reorder them for us:
function BitReversal(SamplesIn, NumberOfSamples)
{
var NewOrder;
var SamplesOut=[];
var NumberOfBits;
var SampleIndex;
NumberOfBits = Math.log(NumberOfSamples) / Math.log(2);
for (SampleIndex=0; SampleIndex<NumberOfSamples; SampleIndex++)
{
NewOrder = 0;
for (BitIndex=0; BitIndex<NumberOfBits; BitIndex++)
{
if ((SampleIndex < Math.pow(2, BitIndex)) == Math.pow(2, BitIndex))
NewOrder = NewOrder + Math.pow(2, (NumberOfBits - 1 - BitIndex));
}
SamplesOut[SampleIndex] = SamplesIn[NewOrder];
}
return SamplesOut;
}
The function accepts at its input, the array of samples to be reordered (I’ve called the array: “SamplesIN”). We also have to tell the function how many samples there are in the array (NumberOfSamples) for reasons that will become clear shortly.
What this function is doing is taking the bits of the index of the current sample and reversing them so that, for example we have 16 samples in all and the the index of the current sample is 8 (which in binary is 0100), the new index of this sample will now be 2 (which in binary is 0010).
The total number of samples in the FFT is important as the order will change depending on how many samples there are. This will affect the number of bits to be reversed. If there are 32 samples in my FFT for example, it takes 5 bits to describe the number 31 (remember the indexes start from 0 so the highest index will be 31 not 32). Therefore, if I bit reverse the index of sample index 8 (which in binary is 00100) in a 32-point FFT, the new index will remain 8 as reversing the bits of the number 8 in 5-bit binary gives (00100) which is the same as before.
The function then uses this new, bit-reversed index (I’ve called it “NewOrder” in the above function) to place the samples into a new array in a new order. I’ve called the new array: “SamplesOut”. This is then returned to the main function and the calculation of the FFT can commence.
So here is the FFT itself:
function FFT(SampleArray, FFTSize)
{
var NumberOfStages = Math.log(FFTSize) / Math.log(2);;
var DFTStage;
var SampleIndex;
var GroupIndex;
var NumberOfSamplesInGroup;
var NumberOfGroups;
var CombinedIndex;
var HalfOfSamplesInGroup;
var TwiddleFactors=[];
var Sample0;
var Sample1;
var TwiddleFactor;
// Reorder the Samples usinge bit reversal technique
SampleArray = BitReversal(SampleArray, FFTSize);
// Main FFT calculation loop
for (DFTStage=0; DFTStage<NumberOfStages; DFTStage++)
{
// Calculate the twiddle factors for this stage of the DFT
TwiddleFactors = CalculateTwiddleFactors(DFTStage);
// Prepare to organize the samples into groups
NumberOfSamplesInGroup = Math.pow(2,DFTStage+1)
NumberOfGroups = FFTSize / NumberOfSamplesInGroup;
HalfOfSamplesInGroup = NumberOfSamplesInGroup / 2;
// Perform the Butterfly calculation on each group
for (GroupIndex=0; GroupIndex<NumberOfGroups; GroupIndex++)
{
for (SampleIndex=0; SampleIndex<(NumberOfSamplesInGroup/2); SampleIndex++)
{
CombinedIndex = NumberOfSamplesInGroup * GroupIndex + SampleIndex;
// Prepare samples and twiddle factor for input into the butterfly
Sample0 = SampleArray[CombinedIndex];
Sample1 = SampleArray[CombinedIndex + HalfOfSamplesInGroup];
TwiddleFactor = TwiddleFactors[SampleIndex];
// Do the butterfly calculation
Results = Butterfly(Sample0, Sample1, TwiddleFactor);
// Place results back into the sample array ready for the next stage
SampleArray[CombinedIndex] = Results[0];
SampleArray[CombinedIndex + HalfOfSamplesInGroup] = Results[1];
}
}
}
return SampleArray;
}
The first thing to do is to work out how many stages we’re going to go through to calculate the FFT using the formula we mentioned at the beginning of this post:
We then setup a counter (DFTStage) to keep track of which stage we are currently calculating. Before we can begin the calculation, we have to reorder our samples so they are in the correct order for the FFT. We do this using the BitReversal function mentioned above. Once the input samples are in the right order, we can begin the main loop which loops through each stage of the FFT calculation.
We calculate all the twiddle factors for this stage then prepare two more loops, the second of which will run inside the first. The first loops through each group and within that loop sits another loop which runs the butterfly itself on each of the pairs of samples within the group.
The results are placed back into the input array ready for inputting into the next stage of the FFT until all the stages have been calculated. This array, which now contains all the results is returned at the end of the function.
The FFT returns an array of complex frequency terms which in itself doesn’t help us much. What we need to do now, is take the results of the FFT calculation and extract Magnitude and Phase data for each frequency index.
To calculate the magnitude we use Pythagoras:
function CalculateMagnitude(SampleArray, FFTSize)
{
var Magnitude=[];
var Real;
var Imaginary;
for (i=0; i<FFTSize; i++)
{
Real = SampleArray[i].Real;
Imaginary = SampleArray[i].Imaginary;
Magnitude[i] = Math.pow(Math.pow(Real, 2) + Math.pow(Imaginary, 2), 0.5);
}
return Magnitude;
}
To calculate the phase we use the inverse Tan function:
function CalculatePhase(SampleArray, FFTSize)
{
var Phase=[];
var Real;
var Imaginary;
for (i=0; i<FFTSize; i++)
{
Real = SampleArray[i].Real;
Imaginary = SampleArray[i].Imaginary;
Phase[i] = 180 * Math.atan(Imaginary / Real) / Math.PI;
}
return Phase;
}
I purposely wrote this algorithm in Javascript so that it could be run on a web browser and this is what we’ll be doing next time in my Javascript FFT calculator allowing you to plug in any signal into the input (so long as the number of samples in your signal is a power of 2) and calculate the FFT for it. The calculator allows you to see the results at each stage of the calculation as well as all the twiddle factors.
