快速傅里叶变换简介
F a s t ? F o u r i e r ? T r a n s f o r m Fast \ Fourier \ Transform Fast?Fourier?Transform ，简称 F F T FFT FFT。这东西似乎在数字通信领域有大用处：（EDA表示很感动）
计算量小的显著的优点，使得FFT在信号处理技术领域获得了广泛应用，结合高速硬件就能实现对信号的实时处理。例如，对语音信号的分析和合成，对通信系统中实现全数字化的时分制与频分制(TDM/FDM)的复用转换，在频域对信号滤波以及相关分析，通过对雷达、声纳、振动信号的频谱分析以提高对目标的搜索和跟踪的分辨率等等，都要用到FFT。可以说FFT的出现，对数字信号处理学科的发展起了重要的作用。——百度百科（快速傅里叶变换词条）
但是本博客只介绍其在算法竞赛中的应用——快速计算两个多项式的乘积。
假设现在有两个多项式： A ( x ) = x 2 + 3 x + 2 A(x)=x^2+3x+2 A(x)=x2+3x+2， B ( x ) = 2 x 2 + 1 B(x)=2x^2+1 B(x)=2x2+1，要求 C ( x ) = A ( x ) × B ( x ) C(x)=A(x) \times B(x) C(x)=A(x)×B(x)。朴素做法是按照乘法分配律，将其中一个多项式的每一项分别与另一个多项式的每一项相乘，如果这两个多项式最高次数为 N N N，这个算法的时间复杂度将会是 O ( n 2 ) O(n^2) O(n2)。然而我们可以用 F F T FFT FFT在 O ( n l o g n ) O(nlogn) O(nlogn)的时间内计算出 C ( x ) C(x) C(x)。 F F T FFT FFT算法的核心，是四个十分精妙的想法，发明这个算法的人应该是天选者。
下面将介绍这四个巧妙的想法，以及 F F T FFT FFT的代码和实现细节。
1. 多项式的另一种表示
一般怎么用计算机存储多项式？
开一个数组 A [ i ] A[i] A[i]，表示多项式 x i x^i xi 项的系数，例如 F ( x ) = 4 x 4 + 6 x 3 + 5 x 2 + 3 x + 2 F(x)=4x^4+6x^3+5x^2+3x+2 F(x)=4x4+6x3+5x2+3x+2，那么 A [ ] = { 4 ， 6 ， 5 ， 3 ， 2 } A[]=\{4，6，5，3，2\} A[]={4，6，5，3，2}。这完全没有问题，形如 F ( x ) = ∑ a i x i F(x)=\sum a_ix^i F(x)=∑ai?xi 的多项式表示方法，叫做“系数表示法”。
然而现在为了更快计算多项式乘积，我们使用多项式的另一种表示法，要理解这种表示方法，可以从最简单的多项式开始。设 A ( x ) = a 0 + a 1 x A(x)=a_0+a_1x A(x)=a0?+a1?x，一次函数就是一种简单的多项式，它的最高幂次为1。
我们知道 A ( x ) = a 0 + a 1 x A(x)=a_0+a_1x A(x)=a0?+a1?x 不但是一个多项式，还是一条直线。由于两点能确定一条直线，我们在平面中随便取两点 ( x 1 ， y 1 ) (x_1，y_1) (x1?，y1?)和 ( x 2 ， y 2 ) (x_2，y_2) (x2?，y2?)，那么系数 a 0 a_0 a0?， a 1 a_1 a1? 也就随之确定了。而我们现在在讲的“多项式的另一种表示方法”，就是用平面上的点来确定多项式。
高次多项式的图像是一条曲线。显然，无论你在平面上取多少点，也无法唯一确定一条曲线，但却可以确定一个多项式。因为多项式“并不是普通的（或者说任意的）曲线”。这里有一个结论：平面上的
N
+
1
N+1
N+1 个点，能唯一确定一个
N
N
N 次多项式。例如：你在平面上随机取四个点，那么有且只有一个三次函数同时穿过这四个点。这个结论的正确性需要用线性代数的知识来证明，详细写可能需要一篇新的博客，所以我只能简单说明一下正确性。其实是不会
设
F
(
x
)
=
∑
i
=
0
i
=
n
a
i
x
i
F(x)=\sum_{i=0}^{i=n}a_ix^i
F(x)=∑i=0i=n?ai?xi，我们选取的点为
{
x
0
，
x
1
，
x
2
，
.
.
.
，
x
n
}
\{x_0，x_1，x_2，...，x_n\}
{x0?，x1?，x2?，...，xn?}。那么：
[
F
(
x
0
)
F
(
x
1
)
?
F
(
x
n
)
]
=
[
1
x
0
x
0
2
?
x
0
n
1
x
1
x
1
2
?
x
1
n
?
?
?
?
?
1
x
n
x
n
2
?
x
n
n
]
[
a
0
a
1
?
a
n
]
\left[ \begin{matrix} F(x_0) \\ F(x_1) \\ \vdots \\F(x_n) \end{matrix} \right]= \left[ \begin{matrix} 1 & x_0 &x^2_0 & \cdots & x_0^n \\ 1 & x_1 &x^2_1 & \cdots & x_1^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n &x^2_n & \cdots & x_n^n \\ \end{matrix} \right] \left[ \begin{matrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{matrix} \right]
??????F(x0?)F(x1?)?F(xn?)???????=??????11?1?x0?x1??xn??x02?x12??xn2???????x0n?x1n??xnn??????????????a0?a1??an????????
易证中间矩阵的行列式不为0，所以方程组只有一组解。好了证毕。
这种表示多项式的方法叫做“点值表示法”。假设对于两个多项式
f
(
x
)
f(x)
f(x)和
g
(
x
)
g(x)
g(x)，我们选取相同的
x
x
x 序列
{
x
0
,
x
1
.
.
.
.
,
x
n
}
\{ x_0, x_1....,x_n\}
{x0?,x1?....,xn?}，那么两个多项式可以表示为：
f
(
x
)
=
(
?
(
x
0
,
f
(
x
0
)
)
,
(
x
1
,
f
(
x
1
)
)
,
.
.
.
,
(
x
n
,
f
(
x
n
)
)
?
)
g
(
x
)
=
(
?
(
x
0
,
g
(
x
0
)
)
,
(
x
1
,
g
(
x
1
)
)
,
.
.
.
,
(
x
n
,
g
(
x
n
)
)
?
)
f(x)=(\ (x_0,f(x_0)),(x_1,f(x_1)),...,(x_n,f(x_n))\ ) \\ g(x)=(\ (x_0,g(x_0)),(x_1,g(x_1)),...,(x_n,g(x_n))\ )
f(x)=(?(x0?,f(x0?)),(x1?,f(x1?)),...,(xn?,f(xn?))?)g(x)=(?(x0?,g(x0?)),(x1?,g(x1?)),...,(xn?,g(xn?))?)
设
F
(
x
)
=
f
(
x
)
×
g
(
x
)
F(x)=f(x) \times g(x)
F(x)=f(x)×g(x)，则
F
(
x
)
=
{
(
x
0
,
f
(
x
0
)
g
(
x
0
)
)
,
(
x
1
,
f
(
x
1
)
g
(
x
1
)
)
,
.
.
.
,
(
x
n
,
f
(
x
n
)
g
(
x
n
)
)
}
F(x)=\{(x_0,f(x_0)g(x_0)),(x_1,f(x_1)g(x_1)),...,(x_n,f(x_n)g(x_n))\}
F(x)={(x0?,f(x0?)g(x0?)),(x1?,f(x1?)g(x1?)),...,(xn?,f(xn?)g(xn?))}
容易发现，如果已知两个多项式的点值表示，求两个多项式乘积多项式的点值表示，可以在线性时间内处理完。如果能找到一种方法，能快速将一个多项式在系数表示和点值表示之间转化，就能实现快速多项式乘法。
2.系数到点值的转换（DFT）
既然任意 N + 1 N+1 N+1 个点就可以确定一个 N N N 阶多项式，那么如果我直接暴力随机选 N + 1 N+1 N+1 个横坐标，带入多项式求值会怎样呢？容易想到这样做的时间复杂度是 O ( n 2 ) O(n^2) O(n2)（找 N + 1 N+1 N+1个点 O ( N ) O(N) O(N)，每个点求值 O ( N ) O(N) O(N)），这样就又回到了原点，没有达到加速的目的，所以需要找另外的方法。
依旧考虑最简单的情况：现在，我想在二次函数 f ( x ) = x 2 f(x)=x^2 f(x)=x2上，任意取 8 个点，并求出这8个点的函数值，有没有什么比较快的方法？容易想到，其实我只需要任意取四个 x x x 就行了。因为如果我知道了 x = x 0 x=x_0 x=x0? 的函数值，我立刻就知道了 x = ? x 0 x=-x_0 x=?x0? 的函数值，因为 f ( x ) = x 2 f(x)=x^2 f(x)=x2 是一个偶函数。奇函数同理，只不过我需要在函数值上再加一个负号。总之：我们可以利用奇偶性来减少选取点的数目。
下面的问题就是，如何将一个一般的多项式，转换成可以利用奇偶性进行优化的多项式。我们可以将原多项式，按照每一项
x
x
x 的幂次分组，偶数次幂一组，奇数次幂一组，然后再把奇数次幂组中的每一项提出一个
x
x
x 来，就像下面这样：
F
(
x
)
=
3
x
5
+
2
x
4
+
x
3
+
7
x
2
+
5
x
+
1
?
F
(
x
)
=
(
2
x
4
+
7
x
2
+
1
)
+
x
(
3
x
4
+
x
2
+
5
)
F(x)=3x^5+2x^4+x^3+7x^2+5x+1 \Rightarrow \\ F(x)=(2x^4+7x^2+1)+x(3x^4+x^2+5)
F(x)=3x5+2x4+x3+7x2+5x+1?F(x)=(2x4+7x2+1)+x(3x4+x2+5)
我们可以把左右两部分视为新的，关于
x
2
x^2
x2 的多项式。
F
1
(
x
2
)
=
2
x
4
+
7
x
2
+
1
F
2
(
x
2
)
=
3
x
4
+
x
2
+
5
F_1(x^2) = 2x^4+7x^2+1 \\ F_2(x^2) = 3x^4+x^2+5
F1?(x2)=2x4+7x2+1F2?(x2)=3x4+x2+5
我们只需要解决这两个小多项式的求值，再利用
F
(
x
)
=
F
1
(
x
2
)
+
x
F
2
(
x
2
)
F(x)=F_1(x^2)+xF_2(x^2)
F(x)=F1?(x2)+xF2?(x2) 和
F
(
?
x
)
=
F
1
(
x
2
)
?
x
F
2
(
x
2
)
F(-x)=F_1(x^2)-xF_2(x^2)
F(?x)=F1?(x2)?xF2?(x2) 这两个公式代回即可。最终取的点的形式是：
{
F
(
x
1
)
,
F
(
?
x
1
)
,
F
(
x
2
)
,
.
.
.
.
,
F
(
x
n
/
2
)
,
F
(
?
x
n
/
2
)
}
\{F(x_1), F(-x_1), F(x_2),....,F(x_{n/2}),F(-x_{n/2})\}
{F(x1?),F(?x1?),F(x2?),....,F(xn/2?),F(?xn/2?)}。而两个小多项式都只有原多项式一半的规模，所以这个方法的复杂度是
O
(
n
log
?
n
)
O(n\log{n})
O(nlogn)。但是这个方法还存在一个问题：我们想利用奇偶性和每次选取
x
=
x
0
x=x_0
x=x0? 和
x
=
?
x
0
x=-x_0
x=?x0? 两个点，来减少选取点的次数，也就是说我们每次取的是成相反数的两个点。然而当问题变为关于
F
1
(
x
2
)
F_1(x^2)
F1?(x2) 和
f
2
(
x
2
)
f_2(x^2)
f2?(x2) 时，可以发现：
x
2
x^2
x2 永远为正，我们不再能取一对相反数了。而解决这个问题的方法，被人们视为
F
F
T
FFT
FFT的精髓。
3.单位复根与快速傅里叶变换
可以看出，问题出在多项式的定义域上。一旦问题的规模缩小，原多项式的定义域会失效。那么如果我们把多项式的定义域扩大到复数域呢？这次不从最简单的情况入手了，这次需要一个逆向思维。下面结合图来说明。
假设递归已经进行到了终点，也就是现在多项式只剩下一次项了。也就是说我们需要选取一个点来求值，那么这个点如何选取？显然我可以选
x
=
1
x=1
x=1。然后考虑上一层递归，我要选两个值，我希望他们俩是相反数，还希望他们俩的平方等于 1。显然我可以取
x
=
1
x=1
x=1 和
x
=
?
1
x=-1
x=?1。再往上一层，左边是 1，和第一层一样，关键是右边的 -1。我同样希望上一层的两个数，他们互为相反数，而他们的平方等于 -1。这两个数显然是
i
i
i 和
?
i
-i
?i。
到目前为止，我们得到了四个数：1，-1，i，-i 。如果我的原多项式是一个三阶多项式，我只需要取
x
x
x 等于这四个数字就可以。但如果我的多项式阶数更高，这四个数还不够。我需要至少
N
N
N 个才行。这
N
N
N 个数需要满足
x
N
=
1
x^N = 1
xN=1，而之前学过的代数知识告诉我们，
x
N
=
1
x^N=1
xN=1 在复数域内恰好有
N
N
N 个解。也就是说我们前面得到的四个数字，其实是
x
4
=
1
x^4=1
x4=1 的四个解。
x
n
=
1
x^n=1
xn=1 在复数域上的
n
n
n 个解，均匀地分布在复平面的单位圆上。例如：
x
8
=
1
x^8=1
x8=1 的八个解，就是单位圆的八个八等分点。读者可以自行验证。这样的点叫做单位复根。将(1，0)点编号为0，其余点依次编号。我们定义
ω
n
k
\omega_{n}^k
ωnk? 是
x
n
=
1
x^n=1
xn=1 的第
k
k
k 号单位复根。
（本来想做个单位圆的图，然而没有比较好的几何工具，就咕咕了）
简单推导可得： ω n k = cos ? ( 2 π k n ) + i × sin ? ( 2 π k n ) \omega_{n}^k=\cos(\frac{2\pi k}{n})+i \times \sin(\frac{2\pi k}{n}) ωnk?=cos(n2πk?)+i×sin(n2πk?)
单位复根有三个重要性质，通过简单推导就可以得到：
- ω n n = 1 \omega_{n}^n=1 ωnn?=1
- ω n k = ω 2 n 2 k \omega_{n}^k=\omega_{2n}^{2k} ωnk?=ω2n2k?
- ω 2 n k + n = ? ω 2 n k \omega_{2n}^{k+n}=-\omega_{2n}^k ω2nk+n?=?ω2nk?
那么原来的公式
F
(
x
)
=
F
1
(
x
2
)
+
x
F
2
(
x
2
)
F(x)=F_1(x^2)+xF_2(x^2)
F(x)=F1?(x2)+xF2?(x2)，利用单位复根的性质可得：
同理可得：
因此我们求出了
D
F
T
(
G
(
ω
n
/
2
k
)
)
DFT(G(\omega^k_{n/2}))
DFT(G(ωn/2k?)) 和
D
F
T
(
H
(
ω
n
/
2
k
)
)
DFT(H(\omega^k_{n/2}))
DFT(H(ωn/2k?)) 后，就可以同时求出
D
F
T
(
f
(
ω
n
k
)
)
DFT(f(\omega^k_n))
DFT(f(ωnk?)) 和
D
F
T
(
f
(
ω
n
k
+
n
/
2
)
)
DFT(f(\omega^{k+n/2}_n))
DFT(f(ωnk+n/2?)) 。于是对
G
G
G 和
H
H
H 分别递归 DFT 即可。考虑到分治 DFT 能处理的多项式长度只能是二的整数次幂，否则在分治的时候左右不一样长，右边就取不到系数了。所以要在第一次 DFT 之前就把序列向上补成长度为
2
m
2^m
2m（高次系数补
0
0
0）、最高项次数为
2
m
?
1
2^m-1
2m?1 的多项式。
所以到现在我们可以发现，最初的原多项式代入的 n n n 个点，实际就是 { ω n 0 , ω n 1 , ω n 2 , . . . , ω n n ? 1 } \{\omega_n^0,\omega_n^1,\omega_n^2,...,\omega_n^{n-1}\} {ωn0?,ωn1?,ωn2?,...,ωnn?1?}，其中 n n n 是二的整数次幂。
得到两个多项式乘积后，还需要把这个结果多项式从点值表示变回系数表示。
4.点值到系数的转换（IDFT）
这里我真写不下去了，因为 DFT 和 IDFT 怎么结合到一起的我真啥也没看懂，就知道一个结论：做 DFT 和 IDFT 的区别就在于一个正负号。
1
ω
k
=
ω
k
?
1
=
e
?
2
π
i
k
=
cos
?
(
2
π
k
)
+
i
×
sin
?
(
?
2
π
k
)
\frac{1}{\omega_k}=\omega_k^{-1}=e^{-\frac{2\pi i}{k}}=\cos(\frac{2\pi}{k})+i \times \sin(-\frac{2\pi}{k})
ωk?1?=ωk?1?=e?k2πi?=cos(k2π?)+i×sin(?k2π?)
因此实际在做
F
F
T
FFT
FFT时，通常会传入一个参数，来决定是做
D
F
T
DFT
DFT 还是
I
D
F
T
IDFT
IDFT。
之后要是弄懂了可能会回来写。
相关文章
精彩导读
热门资讯