15
\$\begingroup\$

Note: If you don't know much about Fourier transform algorithms, a simple review of whether I am doing anything inefficient with C++ in general would be appreciated.

I've been working on implementing an efficient Radix2 Fast Fourier Transform in C++ and I seem to have hit a roadblock. I have optimized it in every possible way I can think of and it is very fast, but when comparing it to the Numpy FFT in Python it is still significantly slower. Note that my FFT is not done in-place, but neither is the Python implementation so I should be able to achieve at least the same efficiency as Numpy.

I've already taken advantage of symmetry when the input is a real signal which allows the use of an N/2 point FFT for an N-length real signal, and I also pre-compute all of the twiddle factors and optimize the twiddle factor calculation so redundant twiddle factors are not re-calculated.

The code:

#include <cassert>
#include <complex>
#include <vector>

// To demonstrate runtime
#include <chrono>
#include <iostream>

static std::vector<std::complex<double>> FFTRadix2(const std::vector<std::complex<double>>& x, const std::vector<std::complex<double>>& W);
static bool IsPowerOf2(size_t x);
static size_t ReverseBits(const size_t x, const size_t n);

std::vector<std::complex<double>> FFT(const std::vector<double>& x)
{
    size_t N = x.size();

    // Radix2 FFT requires length of the input signal to be a power of 2
    // TODO: Implement other algorithms for when N is not a power of 2
    assert(IsPowerOf2(N));

    // Taking advantage of symmetry the FFT of a real signal can be computed
    // using a single N/2-point complex FFT. Split the input signal into its
    // even and odd components and load the data into a single complex vector.
    std::vector<std::complex<double>> x_p(N / 2);
    for (size_t n = 0; n < N / 2; ++n)
    {
        // x_p[n] = x[2n] + jx[2n + 1]
        x_p[n] = std::complex<double>(x[2 * n], x[2 * n + 1]);
    }

    // Pre-calculate twiddle factors
    std::vector<std::complex<double>> W(N / 2);
    std::vector<std::complex<double>> W_p(N / 4);
    for (size_t k = 0; k < N / 2; ++k)
    {
        W[k] = std::polar(1.0, -2 * M_PI * k / N);

        // The N/2-point complex DFT uses only the even twiddle factors
        if (k % 2 == 0)
        {
            W_p[k / 2] = W[k];
        }
    }

    // Perform the N/2-point complex FFT
    std::vector<std::complex<double>> X_p = FFTRadix2(x_p, W_p);

    // Extract the N-point FFT of the real signal from the results 
    std::vector<std::complex<double>> X(N);
    X[0] = X_p[0].real() + X_p[0].imag();
    for (size_t k = 1; k < N / 2; ++k)
    {
        // Extract the FFT of the even components
        auto A = std::complex<double>(
            (X_p[k].real() + X_p[N / 2 - k].real()) / 2,
            (X_p[k].imag() - X_p[N / 2 - k].imag()) / 2);

        // Extract the FFT of the odd components
        auto B = std::complex<double>(
            (X_p[N / 2 - k].imag() + X_p[k].imag()) / 2,
            (X_p[N / 2 - k].real() - X_p[k].real()) / 2);

        // Sum the results and take advantage of symmetry
        X[k] = A + W[k] * B;
        X[k + N / 2] = A - W[k] * B;
    }

    return X;
}

std::vector<std::complex<double>> FFT(const std::vector<std::complex<double>>& x)
{
    size_t N = x.size();

    // Radix2 FFT requires length of the input signal to be a power of 2
    // TODO: Implement other algorithms for when N is not a power of 2
    assert(IsPowerOf2(N));

    // Pre-calculate twiddle factors
    std::vector<std::complex<double>> W(N / 2);
    for (size_t k = 0; k < N / 2; ++k)
    {
        W[k] = std::polar(1.0, -2 * M_PI * k / N);
    }

    return FFTRadix2(x, W);
}

static std::vector<std::complex<double>> FFTRadix2(const std::vector<std::complex<double>>& x, const std::vector<std::complex<double>>& W)
{
    size_t N = x.size();

    // Radix2 FFT requires length of the input signal to be a power of 2
    assert(IsPowerOf2(N));

    // Calculate how many stages the FFT must compute
    size_t stages = static_cast<size_t>(log2(N));

    // Pre-load the output vector with the input data using bit-reversed indexes
    std::vector<std::complex<double>> X(N);
    for (size_t n = 0; n < N; ++n)
    {
        X[n] = x[ReverseBits(n, stages)];
    }

    // Calculate the FFT one stage at a time and sum the results
    for (size_t stage = 1; stage <= stages; ++stage)
    {
        size_t N_stage = static_cast<size_t>(std::pow(2, stage));
        size_t W_offset = static_cast<size_t>(std::pow(2, stages - stage));
        for (size_t k = 0; k < N; k += N_stage)
        {
            for (size_t n = 0; n < N_stage / 2; ++n)
            {
                auto tmp = X[k + n];
                X[k + n] = tmp + W[n * W_offset] * X[k + n + N_stage / 2];
                X[k + n + N_stage / 2] = tmp - W[n * W_offset] * X[k + n + N_stage / 2];
            }
        }
    }

    return X;
}

