代码之家  ›  专栏  ›  技术社区  ›  haleyk

寻求有关FFT模板的帮助

  •  1
  • haleyk  · 技术社区  · 8 年前

    我目前正在研究一个需要FFT进行卷积的问题,然而,当我从存档中引入我的FFT模板时,我意识到输出有问题。

    如:

    输入:(0,0)(0,0)(4166667,0)(1,0)

    模板输出:(4166668,0)(-4166667, -1 1.

    代码:

    #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
    

    我试图检查位反转和单位表的根,但我没有发现这两个部分有任何问题。我还查看了一些在线材料来验证这些步骤,但在我看来没有什么奇怪的地方。

    有人能帮我找出这个模板中的问题吗?

    编辑:我决定在最后依赖另一个模板,谢谢大家的回复。

    1 回复  |  直到 5 年前
        1
  •  4
  •   Paul R    8 年前

    看起来你的权重符号不对(这意味着你可能在做逆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 相反