# Iterative Radix-2 FFT in C

I have a small and weak microcontroller. I also don't have access to the complex library. I wrote this iterative version of the FFT to hopefully get better performance than the recursive version, and also just to learn how the FFT works, and brush up on C.

I'd like to know if the code itself is understandable, and if there are some performance improvements I can make.

#include <stdio.h>
#include <math.h>

#define FALSE 0
#define TRUE 1
#if !defined(M_PI)
#   define M_PI 3.14159265358979323846
#endif

struct complex {
double real;
double imaj;
};
typedef struct complex complex_t;

complex_t complex_init(const double real, const double imaj) {
complex_t temp;
temp.real = real;
temp.imaj = imaj;
return temp;
}

complex_t complex_add(const complex_t a, const complex_t b)
{
return complex_init(a.real + b.real, a.imaj + b.imaj);
}

complex_t complex_subtract(const complex_t a, const complex_t b) {
return complex_init(a.real - b.real, a.imaj - b.imaj);
}

complex_t complex_multiply(const complex_t a, const complex_t b) {
return complex_init(a.real * b.real - a.imaj * b.imaj, a.real * b.imaj + a.imaj * b.real);
}

_Bool is_power_of_2(const unsigned int x) {
return x != 0 && (x & (x - 1)) == 0;
}

_Bool fft(complex_t* input, complex_t* output, const unsigned int size) {
if (!is_power_of_2(size))
return FALSE;

if (size == 1) {
output[0] = input[0];
return TRUE;
}

const unsigned int half_size = size / 2;

// Initial loop. Do the input shuffle and first butterfly at the same time.
// shuffle is the bit reversed representation of i. If i is 11000, then shuffle is 00011.
for (unsigned int skip = size, i = 0, shuffle = 0; i < half_size; ++i) {
const complex_t even = input[shuffle];
const complex_t odd = input[shuffle + half_size];
output[i * 2] = complex_add(even, odd);
output[i * 2 + 1] = complex_subtract(even, odd);

if (i == 0 || is_power_of_2(i + 1)) {
skip /= 2;
shuffle = skip / 2;
} else
shuffle += skip;
}

// Do the rest of the butterfly operations
for (unsigned int even_to_odd = 2; even_to_odd < size; even_to_odd *= 2) {
const double angle = -M_PI / even_to_odd;
const complex_t partial_rotation = complex_init(cos(angle), sin(angle));
complex_t current_rotation = complex_init(1, 0);
for (unsigned int i = 0, to_even = 0; i < half_size; ++i, ++ to_even) {
if (i % even_to_odd == 0) {
to_even = 2 * i;
current_rotation = complex_init(1, 0);
}
complex_t* even = output + to_even;
complex_t* odd = even + even_to_odd;
const complex_t odd_rotated = complex_multiply(*odd, current_rotation);
*odd = complex_subtract(*even, odd_rotated);
current_rotation = complex_multiply(current_rotation, partial_rotation);
}
}
return TRUE;
}

• And, is it faster than the recursive version? If you don't know, benchmark it. – G. Sliepen Apr 15 at 18:54
• @G.Sliepen Seems to be 20-40% faster on my desktop than the recursive version I wrote. (Larger inputs perform better). The microcontroller is a whole other beast though. 1K memory total. I'm not able to run it on there, as it's not actually mine, but my EE coworker's. I'll have him try it when he has time. – IHerdULiekLambdas Apr 15 at 19:46
• What exactly is a "small weak microcontroller"? You simply shouldn't run code like this on anything smaller than a Cortex M4. If you do, then the root of your problem is that you picked the wrong tool for the project. We can't trust the hardware designer with this - it's the programmer's job to specify the correct MCU. Which is something done after writing a specification. During which you should identify the need for trig or imaginary numbers. Professionals do not pick a random MCU that they found in a packet of corn flakes and then attempt to shoehorn completely unsuitable code into it. – Lundin Apr 16 at 7:59
• Also, this is something you'd usually use a DSP for. There will be pre-made solutions that you can just grab from the shelf. So why re-invent the wheel? – Lundin Apr 16 at 8:01
• @KonradRudolph Microcontroller compilers sort under the category of what C calls "freestanding implementations", where they don't actually need to support a lot of standard headers and still be conforming (see C17 4/6). However, even freestanding compilers must support stdbool.h – Lundin Apr 16 at 8:25

