首页 > 代码库 > FFT —— 快速傅里叶变换

FFT —— 快速傅里叶变换

问题:

  已知A[], B[], 求C[],使:

    技术分享

  定义C是A,B的卷积,例如多项式乘法等。

 

  朴素做法是按照定义枚举i和j,但这样时间复杂度是O(n2).

  能不能使时间复杂度降下来呢?

 

点值表示法:

  我们把A,B,C看作表达式。

  即:

    A(x)=a+ a1* x + a2 * x2 +...

  将A={(x1,A(x1)), (x2,A(x2)), (x3,A(x3))...}叫做A的点值表示法。

  那么使用点值表示法做多项式乘法就很简单了:对应项相乘。

  那么,如何将A和B转换成点值表示法,再将C转化回系数表示法(即最初的表示方法)呢?

  如果任取n个点,按照定义计算,那么还是O(n2)的。

  这样就要用到快速傅里叶变换。

 

快速傅里叶变换:

  既然任取n个点,按照定义计算太慢,就要找一些特殊点。

  我们用n个n次单位复数根(1的n次方根,涉及到复数,1的方根不止1和-1)来计算:

    1的n次方根是技术分享,其中i是虚数单位。

  我们定义wn = e^(2∏i)是主n次单位根,那么所有n次单位复数根都是它的次方。

  我们要求出A(wnk),就要采用分治思想。

  我们将奇偶系数分离(先假设n为偶数),即定义

      A1(x)=a+ a2* x + a4 * x2 +...

      A2(x)=a+ a3* x + a5 * x2 +...

  那么A(x)=A1(x2) + xA2(x2)。

  要计算A(wnk)=A1((wnk)2) + wnkA2((wnk)2),

  就要用到(wnk)2 = wn/2k mod (n / 2)(证略)。

  所以A(wnk)=A1(wn/2k mod (n / 2)) + wnkA2(wn/2k mod (n / 2))

  我们发现A1,A2都是n/2项的,且只需要算wn/2k的值,那么这就和开始的问题一样了,可以分治。

  边界也很容易:n=1的时候A1本身就是值。

  合并解。

  A(wnk)=A1(wn/2k mod (n / 2)) + wnkA2(wn/2k mod (n / 2))

  那么可以A(wnk)和A(wnk+n/2)一起算(0<=k<n/2):

    设u = A1(wn/2k), t = wnkA2(wn/2k),

    那么A(wnk) = u + t

      A(wnk+n/2) = A1(wn/2k) + wnk + n/2 A2(wn/2k)

            = A1(wn/2k) + wnk  wnn/2 A2(wn/2k)

            = A1(wn/2k) - wnk A2(wn/2k)

            = u - t

  所以这样就能算出A的点值表示法。

  一个问题:分治要求n是2的幂,不是怎么办? 补0, 直到是2的幂。

 

  剩下的问题:如何把C转化回系数表示法。

快速傅里叶逆变换:

  我们把C做一遍快速傅立叶变换,只是求的是wnn, wnn-1, ..., wn1的值而不是wn0, wn1, ..., wnn-1的值,最后每一项除以n即可。

  证略。

 1 void Rader(complex y[],int len)
 2 {
 3     int i,j,k;
 4     for(i = 1, j = len/2;i < len-1;i++)
 5     {
 6         if(i < j)swap(y[i],y[j]);
 7         k = len/2;
 8         while( j >= k)
 9         {
10             j -= k;
11             k /= 2;
12         }
13         if(j < k)j += k;
14     }
15 }
16 void FFT(complex y[],int len,int on)  //on = 1 快速傅里叶变换, on = 0 快速傅里叶逆变换
17 {
18     Rader(y,len);
19     for(int h = 2;h <= len;h <<= 1)
20     {
21         complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));  //e^ki = cosk + isink
22         for(int j = 0;j < len;j += h)
23         {
24             complex w(1,0);
25             for(int k = j;k < j+h/2;k++)
26             {
27                 complex u = y[k];
28                 complex t = w*y[k+h/2];
29                 y[k] = u+t;
30                 y[k+h/2] = u-t;
31                 w = w*wn;
32             }
33         }
34     }
35     if(on == -1)
36         for(int i = 0;i < len;i++)
37             y[i].r /= len;
38 }
39 //复数实现略

 

FFT —— 快速傅里叶变换