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 &lt;stdio.h&gt;
#include &lt;math.h&gt;

#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 &amp;&amp; (x &amp; (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 &lt; 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 &lt; 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 &lt; 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);
            *even = complex_add(*even, odd_rotated);
            current_rotation = complex_multiply(current_rotation, partial_rotation);
        }
    }
    return TRUE;
}
Bit-reversal permutation bug
Keeping both a normal counter and the &quot;reversed counter&quot; 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 &lt;&lt;&lt;
101 011 &lt;&lt;&lt;
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 &amp; (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.