Return to Answer

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.

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.

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.

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));
}

__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:

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
}

__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: