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);
rev ^= change << (__lzcnt64(change) - (64 - stages));
}
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 is the MSVC intrinsic, you could use some preprocessor tricks to find the right intrinsic 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.