// Returns true if x is a power of 2
static bool IsPowerOf2(size_t x)
{
    return x && (!(x & (x - 1)));
}

// Given x composed of n bits, returns x with the bits reversed
static size_t ReverseBits(const size_t x, const size_t n)
{
    size_t xReversed = 0;
    for (size_t i = 0; i < n; ++i)
    {
        xReversed = (xReversed << 1U) | ((x >> i) & 1U);
    }

    return xReversed;
}

int main()
{
    size_t N = 16777216;
    std::vector<double> x(N);

    int f_s = 8000;
    double t_s = 1.0 / f_s;

    for (size_t n = 0; n < N; ++n)
    {
        x[n] = std::sin(2 * M_PI * 1000 * n * t_s)
            + 0.5 * std::sin(2 * M_PI * 2000 * n * t_s + 3 * M_PI / 4);
    }

    auto start = std::chrono::high_resolution_clock::now();
    auto X = FFT(x);
    auto stop = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::microseconds>(stop - start);

    std::cout << duration.count() << std::endl;
}

Output (running a few times and averaging):

3671677

This was compiled in Visual Studio 2019 in Release mode with the /O2, /Oi and /Ot optimization compiler flags to try and squeeze as much speed as possible out of it.


A comparable snippet of Python code that uses the Numpy FFT is shown below:

import numpy as np
import datetime

N = 16777216
f_s = 8000.0
t_s = 1/f_s

t = np.arange(0, N*t_s, t_s)
y = np.sin(2*np.pi*1000*t) + 0.5*np.sin(2*np.pi*2000*t + 3*np.pi/4)

start = datetime.datetime.now()
Y = np.fft.fft(y)
stop = datetime.datetime.now()
duration = stop - start

print(duration.total_seconds()*1e6)

Output (running a few times and averaging):

2100411.0


As you can see, the Python implementation is still faster by about 43%, but I can't think of any ways my implementation can be improved.

From what I understand, the Numpy version is actually implemented with C code underneath so I'm not terribly disappointed in the performance of my own code, but it still leaves me wondering what I am missing that I could still do better?

\$\endgroup\$
  • 1
    \$\begingroup\$ Can you give a hint about the used C++ standard (11, 14, 17) Can you use additional libraries such as gsl or abseil? \$\endgroup\$ – miscco Aug 17 at 17:26
  • 1
    \$\begingroup\$ @miscco, C++11 is currently used but C++17 is available. Additional libraries are not out of the question. If they are available through the Conan package manager that would be a huge plus. Also it is important that it stays cross platform compatible so OS specific libraries are a no go. \$\endgroup\$ – tjwrona1992 Aug 17 at 17:32
  • \$\begingroup\$ Radix 2 is a slow algorithm by today's standards \$\endgroup\$ – Scott Seidman Aug 18 at 20:45
  • \$\begingroup\$ You're only letting MSVC's auto-vectorize with SSE2? I'd at least try letting it use AVX2 and maybe a fast-math option. Or preferably a compiler like gcc or clang with -O3 -march=native -ffast-math. Or if you have ICC it's well known for good auto-vectorization. NumPy might well be manually vectorized for SSE2 and AVX / AVX2, with intrinsics like _mm_shuffle_ps() for SIMD vectors. And _mm_shuffle_epi8 for bit-reversal using a lookup table for 4-bit chunks. Plain portable ISO C++ can't represent/expose a lot of useful things that modern CPUs can do. \$\endgroup\$ – Peter Cordes Aug 18 at 21:13
  • \$\begingroup\$ @ScottSeidman, can you provide the names of any faster algorithms? As far as I am aware Radix 2 is essentially the fastest you can get but comes with a restriction that your input length is a power of 2. \$\endgroup\$ – tjwrona1992 Aug 18 at 21:19
18
\$\begingroup\$

Putting this through the built-in profiler reveals some hot spots. Perhaps surprisingly: ReverseBits. It's not the biggest thing in the list, but it is significant while it shouldn't be.

