麻煩我用C實現FFT

標準的非遞歸基數-2算法。我意識到我在這裏對於性能並不是最佳的（比如重複的trig調用），我只是想在優化和移植到VHDL之前做好準備。我正在使用complex.h頭文件。麻煩我用C實現FFT

我得到長度爲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

編輯：我意識到我的問題：我的位反轉循環是完全錯誤的。我應該想到這個看似簡單的部分，並且對它進行了更徹底的測試。它爲最後一個案子起作用的原因是因爲它有重複的價值觀，這巧妙地賦予了與工作逆轉相同的向量。這裏是固定循環：

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; 
     } 
    }

來源

2014-09-29 Matthew Hampsey

我不熟悉FFT代碼。你能否提供一個大小爲16或32的樣本數據集，以及你得到的輸出（和打印該答案的代碼），並指出你認爲應該得到的答案。它可能有幫助。 – 2014-09-29 15:17:06

@JonathanLeffler我已經添加了請求的項目 – 2014-09-29 15:27:13

wild guess：int overflow？ – 2014-09-29 15:28:06

這部分2*pi*k/Np讓編譯器決定類型轉換。我也會檢查這些值是否正確計算。例如，我使用2*pi*k/Np; 和2.0*pi*(double) k*(double) Np;時得到不同的結果

（MODS：我寧願將其添加爲評論，但我沒有足夠的聲譽隨意移動我的反應:)）。

來源

2014-09-29 21:21:29 3150

謝謝，我已更改爲明確的演員表。我意識到我的實際問題是：我的位反轉循環不正確。我編輯了正確的代碼。 – 2014-09-30 01:00:46

麻煩我用C實現FFT

回答

相關問題