寻求有关 FFT 模板的帮助

Question

我目前正在处理一个需要 FFT 进行卷积的问题，但是当我从我的档案中引入我的 FFT 模板时，我意识到输出有问题。

例如：

Input: (0, 0) (0, 0) (4166667, 0) (1, 0)

Correct output: (4166668, 0) (-4166667, 1) (4166666, 0) (-4166667, -1)

Template output: (4166668, 0) (-4166667, -1) (4166666, 0) (-4166667, 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

我尝试检查位反转和统一根table，但我没有发现这两部分有任何问题。我也查看了一些在线资料来验证这些步骤，但我看起来没有什么奇怪的。

谁能帮我看看这个模板有什么问题？

提前致谢。

编辑：最后我决定依赖另一个模板，感谢大家的回复。

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 用于您的简单测试用例。我还没有尝试用更苛刻的东西来测试它。

最近巧合的是另一个用户。我猜 - 标志很容易被错过。

另请注意，您的代码效率很低 - 您可能需要考虑使用一个简单的、经过验证的 FFT 库，例如 KissFFT。

寻求有关 FFT 模板的帮助

Looking for help about a FFT template

c++

math

signal-processing

fft