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.