c - 我在 C 中的 FFT 实现有问题

Question

标准非递归 radix-2 算法。我意识到我这里的性能并不是最佳的(比如重复的触发调用)，我只是想在优化和移植到 VHDL 之前把它做好。我正在使用 complex.h header 。

对于长度为 8 的数组，我得到了正确的结果，但对于长度为 16 或 32 的数组，我得到了正确的结果。但是，对于实际信号，我们仍然得到正确的偶对称结果，只是值是错误的。这让我相信排序是正确的，但也许 trig 调用较小的值有问题？谁能帮忙？

很抱歉我的代码没有得到很好的注释，但老实说它对 FFT 算法没有多大帮助。

void FFT(double complex x[], int N){

    //bit reverse

int s = log2(N);

for(int i = 0; i < N/2; i++){
    int h = bitrev(i, s);

    printf("%u ", h);
    double complex temp = x[i];
    x[i] = x[h];
    x[h] = temp;
}

unsigned int Np = 2; //NUM POINTS IN EACH BLOCK. INITIALLY 2
unsigned int Bp = (N/2);
for(int i = 0; i < s; i++){
    int NpP = Np>>1; //num butterflies
    int BaseT = 0;
    for (int j = 0; j < Bp; j++){

        int BaseB = BaseT + NpP;

        for(int k = 0; k < NpP; k++){

            double complex top = x[BaseT + k];
            double complex bot = (ccos(2*pi*k/Np) - I*csin(2*pi*k/Np))*x[BaseB + k];
            x[BaseT + k] = top+bot;
            x[BaseB + k] = top-bot;
        }
        BaseT = BaseT + Np;
    }
    Bp =  Bp>>1;
    Np = Np<<1;
}

输出打印:

 for(int i = 0; i < LENGTH; i++){
    printf("%f + %fj\n", creal(x[i]), cimag(x[i]));
}

这是一个长度为 8 的输入，以及正确的输出:

double complex x[LENGTH] = {1.0, -1.0, 1.0, -1.0, 10.0, -1.0, 5.0, -1.0};

输出:

13.000000 + 0.000000j
-9.000000 + 4.000000j
5.000000 + 0.000000j
-9.000000 + -4.000000j
21.000000 + 0.000000j
-9.000000 + 4.000000j
5.000000 + 0.000000j
-9.000000 + -4.000000j

这是一个长度为 16 的输入，我得到的错误输出，以及它应该是什么:

double complex x[SIZE] = {10.0, -1.0, 5.0, -1.0, 4.0, 4.0, 2.0, 6.0, 9.0, 3.0, 8.0, 4.0, 5.0, 4.0, 3.0, 8.0};

输出(不正确):

73.000000 + 0.000000j
-4.882497 + 10.451032j
12.535534 + 5.020815j
6.975351 + 7.165394j
12.000000 + 7.000000j
-0.732710 + 0.094326j
5.464466 + 19.020815j
2.639856 + 3.379964j
19.000000 + 0.000000j
2.639856 + -3.379964j
5.464466 + -19.020815j
-0.732710 + -0.094326j
12.000000 + -7.000000j
6.975351 + -7.165394j
12.535534 + -5.020815j
-4.882497 + -10.451032j

正确的输出:

73.000000 + 0.000000j
-4.175390 + 10.743925j
13.535534 + 4.020815j
6.268244 + 5.458287j
10.000000 + 7.000000j
-1.439817 + 1.801433j
6.464466 + 20.020815j
3.346963 + 3.087071j
19.000000 + 0.000000j
3.346963 + -3.087071j
6.464466 + -20.020815j
-1.439817 + -1.801433j
10.000000 + -7.000000j
6.268244 + -5.458287j
13.535534 + -4.020815j
-4.175390 + -10.743925j

对于这个长度为 16 的例子，它给出了正确的输出(!？):

double complex x[SIZE] = {30.0, -1.0, 4.0, -6.0, 4.0, 9.0, 2.0, 6.0, 9.0, 3.0, -8.0, 4.0, -5.0, 4.0, 3.0, 8.0};

66.000000 + 0.000000j
22.604378 + -9.862676j
43.535534 + 28.091883j
24.900438 + 4.851685j
37.000000 + -3.000000j
-1.285214 + 2.408035j
36.464466 + 10.091883j
37.780399 + 23.693673j
12.000000 + 0.000000j
37.780399 + -23.693673j
36.464466 + -10.091883j
-1.285214 + -2.408035j
37.000000 + 3.000000j
24.900438 + -4.851685j
43.535534 + -28.091883j
22.604378 + 9.862676j

编辑:我已经意识到我的问题:我的位反转循环是完全错误的。我应该考虑那个看似简单的部分并更彻底地测试它。它对最后一个案例起作用的原因是因为它有重复的值，巧合的是给出了与工作反转相同的 vector 。这是固定循环:

int s = log2(N);

    for(int i = 0; i < N; i++){
        int h = bitrev(i, s);
        if(i < h){
        double complex temp = x[i];
        x[i] = x[h];
        x[h] = temp;
        }
    }

Answer 1

这部分 2*pi*k/Np 让编译器决定类型转换。 我还会检查是否正确计算了这些值。 例如，我在使用 2*pi*k/Np; 时得到不同的结果 和 2.0*pi*(double) k*(double) Np;

(Mods:我宁愿将此添加为评论，但我没有足够的声誉。请随意移动我的回复:))