1

I am currently working on a problem that requires FFT for convolution, however when I brought in my FFT template from my archive I realize that there's something wrong with the output.

eg:

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)

The code:

#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

I tried to check the bit reverse and the root of unity table yet I didn't find any problems in those two parts. I also checked some online materials to verify the steps but there's nothing odd looking to me.

Would someone mind help me to find out what's the problem in this template?

Thanks in advance.

Edit: I decided to rely on another template at the end, thanks for all the reply from everyone.

Community
  • 1
  • 1
haleyk
  • 97
  • 1
  • 8
  • 4
    I doubt, that anyone but you can debug this code. It is very dense and without any comments. – OutOfBound Aug 09 '17 at 08:07
  • 1
    using variable names such as `lll` is fine as long as you are the only person to ever read the code :P – 463035818_is_not_an_ai Aug 09 '17 at 08:12
  • I would suggest to to transform (or inverse transform) unit vectors with only one component non-zero, in that way you might see more easily what coefficients are off – 463035818_is_not_an_ai Aug 09 '17 at 08:16
  • 1
    And `#define op operator` ? Seriously. Not to mention reinventing `std::complex`. That's not even C++11, `std::complex` has been part of every C++ standard. – MSalters Aug 09 '17 at 10:09
  • see [How to compute Discrete Fourier Transform?](https://stackoverflow.com/a/26355569/2521214) you will find mine C++ implementations in the linked answers there for (I)DFT and (I)DFFT both 1D and 2D. I do not see any recursion in your code so I assume it is just DFT not DFFT ... – Spektre Aug 09 '17 at 10:09
  • @MSalters I had to use base instead of std::complex because of its slow runtime performance. The code is intended for solving a problem on an online judge so I didn't follow the naming standards, sorry if that disgusts you. – haleyk Aug 09 '17 at 10:38

1 Answers1

4

It looks like you have the wrong sign for your weights (which means you're probably doing an inverse FFT instead of a forward FFT - remember that for the forward transform the weights are exp(-j * theta)). Try changing:

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

to:

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

This seems to give the correct results for your simple test case. I haven't tried testing it with anything more demanding.

Coincidentally another user recently made exactly the same mistake in a MATLAB implementation. I guess that - sign is easy to miss.

Note also that your code is quite inefficient - you might want to consider using a simple, proven FFT library, like KissFFT instead.

Paul R
  • 208,748
  • 37
  • 389
  • 560
  • Thanks, now the forward pass works fine. I went ahead with @tobi303's suggestion and tried it with convoluting simple unit vectors yet the inverse pass doesn't work now. Would you mind take a further look in it? [code](http://ideone.com/PiezPP) – haleyk Aug 09 '17 at 10:06
  • @haleyk: sorry, I really don't have time to unravel the code - it's rather terse and cryptic. You need to decide whether it's better to spend a few hours debugging it or switch to a tried-and-tested (and much more efficient) third-party library, as suggested above. – Paul R Aug 09 '17 at 10:16