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

#include <stdio.h>
#include <vector>
#include <iostream>
#include <complex>
#include <cmath>
#include <algorithm>

#define N 1048576
#define PI 3.14159265358979323846

/*
Creating the table of all N'th roots of unity.
We use notation omega_k = e^(-2 pi i / n).
*/
template<typename U>
std::vector< std::complex<U> > rootsOfUnityCalculator() {
    std::vector< std::complex<U> > table;

    for (size_t k = 0; k < N; k++) {
        std::complex<U> kthRootOfUnity(std::cos(-2.0 * PI * k / N), std::sin(-2.0 * PI * k / N));
        table.push_back(kthRootOfUnity);
    }

    return table;
}

/*
Fast Fourier transform, T is the precision level, so float or double.
table is a look up table of the roots of unity. Overwrites the input.
For now only works for N a power of 2.
*/
template<typename T>
void FFT(std::complex<T>* input, const std::vector< std::complex<T> >& table, size_t n) {

    if (n % 2 == 0) {
        // Split up the input in even and odd components
        std::complex<T>* evenComponents = new std::complex<T>[n/2];
        std::complex<T>* oddComponents = new std::complex<T>[n/2];

        for (size_t k = 0; k < n/2; k++) {
            evenComponents[k] = input[2 * k];
            oddComponents[k] = input[2 * k + 1];
        }

        // Transform the even and odd input
        FFT(evenComponents, table, n/2);
        FFT(oddComponents, table, n/2);

        // Use the algorithm from Danielson and Lanczos
        for (size_t k = 0; k < n/2; k++) {
            std::complex<T> plusMinus = table[N / n * k] * oddComponents[k]; // omega_n^k = (omega_N^(N/n))^k = omega_N^(Nk/n)
            input[k] = evenComponents[k] + plusMinus;
            input[k + n/2] = evenComponents[k] - plusMinus;
        }

        delete[] evenComponents;
        delete[] oddComponents;

    } else {
        // The Fourier transform on one element does not do anything, so
        // nothing needed here.
    }
}

int main() {
    std::complex<double>* input = new std::complex<double>[N];

    for (size_t k = 0; k < N; k++) {
        input[k] = k;
    }

    const std::vector< std::complex<double> > table = rootsOfUnityCalculator<double>();

    // Overwrites the input with its Fourier transform
    FFT<double>(input, table, N);

    delete[] input;

    return 0;
}

Избегайте чрезмерного использования памяти

На каждом шаге рекурсии вы выделяете два массива, размер которых вместе равен входному. Это означает, что ваш алгоритм использует $ mathcal {O} (N log N) $ Космос. Должна быть возможность переписать ваш код, чтобы вы использовали только $ mathcal {O} (N) $ пространство без существенного изменения алгоритма. Но это также подводит меня к следующему:

Алгоритмы на месте и не на месте

Хотя сигнатура функции делает вид, будто вы выполняете преобразование БПФ на месте, на самом деле вы выполняете преобразование вне места (используя дополнительную память), а затем копируете результат обратно во входной массив. Однако есть и настоящие алгоритмы БПФ на месте которые, кроме нескольких переменных, не нуждаются в дополнительной памяти и по своей природе перезаписывают входной массив. Это означает, что они используют только $ mathcal {O} (1) $ объем памяти. Недостатком алгоритмов на месте является то, что преобразованные данные не в том порядке, который вы ожидаете; вам, возможно, придется использовать бит-инвертированный индекс для доступа к данным. Однако для некоторых приложений это вообще не имеет значения, например, если вы выполняете свертки с использованием БПФ.

Если вы все же используете алгоритм, не относящийся к месту, я предлагаю вам сделать это явным и либо вернуть std::vector с результатами или добавьте к функции параметр, который сообщает ей, где хранить результаты.

Предпочитайте заголовки C ++ (например, <cstdio>) в заголовки совместимости с C (<stdio.h>). Заголовки C ++ определяют идентификаторы в std пространство имен там, где мы хотим.

Кажется, что этот заголовок здесь даже не используется, поэтому мы можем полностью его опустить.

Мы используем только PI как double, поэтому мы можем объявить его как константу, а не как макрос:

constexpr double pi = 3.14159265358979323846;

Мы можем избежать голых new[] и delete[] с помощью std::vector. Если мы используем std::vector::reserve перед добавлением каких-либо элементов не должно быть никакого влияния на производительность, только повышенная безопасность кода.

Использование вектора также избавляет от необходимости явно передавать размер массива.

Полезный метод определения того, является ли число степенью двойки, заключается в том, что побитовое и числа, а его арифметический предшественник равняется нулю для степеней двойки:

constexpr bool is_power_of_two(unsigned n)
{
    return n && (n & (n-1) == 0);
}

Одно значение N наверное, вообще не очень полезно. Мы могли бы сделать его параметром шаблона или простым аргументом функции, чтобы сделать наш код полезным для других размеров ввода.

Модифицированная версия

Я раскололся FFT() в простую общедоступную функцию, которая обеспечивает правильность аргументов и внутреннюю рекурсивную реализацию. Я использовал requires чтобы ограничить сложный тип, который мы принимаем (целочисленные типы здесь бесполезны), но это можно опустить для простого C ++ 17.

#include <complex>
#include <vector>

namespace internal
{
    /*
      Creating the table of all N'th roots of unity.
      We use notation ω_i = e^(-2πi/n).
    */
    template<typename U>
    constexpr auto rootsOfUnity(size_t n)
    {
        constexpr double pi = 3.14159265358979323846;
        std::vector<std::complex<U>> table;
        table.reserve(n);

        for (size_t i = 0;  i < n;  ++i) {
            table.emplace_back(std::cos(-2.0 * pi * i / n), std::sin(-2.0 * pi * i / n));
        }

        return table;
    }

    template<typename V>        // V is vector of complex
    void FFT(V& input, const V& table, size_t total_size)
    {
        if (input.size() <= 1) {
            // The Fourier transform on one element does not do
            // anything, so nothing needed here.
            return;
        }

        const auto n2 = input.size() / 2;

        // Split up the input in even and odd components
        V evenComponents;  evenComponents.reserve(n2);
        V oddComponents;   oddComponents.reserve(n2);

        for (size_t i = 0;  i < input.size();  ) {
            evenComponents.push_back(input[i++]);
            oddComponents.push_back(input[i++]);
        }

        // Transform the even and odd input
        FFT(evenComponents, table, total_size);
        FFT(oddComponents, table, total_size);

        // Use the algorithm from Danielson and Lanczos
        for (size_t i = 0;  i < n2;  ++i) {
             // ω_n^i = (ω_N^(N/n))^i = ω_N^(Nk/n)
            auto const plusMinus = table[total_size / n2 * i] * oddComponents[i];
            input[i] = evenComponents[i] + plusMinus;
            input[i + n2] = evenComponents[i] - plusMinus;
        }
    }
}

/*
  Fast Fourier transform
  T is the precision level: float, double, long double.
  Overwrites the input.
  Only works for input sizes that are a power of 2.
*/
template<typename T>
    requires std::is_floating_point_v<T>
void FFT(std::vector<std::complex<T>>& input)
{
    if (input.size() == 0) {
        throw new std::invalid_argument("input is empty");
    }
    if (input.size() & (input.size() - 1)) {
        throw new std::invalid_argument("input length is not a power of two");
    }

    auto const table = internal::rootsOfUnity<T>(input.size());
    internal::FFT(input, table, input.size());
}

// Test program
#include <iostream>

int main()
{
    static constexpr size_t n = 1048576;
    std::vector<std::complex<double>> input;
    for (size_t i = 0;  i < n;  ++i) {
        input.emplace_back(i);
    }

    // Overwrite input with its Fourier transform
    FFT<double>(input);

    // use the result, to avoid optimising out
    std::cout << input[0] << 'n';
}

Использовать unique_ptr вместо new[]/delete[]

std::complex<double>* input = new std::complex<double>[N];

...

delete[] input;

Эту строку можно заменить соответствующим использованием std::unique_ptr:

auto input_storage = std::make_unique<std::complex<double>[]>(N);  // #include <memory>
std::complex<double>* input = input_storage.get();

Память автоматически восстанавливается один раз input_storage выходит за рамки.