# Bit-reversal permutation bug

Keeping both a normal counter and the "reversed counter" is the right idea, but this implementation isn't quite right. For example, it might result in a sequence such as 0, 4, 2, 6, 1, 3, 5, 7 while the correct one is 0, 4, 2, 6, 1, 5, 3, 7.

A   B
000 000
100 100
010 010
110 110
001 001
011 101 <<<
101 011 <<<
111 111


There are various ways to do this without that bug, but some of them rely on instructions that are unlikely to be efficient (or even exist) on microcontrollers, such as leading zero count and shift-by-variable. I don't know what you can make work, and make efficient, since I don't know your microcontroller. Perhaps a recursive implementation of the bit-reversal permutation is better, but maybe not, given the often high cost of recursion on microcontrollers. I can only recommend trying some options, I don't know how they will work out in your specific case.

Also, see at the bottom for algorithmic changes.

# Periodic test in the inner loop

By that I mean the condition of this if statement:

if (i % even_to_odd == 0) {
to_even = 2 * i;
current_rotation = complex_init(1, 0);
}


While it is the case that a remainder operation with a power of two on the right hand side can be efficient, that only applies when the compiler knows it's dealing with a power of two. Hypothetically, a compiler could detect this case, but there is not much hope for that happening. GCC 10 compiling for x64 (I realize you're not targeting that, but that's not really the point, it's about the compilers ability to reason above the set of values even_to_odd could have and how to use that information) makes this for example:

    mov     eax, ecx
xor     edx, edx
div     ebx
test    edx, edx
jne     .L41


That's not good, and it would be worse on a microcontroller. Since you know that even_to_odd is a power of two, you could write the condition like this: (i & (even_to_odd - 1)) == 0. If GCC 10 cannot do this automatically, it is not likely that the various vendor-specific compilers for some microcontrollers can do it either, as they are usually weaker in that regard.

# Benchmarking note

Seems to be 20-40% faster on my desktop than the recursive version I wrote

By all means test it on your desktop, but when option B is faster than option A on a desktop, that is no guarantee at all that that will still be true on the microcontroller. Even between different desktops there is no such guarantee, if they have processors with a different micro-architecture, for example AMD Zen vs Intel Skylake. Benchmarking the code on your desktop can easily trick you into making decisions that are bad for the microcontroller.

# Algorithmic changes

The Stockham FFT avoids an explicit bit-reversal permutation, removing the need to implement it efficiently. On the other hand, it accesses memory in a different way and needs temporary work-space. Maybe it's good for your application, or maybe not, I'm not sure.

• Thanks for finding the bug. Fixed it. Added some time (but actually gives the right results now). Still outpacing the recursive version. The inner loop test change is great, though. I know my desktop isn't a good benchmark for the microcontroller, but I don't actually have access to it or the tools to build for it, since it's actually my coworker's. – IHerdULiekLambdas Apr 15 at 21:51

# Missing const

I see you sprinkled const almost everywhere. However, you actually missed the one spot where it actually matters most: input should be a const pointer:

_Bool fft(const complex_t* input, complex_t* output, const unsigned int size) {
...
}


# Use the restrict keyword if possible

Since input and output are of the same type, they can alias. This might prevent the compiler from generating optimal code, as it now must assume that any write to a value of output might change a value in input. Annotate these pointers with the restrict keyword if supported by your compiler.

Note that making input a const pointer doesn't prevent aliasing.

# Avoid divisions

Divisions are one of the slowest operations on any CPU, but especially on low end CPUs they might take tens of times more cycles than a multiply. You can avoid the division when calculating angle by writing:

double angle = -M_PI;
for (unsigned int even_to_odd = 2; even_to_odd < size; even_to_odd *= 2) {
angle *= 0.5;
...
}


However, note that the compiler can easily optimize divisions by constants, so there is no need to worry about expressions like skip /= 2.

