It is my DFT implementation, testing with {0, 1, 2, 3}.

#include <vector>
#include <complex>
#include <numbers>
#include <cmath>
#include <iostream>

std::vector<std::complex<double>> DFT(std::vector<std::complex<double>>&& P)
{
    int n = P.size();
    if (n == 1) {
        return P;
    }

    std::vector<std::complex<double>> Pe(n / 2);
    std::vector<std::complex<double>> Po(n / 2);
    for (int i = 0; i < n / 2; ++i) {
        Pe[i] = P[2 * i];
        Po[i] = P[2 * i + 1];
    }

    auto ye = DFT(std::move(Pe)), yo = DFT(std::move(Po));

    // in place algorithm, use input P to store output
    auto wi = std::complex<double>(1.0, 0.0);
    auto wn = std::complex<double>(std::cos(2.0 * std::numbers::pi_v<double> / (double)n),
                                   std::sin(2.0 * std::numbers::pi_v<double> / (double)n));
    for (int i = 0; i < n / 2; ++i) {
        P[i]         = ye[i] + wi * yo[i];
        P[i + n / 2] = ye[i] - wi * yo[i];
        wi           = wi * wn;
    }

    return P;
}

int main()
{
    int                               N = 4;
    std::vector<std::complex<double>> p_input(N);
    for (int i = 0; i < N; ++i)
    {
        p_input[i] = {(double)i, 0.0};
    }

    auto p_output = DFT(std::move(p_input));

    for (int i = 0; i < p_output.size(); ++i)
    {
        std::cout << p_output[i] << std::endl;
    }
}

The testing result is

(6,0)
(-2,-2)
(-2,0)
(-2,2)

But MATLAB answer is

>> fft(0:1:3)

ans =

   6.0000 + 0.0000i  -2.0000 + 2.0000i  -2.0000 + 0.0000i  -2.0000 - 2.0000i

I have tested longer length of input, and my result always has a reverse order in odd index position.

I don’t know which part is wrong? My reference formula is:

P(x): [p_0, p_1, …, p_{n-1}]

w: [w^0, w^1, …, w^{n-1}]

Pe(x^2): [p_0, p_2, …, p_{n-2}]

Po(x^2): [p_1, p_3, …, p_{n-1}]

ye = [Pe(w^0), Pe(w^2), …, Pe(w^{n-2})]

yo = [Po(w^0), Po(w^2), …, Po(w^{n-2})]

P(w^j) = ye[j] + w^j yo[j]

P(w^{j+n/2}) = ye[j] – w^j yo[j]

y = [P(w^0), P(w^1), …, P(w^{n-1})]