Toggle nav

Why Should We Even Care?

The fast fourier transform is a powerful tool, and I believe everyone should at least understand the basics. The theory behind the fourier transforms is daunting and only a few require to understand the mechanisms of such complexity, however we can enjoy the benefits quite easily.

This post will not have too much text, but rather programs and simple short explanations on how it works. Seeing the programs working, and analyzing the results will give enough experience to play around with it, maybe even making digital filters.

First Things First

Let´s see, what does the FFT do? Well, it takes a signal in the time domain and converts it to the frequency domain. Not understood yet? In other words, imagine that we have a sinewave with a frequency of 1000Hz. So, if we calculate the FFT, the result will be a bunch of complex numbers that apparently don´t have any meaning, but in reality the new values are amplitudes of the signal for each frequency range defined by a single value. Therefore a pure 1000Hz sinewave will have a single non-zero value in the 1000Hz "bin" while others shall be zero.

Bin? Bin or bucket, I already saw different designations to it, but let´s go for bin. The frequency spectrum will be shown subdivided in equal segments, and each segment is a bin that contains a certain range of frequencies. In practical terms, to calculate the FFT of one signal it´s require to sample it, for example using an micro ADC, with a fixed time period which is twice as fast as the maximum frequency signal we want to sample. Remember the 1000Hz signal? Well, we need to sample the signal with at least a frequency of 2000Hz (Nyquist Rate). This means the signal shall be continuously sampled each 0.5ms (T=1/2000Hz).

Ok, so far we understood that a signal is converted from a time domain to frequency domain and the maximum frequency that we can represent is half the rate we can sample. Now let´s talk about length. The number of measured points is the FFT length and normally this value is a power of 2, for example, 64, 128, 256... The number of points we use to calculate the frequency spectrum will be directly connected to the size of each bin. Thus, the range of frequencies that we can differentiate.

A lot more can be said but those are the fundamentals, let´s now get to the programming bit. Note that FFT analysis can be done in Ms Excel, which is really useful to validate results. In fact, attached on the right bar I´ve included a sample file that will then be used for comparison. Our 64 point "dummy" data (C2:C65) has two frequency components, one at 0.1 Hz (G1) and the other at 0.25 Hz (H1). The data was "recorded" for 40s (J1) which gives us a sampling frequency of 1,6Hz (I1). The complex data (output of the FFT) is on cells F2:F64 and the magnitude is calculated and shown on cells E2:E65.

The data on the right of the spreadsheet/graph is just some fixed point FFT tests (no decimals). In microcontrollers due to their reduced calculation capacity it´s required to simplify all the math and use tricks to allow it to run in "real time". In micros suchs as PIC's, ATMEGA and others, the cosines and sines are typically not calculated, they are picked from a lookup table.

Discrete Fourier Transform

The FFT is just an incredibly faster way of calculating the DFT. Below is the code in C# (on the right tab the solution file) to calculate the complex DFT by correlation. This code in my laptop (Intel Core i7-2630QM - 2Ghz) takes roughly 100ms to calculate a 1024 points DFT and 1900ms for 4096 points.

DateTime currentDate = DateTime.Now;
int tickStart = Environment.TickCount & Int32.MaxValue;
UInt16 M = 10;
UInt16 N = (UInt16)Math.Pow(2, M); 
double[] XR;
XR = new double[N];
double[] XI;
XI = new double[N];
for (UInt16 i = 0; i < N; i++)
{
    // See FFT.xls to confirm values
    double ti = (double)i * (double)40 / (double)N;
    XR[i] = 3 * Math.Sin((2 * Math.PI * ti) / 10) + Math.Sin((2 * Math.PI * ti) / 4);
    XI[i] = 0;
}
double[] REX;
REX = new double[N];
double[] IMX;
IMX = new double[N];
for (UInt16 i = 0; i < N; i++)
{
    REX[i] = 0;
    IMX[i] = 0;
}
for(int k=0; k< N ; k++)
{
    for (int i=0; i< N; i++)
    {
        double SR = Math.Cos(2 * Math.PI * k * i / N);
        double SI = -Math.Sin(2 * Math.PI * k * i / N);
        REX[k] = REX[k] + XR[i] * SR - XI[i] * SI;
        IMX[k] = IMX[k] + XR[i] * SI + XI[i] * SR;
    }
}
int tickEnd = Environment.TickCount & Int32.MaxValue;
int elapsedTime = tickEnd - tickStart;
for (UInt16 i = 0; i < N; i++)
{
    Console.WriteLine("Value[{0}]={1:0.00} {2:0.00}i", i, REX[i], IMX[i]);
}
Console.WriteLine("Elapsed time= {0}", elapsedTime);
Console.WriteLine("Number of points= {0}", N);
Console.ReadLine();

