C++实现FFT算法结果异常求助：按流程图编写但输出不符

我正尝试按照流程图用C++实现数组的FFT计算，本以为按流程图实现会很简单，但输出结果始终与预期不符。已重新检查所有复数乘法步骤，确认符合算法要求，但计算结果仍不正确。

我知道可以使用complex等模块通过更现代、简便的方法实现FFT，但我想弄清楚当前代码出错的原因。

测试用例 #

测试数组：{1, 2, 3, 4, 4, 3, 2, 1, 1, 2, 3, 4, 4, 3, 2, 1}

输出对比 #

实际输出 #

实数部分：40 0 0 0 -1.75736 -6.54712 -8.61313 -7.08239 0 0 0 0 0 0 0 0
虚数部分：0 0 0 0 4.24264 -2.71191 -8.61313 7.08239 0 0 0 0 0 0 0 0

预期结果 #

(40,0) (0,0) (-11.6569,-4.82843) (0,0) (0,0) (0,0) (-0.343146,-0.828427) (0,0) (0,0) (0,0) (-0.343146,0.828427) (0,0) (0,0) (0,0) (-11.6569,4.82843) (0,0)

实现代码 #

#include<iostream>
using namespace std;


double xreal[100] = {1, 2, 3, 4, 4, 3, 2, 1, 1, 2, 3, 4, 4, 3, 2, 1}, ximag[100];
int gamma = 4;

int itbtr(int m);

void fft(int gamma, double xreal[100], double ximag[100]){
    int n2, nu1;
    int N = 1<<gamma;
 
    int l = 1;
    n2 = N/2;
    nu1 = gamma - 1;
    int k = 0;


    while(l <= gamma) {
        while(k < N) {
            int I = 1;
            while(I <= n2)
            {
                int M = int(k >> nu1);
                int P = itbtr(M);
  
                double arg = 6.283185 * P / N;

                double cosarg = cos(arg);
                double sinarg = sin(arg);

                double treal = cosarg * xreal[k + n2] - ximag[k + n2] * sinarg;
                double timag = sinarg * xreal[k + n2] + cosarg * ximag[k + n2];

                xreal[k + n2] = xreal[k] - treal;
                ximag[k + n2] = ximag[k] - timag;

                xreal[k] = xreal[k] + treal;
                ximag[k] = ximag[k] + timag;

                k++;
                I++;
            }
            k += n2;
        }
        l++;
        n2 /= 2;
        nu1 --;
        k = 0;
    }

    if(l > gamma){
        while(k <= N - 1){
            int i = itbtr(k);
            if(i > k) swap(xreal[k], xreal[i]), swap(ximag[k], ximag[i]);
            k ++ ;
        }
    }
}

//bit reversal
int itbtr(int k){
    int m = k;
    int newn = 0;
    while(m){
        newn *= 2;
        newn += m % 2;
        m /= 2;
    }
    return newn;
}


int main(){
    fft(gamma, xreal, ximag);
    for(int i = 0; i < pow(2, gamma); i++)
    cout<<xreal[i]<<" ";
    cout<<endl;
    for(int i = 0; i < pow(2, gamma); i++)
    cout<<ximag[i]<<" ";
}

我严格将流程图附带的Pascal代码转换为C++，同时对照了相关公式，还实现了基于矩阵乘法的DFT程序来对比结果，恳请帮忙排查问题。

内容的提问来源于stack exchange，提问作者anonymous_wolf