FFT（快速傅里叶变换）——从入门到放弃

FFT（Fast Fourier Transformation）就是“快速傅里叶变换”的意思，它是一种用来计算DFT（离散傅里叶变换）和IDFT（离散傅里叶逆变换）的一种快速算法。这种算法运用了一种高深的数学方式、把原来复杂度为$O(n^2)$的朴素多项式乘法转化为了$O(nlogn)$的算法。 最常用的FFT算法是Cooley和Turkey提出的Radix-2 FFT（基-2FFT），把一个DFT拆成两个一半大小的DFT。这篇文章只介绍Radix-2 FFT。此方法因为实现简单倍受欢迎。Radix-2 FFT适用于N为2的整数次方的FFT。 此外还有Radix-4 FFT（基-4FFT），Radix-8 FFT（咳咳FFT），它们效率更高，但适用范围太窄。 综合了Radix-2的适用范围和Radix-4的效率的算法是SplitRadix FFT（分裂基FFT），名字听上去就很帅，实际玩起来也挺帅的……效率比Radix2快一点吧大概，但实现起来复杂很多。 ## 一、朴素的多项式乘法 设$f(x)$和$g(x)$各表示一个$n-1$次多项式即： $$f(x)=\sum^{n-1}_{i=0}a_i*x^i$$ $$g(x)=\sum^{n-1}_{i=0}b_i*x^i$$ 将所有的两项系数相乘后再相加，利用这种方法计算多项式乘法复杂度为$O(n^2)$。对于1e6的数据范围来说，它的时间复杂度是无法接受的。 ## 二、多项式的点值表示 将$n$个互不相同的$x$带入多项式，会得到$n$个不同的取值$y$。则该多项式被这$n$个点$(x1,y1)$,$(x2,y2)$,…,$(xn,yn)$唯一确定。其中 $$y_i= \sum ^{n-1}_{j=0}a_j*x^j_i$$ 已知$f(x)$与$g(x)$的点值表示，可以在$O(n)$的时内求出$f(x)*g(x)$的点值表示。 但是对于系数表示的多项式，在转换为点值表示时，由于要选$n$个点，每个点$O(n)$，这种表示法计算多项式乘法的时间复杂度仍然为$O(n^2)$，而如果利用高斯消元转换回系数表示法，所需的时间还要更多。 可以优化吗？ ## 三、复数 可以借助复数优化。 ### i的定义 $i=\sqrt{-1}$，这就是$i$的定义。虚数的出现，把实数数系进一步扩张，扩张到了复平面，也使得开方运算终于封闭。实数轴已经被自然数、整数、有理数、无理数塞满了，虚数只好向二维要空间了。 ### 复数的表示 形如$a+bi$的数叫复数。在复平面中，$x$轴代表实数，$y$轴（除原点外的点）代表虚数，从原点$(0,0)$到$(a,b)$的向量表示复数$a+bi$。这个向量的模长就是这个复数的模长。以逆时针为正方向，从$x$轴正半轴到已知向量的转角的有向角叫做幅角。 ``` struct complex { double x,y; complex (double a=0,double b=0){x=a,y=b;} }; ``` ### 复数的运算 #### 加法、减法 与向量相同。 ``` complex operator + (complex a,complex b) { return complex(a.x+b.x , a.y+b.y); } complex operator - (complex a,complex b) { return complex(a.x-b.x , a.y-b.y); } ``` #### 乘法 由定义知，两个复数$a+bi$与$c+di$的乘积： $$(a+bi)*(c+di)$$ $$=ac+adi+bci-bd$$ $$=(ac-bd)+(bc+ad)i$$ 所以： ``` complex operator * (complex a,complex b) { return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x); } ``` ### 单位根 #### 单位根的定义 在复平面上，以原点为圆心，$1$为半径作圆，所得的圆叫单位圆。可以发现，单位圆上的两个复数相乘，结果仍在单位圆上，且幅角等于原来两个复数的幅角之和。 证明：设复数$cos\theta+sin\theta i$与$cos\lambda+sin\lambda i$ $$(cos\theta+sin\theta i)*(cos\lambda+sin\lambda i)$$ $$=cos\theta cos\lambda-sin\theta sin\lambda+cos\theta sin\lambda i+sin\theta cos\lambda i$$ $$=cos(\theta+\lambda)+sin(\theta+\lambda)i$$ #### 单位根的求法 所以，以原点为起点，圆的$n$等分点为终点，做$n$个向量，设幅角为正且最小的向量对应的复数为$\omega _n$，称为$n$次单位根。显然，其余$n-1$个复数为$\omega^2_n$,$\omega^3_n$,…,$\omega^n_n$($\omega^0_n=\omega^n_n=1$)。 对三角函数稍有常识的人就能看出，$\omega^k_n=cos\frac{2\pi*k}{n}+sin\frac{2\pi*k}{n}i$。 #### 单位根的性质 1.$\omega^{2k}_{2n}=\omega^k_n$ 证明： $$\omega^{2k}_{2n}=cos\frac{2\pi*2k}{2n}+sin\frac{2\pi*2k}{2n}i$$ $$=\omega^k_n=cos\frac{2\pi*k}{n}+sin\frac{2\pi*k}{n}i=\omega^k_n$$ 2.$\omega^{k+\frac{n}{2}}_n=-\omega^k_n$ 证明： $$\omega^{k+\frac{n}{2}}_n=cos\frac{2\pi*(2k+n)}{2n}$$ $$=cos(\frac{2\pi*2k}{2n}+\pi)=-\omega^k_n$$ ## 四、DFT（离散傅里叶变换） 我们可以把$n$次单位根代入两个多项式，获得两个多项式的点值表示。 且慢，一共$n$个单位根，每个单位根代入$n$次，这岂不是$O(n^2)$的时间复杂度吗？ 其实是可以优化的： $$f(x)=a_0+a_1*x+a_2*x^2+......+a_{n-2}*x^{n-2}+a_{n-1}*x^{n-1}$$ 按奇偶分组： $$f(x)=(a_0+a_2*x^2+......+a_{n-2}*x^{n-2})+(a_1*x+a_3*x^3+......+a_{n-1}*x^{n-1})$$ 设： $$f_0(x)=a_0+a_2*x+......+a_{n-2}*x^{\frac{n-2}{2}}$$ $$f_1(x)=a_1+a_3*x+......+a_{n-1}*x^{\frac{n-2}{2}}$$ 显然， $$f(x)=f_0(x^2)+xf_1(x^2)$$ 将$\omega^k_n$代入： $$f(\omega^k_n)=f_0(\omega^{2k}_n)+\omega^k_nf_1(\omega^{2k}_n)$$ 再将$\omega^{k+\frac{n}{2}}_n$代入： $$f(\omega^{k+\frac{n}{2}}_n)=f_0(\omega^{2k+n}_n)+\omega^{k+\frac{n}{2}}_nf_1(\omega^{2k+n}_n)$$ $$=f_0(\omega^{2k}_n)-\omega^k_nf_1(\omega^{2k}_n)$$ 两者实部相同，只有虚部相反。因此在算第一个时可以顺便$O(1)$算出第二个。问题就这样缩小了一半。而由于缩小后的问题仍能这样解决，我们就可以递归下去。这样我们就获得了一个$O(nlogn)$的做法。需要注意的是，我们要把多项式的次数补成$2^n-1$次，使得分治可以正常进行。 可以改成非递归版吗？ 观察分组过程： | 原数组 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | | :----------: | :----------: | :----------: | :----------: | :----------: | :----------: | :----------: | :----------: | :----------: | | 分组 | 0 | 2 | 4 | 6 | 1 | 3 | 5 | 7 | | 再分组 | 0 | 4 | 2 | 6 | 1 | 5 | 3 | 7 | 转换成二进制： | 000 | 001 | 010 | 011 | 100 | 101 | 110 | 111 | | :----------: | :----------: | :----------: | :----------: | :----------: | :----------: | :----------: | :----------: | | 000 | 010 | 100 | 110 | 001 | 011 | 101 | 111 | | 000 | 100 | 010 | 110 | 001 | 101 | 011 | 111 | 在这个过程结束后，下标的二进制恰好与原来相反！ 那么我们就可以用下面的方法$O(n)$求出所有反转后的二进制下标了： ``` for(int i=0;i<limit;i++) r[i]= ( r[i>>1]>>1 )| ( (i&1)<<(l-1) ) ; ``` 这也就是所谓的蝴蝶运算。 原理： 我们可以把一个二进制数看成两部分，它的前bit-1位是一部分，它的最后一位是一部分。 一个数的二进制翻转就相当于是把它的最后一位当成首位，然后在后面接上它前bit-1为的二进制翻转。而且在这个循环中我们能保证，在计算“i”的二进制翻转之前1~i-1中的所有数的二进制翻转都已经完成。“i”的前bit-1位的数值其实就是i>>1的值，直接调用i>>1的二进制翻转的结果就相当于调用了“i”的前bit-1位二进制翻转的结果。 至此，我们已经求出了两个多项式的点值表示。然后即可在$O(n)$时间内求出两者乘积的点值表示。 ## 五、IDFT（离散傅里叶逆变换） 可是，我们要的是系数啊！ 所以我们要考虑如何把点值表示法转换为系数表示法，这个过程就是IDFT。利用单位根的一些性质，可以加速这个转换的过程。 先给出结论： $$a_k=\frac{\sum^{n-1}_{i=0}f(\omega^i_n)*(\omega^{-k}_n)^i}{n}$$ 其中$a_k$是$k$次项系数，$f(x)$是多项式在$x$处的点值表示。 我们试着来证明一下：考虑式子右上边的 $$\sum^{n-1}_{i=0}f(\omega^i_n)*(\omega^{-k}_n)^i$$ 将$f(x)$展开：原式$=$ $$\sum^{n-1}_{i=0}(\sum^{n-1}_{j=0}a_j(\omega^i_n)^j)(\omega^{-k}_{n})^i$$ $$=\sum^{n-1}_{i=0}(\sum^{n-1}_{j=0}a_j(\omega^j_n)^i)(\omega^{-k}_{n})^i$$ 显然$(\sum^{n-1}_{j=0}a_j(\omega^j_n)^i)(\omega^{-k}_{n})^i=\sum^{n-1}_{j=0}(a_j(\omega^j_n)^i(\omega^{-k}_{n})^i)$。 所以原式$=$ $$\sum^{n-1}_{i=0}\sum^{n-1}_{j=0}(a_j(\omega^j_n)^i(\omega^{-k}_{n})^i)$$ $$=\sum^{n-1}_{i=0}\sum^{n-1}_{j=0}(a_j(\omega^{j-k}_n)^i)$$ $$=\sum^{n-1}_{j=0}(\sum^{n-1}_{i=0}(\omega^{j-k}_n)^i)a_j$$ 设$S(x)=\sum^{n-1}_{i=0}x_i$ 首项为$1$，公比$\omega^k_n$，由等比数列求和公式得： $$S(\omega^k_n)=\frac{(1-(\omega^k_n)^n)}{1-\omega^k_n}$$ 显然$(\omega^k_n)^n=(\omega_n^1)^{nk}=(\omega^n_n)^k=1$，于是分子$=1-1=0$ 我们惊讶地发现：$S(\omega^k_n)=0$！ 但是别忘了：我们没有考虑$k=0$时的情况。此时的每一项都是$1$，因此和为$n$。 将这个结果回带到$\sum^{n-1}_{j=0}(\sum^{n-1}_{i=0}(\omega^{j-k}_n)^i)a_j$里面： 只有$j=k$时$\sum^{n-1}_{i=0}(\omega^{j-k}_n)^i$值为$n$，否则值为$0$。也就是说原式$=na_k$。 因此： $$a_k=\frac{\sum^{n-1}_{i=0}f(\omega^i_n)*(\omega^{-k}_n)^i}{n}$$ 只要把多项式的点值表示当做系数，代入相应单位根的倒数，用DFT的方法求出结果，所得的结果再除以$n$，就可以得到多项式的系数表示！ 如何求出单位根的倒数呢？只要找到另一个复数，让两个复数模长相同，幅角相加等于$2\pi$即可。画图可知，若单位根$\omega^k_n=cos\theta+sin\theta i$，那么它的倒数就等于$cos\theta-sin\theta i$。 ## 六、总结  ## 七、扩展：傅里叶变换应用初探 MP3作为一种有损压缩格式，通过剔除掉音乐中特定人耳听觉不敏感的频段，减少数据量，获得更大的压缩比。这个“剔除”的过程就用了DFT。 基于Vocaloid歌声合成软件的初音未来，其合成过程用到了大量的DFT。Vocaloid的合成原理基于STFT(短时傅立叶变换)，单音轨合成，每秒至少要进行344次2048点FFT，从而修改信号的频谱特征，Vocaloid合成时大部分时间花在FFT上。所以，你的老婆就是DFT！ 汤姆猫之类的变声器使用DFT改变声音的频谱特征。 均衡器，它通过DFT增强或削弱声音中某个特定频率范围的强度。提高高频和低频部分可以让音乐更明亮有节奏感。 降噪器通过DFT减弱声音中含有噪音的频段，从而实现降噪。我们的手机、电脑的声卡也一般具有降噪器，让你的麦克风输出质量更高。 DFT不仅可以用于音频处理，还可以用于图像处理：给图片加一个水印，直接贴在图片上，很容易被抹掉。一种合理的方法是DFT，在频谱上添加高频水印，再IDFT。 ## 八、代码 ```cpp #include<cstdio> #include<cmath> #include<cctype> const int MAXN=1e7+10; inline int Read() { char c=getchar(); int x=0; while(!isdigit(c)) { c=getchar(); } while(isdigit(c)) { x=(x<<3)+(x<<1)+(c^48); c=getchar(); } return x; } const double Pi=acos(-1.0); struct complex { double x,y; complex (double xx=0,double yy=0){x=xx,y=yy;} }a[MAXN],b[MAXN]; void swap(complex &a,complex &b) { complex c=a; a=b; b=c; } complex operator + (complex a,complex b) { return complex(a.x+b.x , a.y+b.y); } complex operator - (complex a,complex b) { return complex(a.x-b.x , a.y-b.y); } complex operator * (complex a,complex b) { return complex(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x); } int N,M; int l,r[MAXN]; int limit=1; void FFT(complex *A,int opt) { for(int i=0;i<limit;i++) { if(i<r[i]) { swap(A[i],A[r[i]]); } } for(int mid=1;mid<limit;mid<<=1) { complex Wn( cos(Pi/mid) , opt*sin(Pi/mid) ); for(int R=mid<<1,j=0;j<limit;j+=R) { complex w(1,0); for(int k=0;k<mid;k++,w=w*Wn) { complex x=A[j+k],y=w*A[j+mid+k]; A[j+k]=x+y; A[j+mid+k]=x-y; } } } } int main() { int N=Read(),M=Read(); for(int i=0;i<=N;i++) a[i].x=Read(); for(int i=0;i<=M;i++) b[i].x=Read(); while(limit<=N+M) { limit<<=1; ++l; } for(int i=0;i<limit;i++) { r[i]= ( r[i>>1]>>1 )| ( (i&1)<<(l-1) ) ; } FFT(a,1); FFT(b,1); for(int i=0;i<=limit;i++) { a[i]=a[i]*b[i]; } FFT(a,-1); for(int i=0;i<=N+M;i++) { printf("%d ",(int)(a[i].x/limit+0.5)); } return 0; } ```