快速傅里叶变换C++

Question

我正在尝试用 C++ 重现此程序流程图 来计算数组的 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 <- real
0 0 0 0 4.24264 -2.71191 -8.61313 7.08239 0 0 0 0 0 0 0 0 <- imaginary

当我应该得到这个时：

(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 编写的代码，因此我尝试严格地将其转换为 cpp，同时也检查我给出的公式。我还实现了一个程序，通过简单的矩阵乘法计算 DFT，以比较结果。

Answer 1

我显然忽略了一个事实，即位反转应该发生在

gamma

位上，这就是整个问题。 itbtr 函数的正确版本是：

int itbtr(int k, int gamma){
    int m = k;
    int newn = 0;
    for(int i = 0; i < gamma; i++){
        newn *= 2;
        newn += m % 2;
        m /= 2;
    }
    return newn;
}