У меня маленький и слабенький микроконтроллер. У меня тоже нет доступа к сложной библиотеке. Я написал эту итеративную версию БПФ, чтобы, надеюсь, получить лучшую производительность, чем рекурсивная версия, а также просто чтобы узнать, как работает БПФ, и освежить в памяти C.

Я хотел бы знать, понятен ли сам код и можно ли улучшить производительность.

#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;
}

Ошибка перестановки разворота битов

Правильная идея — сохранить как обычный счетчик, так и «обратный счетчик», но эта реализация не совсем верна. Например, это может привести к такой последовательности, как 0, 4, 2, 6, 1, 3, 5, 7 в то время как правильный 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

Есть различные способы сделать это без этой ошибки, но некоторые из них полагаются на инструкции, которые вряд ли будут эффективными (или даже существовать) на микроконтроллерах, такие как начальный нулевой счет и сдвиг по переменной. Я не знаю, что вы можете заставить работать и делать эффективными, так как я не знаю вашего микроконтроллера. Возможно, рекурсивная реализация перестановки с инверсией битов лучше, но, возможно, нет, учитывая часто высокую стоимость рекурсии на микроконтроллерах. Могу только порекомендовать попробовать некоторые варианты, не знаю, как они сработают в вашем конкретном случае.

Также смотрите внизу алгоритмические изменения.

Периодический тест во внутреннем цикле

Под этим я подразумеваю состояние этого if утверждение:

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

Хотя это случай, когда операция остатка с степенью двойки в правой части может быть эффективной, это применимо только тогда, когда компилятор знает он имеет дело со степенью двойки. Гипотетически компилятор мог бы обнаружить этот случай, но на это мало надежды. Компиляция GCC 10 для x64 (я понимаю, что вы не нацелены на это, но на самом деле не в этом суть, дело в способности компилятора рассуждать над набором значений even_to_odd мог бы иметь и как использовать эту информацию), например:

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

Это нехорошо, а на микроконтроллере было бы хуже. Поскольку вы знаете, что even_to_odd является степенью двойки, вы можете записать условие следующим образом: (i & (even_to_odd - 1)) == 0. Если GCC 10 не может сделать это автоматически, маловероятно, что различные компиляторы, зависящие от поставщика для некоторых микроконтроллеров, тоже могут это сделать, поскольку они обычно слабее в этом отношении.

Примечание по сравнительному анализу

Кажется, на моем рабочем столе на 20-40% быстрее, чем написанная мной рекурсивная версия

Обязательно протестируйте его на своем рабочем столе, но когда вариант B работает быстрее, чем вариант A на рабочем столе, то это нет никакой гарантии что это все еще будет актуально для микроконтроллера. Даже между разными десктопами такой гарантии нет, если у них есть процессоры с разной микроархитектурой, например AMD Zen против Intel Skylake. Тестирование кода на рабочем столе может легко заставить вас принять решения, которые вредны для микроконтроллера.

Алгоритмические изменения

БПФ Стокхэма позволяет избежать явной перестановки с обращением битов, устраняя необходимость в его эффективной реализации. С другой стороны, он обращается к памяти другим способом и требует временного рабочего пространства. Может быть, это хорошо для вашего приложения, а может и нет, я не уверен.

Отсутствующий const

Я вижу тебя окропили const почти всюду. Однако вы на самом деле упустили одно место, где это действительно важно: input должен быть const указатель:

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

Использовать restrict ключевое слово, если возможно

С input а также output однотипные, они могут быть псевдонимами. Это может помешать компилятору сгенерировать оптимальный код, поскольку теперь он должен предполагать, что любая запись в значение output может изменить значение в input. Аннотируйте эти указатели с помощью restrict ключевое слово, если оно поддерживается вашим компилятором.

Обратите внимание, что создание input а const указатель не предотвращает сглаживание.

Избегайте разделений

Деление — одна из самых медленных операций на любом ЦП, но особенно на ЦП низкого уровня они могут занимать в десятки раз больше циклов, чем умножение. Вы можете избежать деления при расчете angle написав:

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

Однако обратите внимание, что компилятор может легко оптимизировать деление с помощью констант, поэтому нет необходимости беспокоиться о таких выражениях, как skip /= 2.

Но, возможно, даже лучше:

Рассмотрите возможность использования справочных таблиц sin / cos

cos() а также sin() являются трансцендентными функциями, оценка которых может занять много времени, особенно на слабых микроконтроллерах. С другой стороны, микроконтроллеры обычно имеют очень низкие задержки доступа к памяти. Поэтому имеет смысл предварительно рассчитать все возможные значения sin(angle) а также cos(angle) которые вы можете встретить, и сохраните их в справочной таблице.

В других ответах много хорошей информации, поэтому вот несколько незначительных моментов.

• Использовать <stdbool.h>, если доступно.

Производительность вещей.

• Определенно избегайте оператора мода. В этом случае, учитывая, что ваш делитель — степень 2, предложение @harold отлично. В общем, вы можете использовать один и тот же эффект (каждый N-й элемент) с помощью счетчика:

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

• Вы можете повысить производительность, передав указатели вместо сложных типов, чтобы избежать временных затрат и копирования.

   void complex_init(complex_t* self, double re, double im} {
       *self->real = re; *self->imaj = im; 
    }
    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);

• Что касается sin и cos, помимо исследования таблицы поиска, проверьте, предоставляет ли ваша платформа sincos() в math.h это может быть быстрее.

Стиль вещей

• Остерегаться #define TRUE 1. Здесь все нормально, если вы используете его только в return TRUE;. Но это может соблазнить вас написать позже if (something() == TRUE) который кусает тебя, когда something использует истинное значение, отличное от 1.

• Кстати говоря, гораздо более стандартно возвращаться 0 для отсутствия ошибок и ненулевое значение для сбоев, что позволяет отказаться от использования стандартного или самодельного bool.

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

• Я бы только инициализировал переменные цикла в операторе for, переместите skip а также shuffle выше. Но это личное предпочтение.

    } else
       shuffle += skip;
  

• Никогда не оставляйте такой без скобок