Wednesday 2 January 2013

DSP Tech Brief : I Have Downloaded xyz FFT (Fourier) Transform Algorithm But How Do I Know I Have Implemented It Correctly And There Are No Bugs ?


Implementing an FFT on any RISC/VLIW architecture is a relatively simple task but getting it optimized for a particular architecture can be really time consuming which is why most of the time the first place I look is the device manufacturers web site to see what their Apps team have come up with.

When ever I use a new FFT function for the first time I always verify the results of processing both a sine and cosine wave against what is expected from the underlying maths.

If you perform a DFT on a cosine wave then the results should be purely real, as shown in the following diagram. The top graph shows 32 samples of a cosine, with the middle and bottom graphs showing the real and imaginary components of the frequency domain output from the FFT. Notice that there is no energy in the imaginary component and the real energy is located in two impulses representing the magnitudes of the two phasors that are the components of the original cosine waveform.


Similarly if you do the same for a cosine wave then the results should be purely imaginary, as follows :

Note the magnitude of the two phasors that form the sine wave.

If this is the first time you have seen graphs like this then can you explain the magnitude (16) of all four impulses in the frequency domain ? If not then you might like to check the following blog entry : http://realgonegeek.blogspot.com/2012/12/i-am-performing-fast-fourier-transform.html.

The above input datasets use multiple cycles of sine and cosine waves that have 8 points in their fundamental. This has been chosen because of the ease of creating the wave forms with just the values 0, 0.7071, 1 and their negatives. Here are the base lookup tables that you can use to create the larger input waveforms by simply duplicating them the required number of times :

int sin8[] = {Q28(0.0), Q28(0.70710678118654752440084436210485), Q28(1.0), Q28(0.70710678118654752440084436210485), Q28(0.0), Q28(-0.70710678118654752440084436210485), Q28(-1.0), Q28(-0.70710678118654752440084436210485)};
int cos8[] = {Q28(1.0), Q28(0.70710678118654752440084436210485), Q28(0.0), Q28(-0.70710678118654752440084436210485), Q28(-1.0), Q28(-0.70710678118654752440084436210485), Q28(0.0), Q28(0.70710678118654752440084436210485)};

Note the Qnn macros are used to convert to a Q Format appropriate for a fixed point device, these are not required on a floating point device. If you are not familiar with the Q Format for fixed point numbers then you might like to checkout this blog entry : http://blog.numerix-dsp.com/2015/01/introduction-to-q-format-number-system.html.

If you see anything different to the above then the following points should help find the problem.

Has the code implemented the bit-reversed addressing functionality required in the FFT ?
If so, is this in the right place ?
The rule of thumb is if it is decimation-in-time then bit-reverse in time. If DIF then bit-reverse in the frequency domain.
Has the code scaled the results to allow for the 1/N scaling in the DFT ?
Have you used the forward or reverse DFT function ?
These can be interchanged depending on whether you use e^-jwt or e^jwt to represent the forward or inverse operation. There is no right or wrong answer to which to use. At it is most simple, it depends on whether you and/or the person who wrote the code is an engineer or a mathmatician. If this sounds weird it is a bit like whether to use i or j for imaginary numbers. I'm an engineer so i always refers to current ;-)
Is this code written for a traditional DSP architecture ? If so, does it require that the input and output arrays (including twiddle factors) are correctly aligned on the correct modulo memory address.
Has the twiddle factor (FFT coefficient) look-up-table been implemented correctly ?
This usually takes the form of a 3/4 sine look-up-table.

These are the most common reasons for FFTs not working.
Once you get your FFT working for the simple sine and cosine waveforms then barring fixed-point overflow issues then in general it will work for any input data sequence.

If you have found this solution useful then please do hit the Google (+1) button so that others may be able to find it as well.

Numerix-DSP Libraries : http://www.numerix-dsp.com/eval/

No comments:

Post a Comment