Fast Fourier
Transform
Introduction:
In signal analysis, itˇ¦s often more convenient to do computations of the signal in frequency domain instead of time domain. Fourier Transform is one of the ways used to take a signal in time domain and transform it into frequency domain. Fast Fourier Transform is an algorithm that is created to speed up the calculations to transform the signal to frequency domain.
Algorithm:
*Here are some basic notations:
x (t) = time domain
X (f) = frequency domain
x (t) à X (f) = Fourier transform from time to frequency domain
The basic Fast Fourier transform looks like this:
Where
is cos(
) + j*sin()
A basic two input Fast Fourier transform looks like this:
So X (0) = x (0) + x (1), and X (1) = x (0) ˇV x (1), because cos(0)=1 and sin(0)=0.
For example:
Input signal {1, 2} à output signal {1+2, 1-2} = {3, -1}
or Input signal {0, 1} à output signal {0+1, 0-1} = {1, -1}
The 4-point FFT algorithm cascades two 2 point FFT algorithms and adds on a second stage called the combine 2 point FFT.
*NOTE: The graph says DFT (Discrete Fourier Transform), FFT (Fast Fourier Transform) is just a quicker algorithm
to do DFT, so itˇ¦s essentially the same thing.
The actual butterfly diagram looks like this:
Notice that when the input signal goes through the butterfly algorithm, it doesnˇ¦t go in order. There is a reordering that occurs before the signal is sent through the mathematical algorithm, and we will discuss that reordering later, but in short, you do a bit reversal. For example, the corresponding binary bits to 0-3 are 00, 01, 10, and 11. Doing bit reversals to them results in 00, 10, 01, and 11, which is 0, 2, 1, and 3. Thatˇ¦s the order you put in for the FFT. Here is an example of a 4 point FFT:
Input signal = {1, 2, 3, 4}
x (0) = 1 x (1) = 2 x (2) = 3 x (3) = 4
To get the signal at point A and point B (from the graph above), we do a 2 point FFT on x(0) and x(2):
{1, 3} à {1+3, 1-3}= {4, -2}
To get the signal at point C and D, we do another 2 point FFT on x(1) and x(3):
{2, 4} à {2+4, 2-4} = {6, -2}
At this point, our diagram should look like this:
Now we do the second stage of this FFT, the combined 2 point
FFT. From the diagram, we see that we take 4 and 6, and do a 2 point FFT with
, and then we take -2 and -2 and do a 2 point FFT with
.
Putting it all together:
Hereˇ¦s another example of a 4 point FFT:
Input signal = {10, 5, 8, 1}
Now we take a look at 8 point FFT, from the concept of a 4 point FFT, an 8 point FFT cascades two 4 point FFTs and then adds another stage called the 4 point combining stage:
Or, to break down the two 4 point FFTs, it would look like this:
Notice the bit reversal ordering that occurs to the input before itˇ¦s processed.
0 -> 000 -> 000 -> 0
1 -> 001 -> 100 -> 4
2 -> 010 -> 010 -> 2
3 -> 011 -> 110 -> 6
4 -> 100 -> 001 -> 1
5 -> 101 -> 101 -> 5
6 -> 110 -> 011 -> 3
7 -> 111 -> 111 -> 7
The butterfly diagram looks like this:
For Example:
Input signal = {1, 2, 3, 4, 5, 6, 7, 8}
After bit ordering:
After stage one:
After stage two:
Final stage:
So the final signal is:
X(0)=36
X(1)=-4+9.656j
X(2)=-4+4j
X(3)=-4+1.656j
X(4)=-4
X(5)=-4-1.656j
X(6)=-4-4j
X(7)=-4-9.656j
MIPS
Implementation:
The goal of the project is to implement the FFT function in MIPS assembly code and test it on a set of input samples. The number of input samples have to be a power of 2. To make it easy to start writing the project, we are providing skeleton code, so that you only have to fill in the functions for reordering the input signals (FFT_prep) and the actual function performing the FFT. The routines for performing fixed point multiplication and calculating Sin and Cos values and the functions for printing the result on the console are already provided.
The following are the files that are provided to you:
1. turnin.s: This is the main skeleton file where you need to fill in the FFT_prep and the FFT functions.
2. testcase1, testcase2, testcase3, testcase4: Skeleton codes for testing the program for different number of inputs.
3. testcase1_reorder, testcase2_reorder, testcase3_reorder, testcase4_reorder: Skeleton codes for testing only the reordering of the input signals
You need to write both the FFT_prep and the FFT functions in the file turnin.s for full credit.
Before we actually start implementing this function, here are a couple things we have to note:
First, the signal contains real and imaginary coefficients. How do we implement real and imaginary signal in MIPS? Here we take the approach of creating two arrays, one storing the real coefficients, and one storing imaginary coefficients. But if we do that, we have to make sure to synch the two arrays when we do operations on it. Then, after we do the math operations on the coefficients, we simply print out the ˇ§jˇ¨ onto the console.
Second, Sin/Cos involve decimal
numbers like .5, or .707 (
), how do we implement decimal numbers in MIPS where
registers only store positive integers? One approach is to use the floating
point registers provided by MIPS. However, the floating point registers have
their own floating point operations, which means we have learn a whole new set
of instructions to deal with floating point registers, so this approach is not
efficient to be implemented. The other approach is to use fixed point notation,
which we assume there is a decimal point somewhere in our 32 bit register. In
this case, we would do add/subtract exactly the same as normal. For example:
00010000 à (assume) 0001.0000
+) 00030000 à (assume) 0003.0000
----------------------------------------------
00040000 à 0004.0000
However, if we are to do a multiply, things get a little bit uglier:
00010000 à (assume) 0001.0000
x) 00030000 à (assume) 0003.0000
----------------------------------------------
00000004.00000000
Because of the fixed point is moved when we do a multiply, we have to find a way to shift the registers and figure out what number it is. The approach that is used is to assume the fixed point to be at the 4th spot just like above. When we multiply, instead of using ˇ§mulˇ¨ which takes 3 registers, we use ˇ§multˇ¨, which takes the two registers that are to be multiplied, and stores the result in two registers high and low. Then, we simply combine the first 4 digits from the high register, and the last 4 digits of the low register, and put it together to make our number. For example, if the above product 00000004.00000000 was split and put together, it would just become 00040000 stored in the register.
Notice in the problem however, we need to multiply two unreal numbers. Itˇ¦s not just 3 x 4, but itˇ¦s (3 + 4j ) ( 5 ˇV 6j ), how do we determine what the answer is since weˇ¦re not really keeping track of real and unreal numbers? If you break it down, the equation is like this (a+bj)(c+dj)= (ac-bd) + (ad+cb)j. So mixing together the shifting and combing registers for multiply, and breaking down the unreal multiplication, we have can multiply unreal numbers using fixed point notation. A function has been provided for you in assembly code.
Now
back to the original problem, how do we determine the sin and cos of a value?
In MIPS there are no built in sin and cos functions like C or Java. There has
been several approaches that people have used, and they all have advantages and
disadvantages. For example, we could build a lookup table containing all the
values from 0 to 90 for sin and cos, and every time the user calls sin or cos,
it just goes to that lookup table and finds the value and returns it. This
case, it uses the most space because the program needs to create memory to
store all the values of sin and cos. However, itˇ¦s the fastest to process,
because the CPU doesnˇ¦t need to execute extra instructions, it just needs to
find the value and return it. Another way is to use the
Our third concern is how we convert the fixed point notation back to decimal notation? The answer is you donˇ¦t need to. However, we do provide you with a function that prints out your signal onto the console so you donˇ¦t need to go on decoding Hex numbers and figuring out the fixed point notations. Once you finish the program, you may check your outputted values with matlab, the function for FFT is fft(). Hereˇ¦s a matlab example that does the FFT of the signal [1,2,3,4]:
f=[1,2,3,4]
fft(f)
You should see the answer come out, and you will be able to confirm your program to see if it works.