But perhaps even better:

# Consider using sin/cos look-up tables

cos() and sin() are transcedental functions that can take a long time to evaluate, especially on weak microcontrollers. On the other hand, microcontrollers usually have very low memory access latencies. So it might make sense to pre-calculate all the possible values of sin(angle) and cos(angle) that you might encounter, and store them in a look-up table.

• @IHerdULiekLambdas: Most mainstream C++ implementations support __restrict as an extension. (MSVC, and GNU-compatible including clang and ICC). Not surprised that avoiding division helped a lot on a modern x86-64 CPU; it's kind of the only thing that's still slow for throughput. Although LUTs for sin/cos can still help on modern x86 if they don't have to be huge, otherwise you're just costing as much in cache misses as you save in FP ALU work, unless the access pattern is bad. (Which wouldn't be the case on a microcontroller). Also, LUTs may hurt a compiler's ability to auto-vectorize. – Peter Cordes Apr 16 at 5:06
• @IHerdULiekLambdas: Floating point division vs floating point multiplication compares things for modern x86 CPUs. – Peter Cordes Apr 16 at 5:08
• "but especially on low end CPUs they might take tens of times more cycles than a multiply" Well, code like this shouldn't be executed on such a low-end part. That's just bad design. You can get a Cortex M4 for a dollar, so why would you pick some dysfunctional 1990s 8-bitter architecture when your project needs floating point and imaginary numbers. That's like showing up at the start of a Formula 1 race in a rusty tractor. After which the problem isn't to adapt your driving style to suit the tractor, but that you showed up in a tractor in the first place. – Lundin Apr 16 at 8:06
• @Lundin The answer already addresses that: in this context, const* and *restrict are not the same thing, since the compiler can’t prove that they aren’t aliasing (and in fact they could be). – Konrad Rudolph Apr 16 at 8:22
• @Lundin A Cortex M4 doesn't have hardware support for 64-bit floating point numbers. But regardless of whether you are using an 8-bitter or the fastest CPU on earth, optimizing the code and avoiding expensive operations is a worthy goal. Also, if the code can be made to run fast enough for its intended purpose on an 8-bit CPU, why not? – G. Sliepen Apr 16 at 16:10

Lots of good info in the other answers, so here are a few minor points.

• Use <stdbool.h>, if available.

Performance things.

• Definitely avoid the mod operator. In this case, given your divisor is a power of 2, @harold's suggestion is great. In general, you can the same effect (every N'th item) with a counter:

 if (--remaining > 0) { remaining = even_to_odd; to_even = 2*i; ... }

• You may get a performance speedup by passing pointers instead of complex types, to avoid temps and copying
   static inline void complex_init(complex_t* self, double re, double im} {
*self->real = re; *self->imaj = im;
}
static inline void complex_add(complex_t* result, const complex_t* a, const complex_t* b){
complex_init(result, a->real + b->real, a->imaj + b->imaj);
}
//usage:
complex_t partial_rotation;
complex_init(&partial_rotation, cos(angle), sin(angle));
...


• Regarding sin and cos, in addition to investigating the lookup table, check if your platform provides sincos() in math.h, it may be faster.

Style things

• Beware of #define TRUE 1. It's OK here, where you are only using it in return TRUE;. But it can tempt you into later writing if (something() == TRUE) which bites you when something uses a truthy value other than 1.

• Speaking of which, it's far more standard to return 0 for no errors, and non-zero for failures, which lets you get away from using either standard or homemade bool.

for (unsigned int skip = size, i = 0, shuffle = 0; i < half_size; ++i) {

• I'd only initialize the loop variables in the for statement, move skip and shuffle above. But that's personal preference.

  } else
shuffle += skip;

• Never leave an unbracketed else like this

• You want a trivial function like complex_init to just inline, so declare it inline. Setting up the args is going to cost as much space as just inlining the assignments. Possibly not for complex_add on a system without hardware FP (so each double + costs a function call anyway, perhaps with pointers to memory), but otherwise you want that to inline as well. – Peter Cordes Apr 17 at 1:59