尋找關於FFT模板的幫助

Question

我目前正在研究需要FFT進行卷積的問題，但是當我從我的存檔中引入FFT模板時，我發現輸出存在問題。

例如：

輸入：（0,0）（0,0）（4166667，0）（1，0）

正確的輸出：（4166668，0）（-4166667， 1）（4166666，0）（-4166667，-1）

模板輸出：（4166668，0）（-4166667，-1）（4166666，0）（-4166667，）

代碼：

#define MAXN 
#define ld long double 
#define op operator 

struct base { 
    typedef ld T; T re, im; 
    base() :re(0), im(0) {} 
    base(T re) :re(re), im(0) {} 
    base(T re, T im) :re(re), im(im) {} 
    base op + (const base& o) const { return base(re + o.re, im + o.im); } 
    base op - (const base& o) const { return base(re - o.re, im - o.im); } 
    base op * (const base& o) const { return base(re * o.re - im * o.im, re * o.im + im * o.re); } 
    base op * (ld k) const { return base(re * k, im * k); } 
    base conj() const { return base(re, -im); } 
}; 

base w[MAXN]; //omega lookup table 
int rev[MAXN]; //reverse lookup table 

void build_rev(int k) { 
    static int rk = -1; 
    if(k == rk)return ; rk = k; 
    for(int i = 1; i < (1<<k); i++) { 
     int j = rev[i-1], t = k-1; 
     while(t >= 0 && ((j>>t)&1)) { j ^= 1 << t; --t; } 
     if(t >= 0) { j ^= 1 << t; --t; } 
     rev[i] = j; 
    } 
} 

void fft(base *a, int k) { 
    build_rev(k); int n = 1 << k; 
    for(int i = 0; i < n; i++) if(rev[i] > i) swap(a[i], a[rev[i]]); 
    for(int l = 2, lll = 1; l <= n; l += l, lll += lll) { 
     if(w[lll].re == 0 && w[lll].im == 0) { 
      ld angle = PI/lll; 
      base ww(cosl(angle), sinl(angle)); 
      if(lll > 1) for(int j = 0; j < lll; ++j) { 
       if(j & 1) w[lll + j] = w[(lll+j)/2] * ww; 
       else w[lll + j] = w[(lll+j)/2]; 
      } else w[lll] = base(1, 0); 
     } 
     for(int i = 0; i < n; i += l) 
     for(int j = 0; j < lll; j++){ 
      base v = a[i + j], u = a[i + j + lll] * w[lll + j]; 
      a[i + j] = v + u; a[i + j + lll] = v - u; 
      } 
    } 
} 

//ideone compiled example: http://ideone.com/8PTjW5 

我試圖檢查位反向統一表的根但我沒有發現在這兩個部分的任何問題。我也檢查了一些在線資料來驗證步驟，但沒有什麼奇怪的東西給我看。

會有人介意幫我找出這個模板有什麼問題嗎？

在此先感謝。

編輯：我決定在最後依靠另一個模板，感謝所有人的回覆。

Answer 1

它看起來像你的權重有錯誤的符號（這意味着你可能做了反FFT而不是正向FFT - 記住正向變換the weights are exp(-j * theta)）。嘗試改變：

base ww(cosl(angle), sinl(angle)); 

到：

base ww(cosl(angle), -sinl(angle)); 

這seems to give the correct results爲您簡單的測試案例。我還沒有嘗試過任何更苛刻的測試。

巧合的是另一位用戶最近made exactly the same mistake in a MATLAB implementation。我猜-標誌很容易錯過。

還要注意，你的代碼效率很低 - 你可能要考慮使用一個簡單的，經過驗證的FFT庫，如KissFFT。