You could use one of the many alternate ways to implement ReverseBits, or the sequence of bit-reversed indexes (which does not require reversing all the indexes), or the overall bit-reversal permutation (which does not require bit reversals).

For example here is a way to compute the sequence of bit-reversed indexes without explicitly reversing any index:

for (size_t n = 0, rev = 0; n < N; ++n)
{
    X[n] = x[rev];
    size_t change = n ^ (n + 1);
#if _WIN64
    rev ^= change << (__lzcnt64(change) - (64 - stages));
#else
    rev ^= change << (__lzcnt(change) - (32 - stages));
#endif
}

On my PC, that reduces the time from around 2.8 million microseconds to 2.3 million microseconds.

This trick works by using that the XOR between adjacent indexes is a mask of ones up to and including the least significant zero (the +1 carries through the least significant set bits and into that least significant zero), which has a form that can be reversed by just shifting it. The reversed mask is then the XOR between adjacent reversed indexes, so applying it to the current reversed index with XOR increments it.

__lzcnt64 and _WIN64 are for MSVC, you could use more preprocessor tricks to find the right intrinsic and bitness-detection for the current compiler. Leading zero count can be avoided by using std::bitset and its count method:

size_t change = n ^ (n + 1);
std::bitset<64> bits(~change);
rev ^= change << (bits.count() - (64 - stages));

count is recognized by GCC and Clang as an intrinsic for popcnt, but it seems not by MSVC, so it is not reliable for high performance scenarios.

Secondly, there is a repeated expression: W[n * W_offset] * X[k + n + N_stage / 2]. The compiler is often relied on to remove such duplication, but here it didn't happen. Factoring that out reduced the time to under 2 million microseconds.

Computing the twiddle factors takes a bit more time than it needs to. They are powers of the first non-trivial twiddle factor, and could be computed iteratively that way. This suffers from some build-up of inaccuracy, which could be improved by periodically resetting to the proper value computed by std::polar. For example,

auto twiddle_step = std::polar(1.0, -2.0 * M_PI / N);
auto twiddle_current = std::polar(1.0, 0.0);
for (size_t k = 0; k < N / 2; ++k)
{
    if ((k & 0xFFF) == 0)
        twiddle_current = std::polar(1.0, -2.0 * M_PI * k / N);
    W[k] = twiddle_current;
    twiddle_current *= twiddle_step;

    // The N/2-point complex DFT uses only the even twiddle factors
    if (k % 2 == 0)
    {
        W_p[k / 2] = W[k];
    }
}

On my PC that reduces the time from hovering around 1.95 million µs to around 1.85 million µs, not a huge difference but easily measurable.

More advanced: use SSE3 for the main calculation, for example (not well tested, but seems to work so far)

__m128d w_real = _mm_set1_pd(W[n * W_offset].real());
__m128d w_imag = _mm_set1_pd(W[n * W_offset].imag());
__m128d z = _mm_loadu_pd(reinterpret_cast<double*>(&X[k + n + N_stage / 2]));
__m128d z_rev = _mm_shuffle_pd(z, z, 1);
__m128d t = _mm_addsub_pd(_mm_mul_pd(w_real, z), _mm_mul_pd(w_imag, z_rev));

__m128d x = _mm_loadu_pd(reinterpret_cast<double*>(&X[k + n]));
__m128d t1 = _mm_add_pd(x, t);
__m128d t2 = _mm_sub_pd(x, t);
_mm_storeu_pd(reinterpret_cast<double*>(&X[k + n]), t1);
_mm_storeu_pd(reinterpret_cast<double*>(&X[k + n + N_stage / 2]), t2);

That takes it from 1.85 million µs down to around 1.6 million µs on my PC.


Using a different algorithm, Stockham algorithm the version from List-8 and some miscellaneous things, the time goes down to 0.9 million µs. It's a huge win already and this is not the best version of the algorithm. The linked website has faster versions with fancier tricks and SIMD too, so it's there if you want it. As a bonus, no bit reversing is used at all, so no need for a compiler-specific intrinsic.

The real work happens here: (taken from the linked website)

