Iterative Radix-2 FFT in C

Question

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);
            *even = complex_add(*even, odd_rotated);
            current_rotation = complex_multiply(current_rotation, partial_rotation);
        }
    }
    return TRUE;
}

Answer 1

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.

Answer 2

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.

Answer 3

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));
   ...
   complex_add(even, even, &odd_rotated);

• 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