Custom DFT implementation return reverse order of odd index answer

Question

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

Answer 1

You're just using a different nth root of unity than their implementation, namely

cos(2*pi / n) + sin(2*pi / n) * i

instead of

cos(2*pi / n) - sin(2*pi / n) * i

As far as FFT is concerned, which one you use is irrelevant as long as the inverse FFT is consistent with it. I'm not sure if one is preferred over another by convention though.