Xjyuk

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)
 ((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?

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 borrows 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.

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