void fft0(int n, int s, bool eo, complex_t* x, complex_t* y)
// n  : sequence length
// s  : stride
// eo : x is output if eo == 0, y is output if eo == 1
// x  : input sequence(or output sequence if eo == 0)
// y  : work area(or output sequence if eo == 1)
{
    const int m = n / 2;
    const double theta0 = 2 * M_PI / n;

    if (n == 2) {
        complex_t* z = eo ? y : x;
        for (int q = 0; q < s; q++) {
            const complex_t a = x[q + 0];
            const complex_t b = x[q + s];
            z[q + 0] = a + b;
            z[q + s] = a - b;
        }
    }
    else if (n >= 4) {
        for (int p = 0; p < m; p++) {
            const complex_t wp = complex_t(cos(p*theta0), -sin(p*theta0));
            for (int q = 0; q < s; q++) {
                const complex_t a = x[q + s * (p + 0)];
                const complex_t b = x[q + s * (p + m)];
                y[q + s * (2 * p + 0)] = a + b;
                y[q + s * (2 * p + 1)] = (a - b) * wp;
            }
        }
        fft0(n / 2, 2 * s, !eo, y, x);
    }
}

void fft(int n, complex_t* x) // Fourier transform
// n : sequence length
// x : input/output sequence
{
    complex_t* y = new complex_t[n];
    fft0(n, 1, 0, x, y);
    delete[] y;
    // scaling removed because OP doesn't do it either
    //for (int k = 0; k < n; k++) x[k] /= n;
}

And here is that wrapper to do a Real FFT with a Complex FFT with half the number of points,

std::vector<std::complex<double>> FFT2(const std::vector<double>& x)
{
    size_t N = x.size();

    // Radix2 FFT requires length of the input signal to be a power of 2
    // TODO: Implement other algorithms for when N is not a power of 2
    assert(IsPowerOf2(N));

    // Taking advantage of symmetry the FFT of a real signal can be computed
    // using a single N/2-point complex FFT. Split the input signal into its
    // even and odd components and load the data into a single complex vector.
    std::vector<std::complex<double>> x_p(N / 2);
    std::copy(x.data(), x.data() + x.size(), reinterpret_cast<double*>(x_p.data()));

    fft(N / 2, x_p.data());

    // Extract the N-point FFT of the real signal from the results 
    std::vector<std::complex<double>> X(N);
    X[0] = x_p[0].real() + x_p[0].imag();
    auto twiddle_step = std::polar(1.0, -2.0 * M_PI / N);
    auto twiddle_current = twiddle_step;
    for (size_t k = 1; k < N / 2; ++k)
    {
        auto Wk = twiddle_current;
        // Extract the FFT of the even components
        auto A = std::complex<double>(
            (x_p[k].real() + x_p[N / 2 - k].real()) / 2,
            (x_p[k].imag() - x_p[N / 2 - k].imag()) / 2);

        // Extract the FFT of the odd components
        auto B = std::complex<double>(
            (x_p[N / 2 - k].imag() + x_p[k].imag()) / 2,
            (x_p[N / 2 - k].real() - x_p[k].real()) / 2);

        // Sum the results and take advantage of symmetry
        X[k] = A + Wk * B;
        X[k + N / 2] = A - Wk * B;

        twiddle_current *= twiddle_step;
    }

    return X;
}

Using std::copy was faster than a manual loop, and not storing the twiddles was also faster. Of course I used the fast twiddle factor generation scheme (without resets this time as per the comments, of course that's easy to put back in). Avoiding the copy altogether would obviously be better, but then the input data will be turned into its FFT instead of leaving it read-only, it's not a drop-in replacement.

Extracting the Real FFT takes a significant portion of the total time by the way.

\$\endgroup\$
  • \$\begingroup\$ Is there a way that doesn't involve preprocessor tricks to make this cross platform compatible? I'd prefer to avoid using the preprocessor although it's not completely out of the question. \$\endgroup\$ – tjwrona1992 Aug 17 at 17:46
  • 1
    \$\begingroup\$ @tjwrona1992 I added an alternative. You could try one of the other permutation algorithms to avoid this problem altogether. \$\endgroup\$ – harold Aug 17 at 17:59
  • \$\begingroup\$ Thanks this all looks like very useful advice! I'll try to take advantage of some of it and let you know how it goes. Out of curiosity, how did you know the compiler didn't factor out the duplicated code? \$\endgroup\$ – tjwrona1992 Aug 17 at 18:00
  • 2
    \$\begingroup\$ @tjwrona1992 by proof-reading the assembly code, painful as it was \$\endgroup\$ – harold Aug 17 at 18:05
  • 2
    \$\begingroup\$ @tjwrona1992 yes good point, would you accept working with uint32_t? It can be made size-adaptive but at the cost of more verbosity.. \$\endgroup\$ – harold Aug 18 at 1:16
7
\$\begingroup\$
  • These lines:

    size_t N_stage = static_cast<size_t>(std::pow(2, stage));
    size_t W_offset = static_cast<size_t>(std::pow(2, stages - stage));
    

    should not use floating-point math because they can be inaccurate. Instead use pure integer arithmetic:

    size_t N_stage = static_cast<size_t>(1) << stage;
    size_t W_offset = static_cast<size_t>(1) << (stages - stage);
    
  • Each iteration of your loop for (size_t stage = 1; stage <= stages; ++stage) will linearly traverse the entire vector once. But this is not optimal if you consider the memory hierarchy. It would be a large change in your code, but you could rework the memory access pattern so that you transform increasingly larger blocks. This technique is known as a cache-oblivious algorithm.

  • As a matter of personal taste, I would do using std::vector and using std::complex because they are referred to so many times in the code.

  • Size 2 and size 4 DFTs have trivial integer twiddle factors (without irrationals), so you could special-case them to save some multiplications. You can use the formulas above to special-case your outer loop when stage = 1 (length-2 DFT) and stage = 2 (length-4 DFT). I have a working example on another page.

    The DFT of the length-2 complex vector [x0, x1] looks like this:

    X0 = x0 + x1
    X1 = x0 - x1
    

    The DFT of the length-4 complex vector [x0, x1, x2, x3] looks like this:

    X0 = x0 + x1 + x2 + x3
    X1 = x0 - i*x1 - x2 + i*x3
    X2 = x0 - x1 + x2 - x3
    X3 = x0 + i*x1 - x2 - x*x3
    

    If a complex number is represented as a pair of real numbers in rectangular form, then multiplication by i is just a matter of swapping the real/imaginary parts and negating the correct part - so no multiplication or addition is needed for this operation.

  • The famous FFTW library has a bunch of speedup techniques, and there are articles you can find online that talk about how they work.

  • Overall your FFTRadix2() looks quite similar to my FftComplex.cpp Fft::transformRadix2().

\$\endgroup\$
  • \$\begingroup\$ Thanks @Nayuki! There is a lot here so I will reply with a comment per bullet. I have done some refactoring already which took care of bullet 1. \$\endgroup\$ – tjwrona1992 Aug 18 at 21:31
  • \$\begingroup\$ bullet 2 is interesting and I'll have to do more research on that. I'm still fairly inexperienced when it comes to DSP algorithms, but I'm working my way through a book and learning a lot. \$\endgroup\$ – tjwrona1992 Aug 18 at 21:32
  • \$\begingroup\$ bullet 3 I agree, I will likely add using statements because it is a bit cumbersome haha \$\endgroup\$ – tjwrona1992 Aug 18 at 21:32
  • \$\begingroup\$ bullet 4, could you possibly elaborate on this? I'm not sure I fully understand what you mean by this but it sounds interesting \$\endgroup\$ – tjwrona1992 Aug 18 at 21:32
  • 2
    \$\begingroup\$ bullet 5, I actually stumbled across that library when I was part-way through writing mine and used it as a reference. It was VERY helpful! Thank you! :) \$\endgroup\$ – tjwrona1992 Aug 18 at 21:33
4
\$\begingroup\$

You should be able to work inplace by utilizing the fact that std::complex has a predefined layout. Therefore you can actually cast between an array of double and an array of (half as many) complex numbers.

std::vector<std::complex<double>> a(10);
double *b = reinterpret_cast<double *>(a.data());

EDIT:

To be more clear I would write

span<std::complex<double>> x_p(reinterpret_cast<std::complex<double>*>(x.data()), x.size() / 2);

This works in both ways. To enable safe and modern features you should use a span object. Unfortunately std::span is only available in C++20 so you should either write your own (which is a nice exercise) or have a look at abseil::span or gsl::span.

The code to implement those is rather minimal. With that you can remove two copies from your code

\$\endgroup\$
  • \$\begingroup\$ Are you saying replace vector with span? I haven't used span before. I'll look into it! \$\endgroup\$ – tjwrona1992 Aug 17 at 18:01
  • 1
    \$\begingroup\$ No, A span is non-owning. But rather than copying the data into x_p you should create a span that reinterprests the data in x as an array of std::complex \$\endgroup\$ – miscco Aug 17 at 18:22
  • \$\begingroup\$ Ahh I see, maybe I will implement another one that does it in place like that. But I would also like to have the option to do it out of place as well \$\endgroup\$ – tjwrona1992 Aug 17 at 23:09
  • \$\begingroup\$ Do any compilers already support C++20? I didn't think it was officially released yet \$\endgroup\$ – tjwrona1992 Aug 18 at 1:34
  • 2
    \$\begingroup\$ @PeterCordes there is the array-oriented access thing \$\endgroup\$ – harold Aug 19 at 1:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.