Fast Fourier Transform

The FFT shall be faster than the DFT, so let´s put it to the test. Below the code and on the right tab the solution file. This code is super fast as expected, less than 1ms for 1024 and 4096 points. It only starts to be slower than a milisecond for 32768 points (2^15).

DateTime currentDate = DateTime.Now;
int tickStart = Environment.TickCount & Int32.MaxValue;
UInt16 M = 15;
UInt16 N = (UInt16)Math.Pow(2, M);
//UInt16 N = (UInt16)(1 << M);
UInt16 K = 0;
            
double[] REX;
REX = new double[N];
double[] IMX;
IMX = new double[N];
for (UInt16 i=0; i< N; i++)
{
    // See FFT.xls to confirm values
    double ti = (double)i * (double)40 / (double)N;
    REX[i] = 3*Math.Sin((2 * Math.PI * ti)/ 10)+ Math.Sin((2 * Math.PI * ti) / 4);
    //Console.WriteLine("REX[{0}]={1:0.0000} - {2:0.000}",i, REX[i], ti);
}
//Console.ReadLine();
//Console.Clear();
UInt16  NM1 = (UInt16)(N - 1);
UInt16 ND2 = (UInt16)(N >> 1);
UInt16 J = ND2;
double TR = 0;
double TI = 0;
for (UInt16 I=1; I> 1);
    }
    J = (UInt16)(J + K);
}

for (UInt16 L=1; L<= M; L++)
{
    UInt16 LE = (UInt16)Math.Pow(2, L);
    UInt16 LE2 = (UInt16)(LE >> 1);
    double UR = 1;
    double UI = 0;
    double SR = Math.Cos(Math.PI / LE2);
    double SI = -Math.Sin(Math.PI / LE2);
    for (J=1; J<= LE2; J++)
    {
        UInt16 JM1 = (UInt16)(J - 1);
        for (UInt16 I=JM1; I<= NM1; I= (UInt16)(I + LE))
        {
            UInt16 IP = (UInt16)(I + LE2);
            TR = REX[IP] * UR - IMX[IP] * UI;
            TI = REX[IP] * UI + IMX[IP] * UR;
            REX[IP] = REX[I] - TR;
            IMX[IP] = IMX[I] - TI;
            REX[I] = REX[I] + TR;
            IMX[I] = IMX[I] + TI;
        }
        TR = UR;
        UR = TR * SR - UI * SI;
        UI = TR * SI + UI * SR;
    }
}

int tickEnd = Environment.TickCount & Int32.MaxValue;
int elapsedTime = tickEnd - tickStart;
for (UInt16 i = 0; i < N; i++)
{
    Console.WriteLine("Value[{0}]={1:0.00} {2:0.00}i", i, REX[i], IMX[i]);
}
Console.WriteLine("Elapsed time= {0}",elapsedTime);
Console.WriteLine("Number of points= {0}",N);
Console.ReadLine();

// End of main

FFT for Microcontrollers

Fast Fourier Transforms in microcontrollers is very demanding for the processor due to the amount of calculations required. I made a FFT to run on PIC18F45K20 (which can run in other PIC18 series microcontrollers) using the multiplication instruction available on these type of MCU's, able to calculate the 64 points FFT in approximately 1.13ms. The attached code is a good base for development. The fastest way to use this code is to connect a PICkit 3 to a PIC18F45k20, build the solution, program it and run it. Go to the file "user.c" and add a breakpoint on line 105 test=0;. Run the program and check the output of the variable realOutput. Remember to see in the main.c file the instructions, connect the input signal on pin RA1.

Alternatively, to avoid using a real MCU, a real signal or even only to test the results, comment the first bit of code in main.c, namely from ConfigureOscillator(); all the way down up to calcFFT(); and read the instruction on the commented code right below. This test code uses a calculated wave and throw it to the FFT.