声明
首先,我需要声明,本文是在转载的基础上稍微修饰的,经过原创作者 ZLH_HHHH(佐理慧学姐) 的许可方才转载并修饰的,由于我就是初学者,并且是数学渣滓,所以我学姐建议我写一下残疾人手册,我当然是欣然接受!!!
正文:
文章写的有点急。有错误的地方望指出
我学习 FFT 是一个比较慢的过程。 期间反反复复。 我写这篇博文只是一个非常浅显的理解。同时也可以帮助初学者在学习FFT的时候。有所偏重。避免太多思维上的负担。
直接正题吧:
首先: DFT 本身并不负责多项式之间的乘法。
DFT 只是一种变换。
FFT 则是DFT的快速算法。(分治 提高效率)
利用 FFT 。我们快速的将多项式变换为利于计算的形式
用这种方便计算的形式计算出来两个多项式的乘积。
这时候我们虽然已经得到目标多项式。但其形式并不是我们想要的
所以 之后利用 FFT 的 逆运算 又快速的变换回去
我们记: FFT 的逆运算为 FFT− 1
对于一个
n
次多项A的表示。最常见的形式(系数表达):
A(x)=∑i=0n−1aixi
这里的n次多项式是指最高次项指数为 n−1 的多项式
2
小学生手算 系数形式 的 多项式乘法 的 复杂度是 O(n2)
如果我们知道了一个
n
次多项式 的曲线上的n 个不同的的点。
我们是可以计算出来这个多项式的
why?
把系数看做未知数。列出来n个方程组。解出来这n个系数。
也就是说,给出曲线上的n个点:
<(x0,y0),(x1,y1),(x2,y2),.....(xn−1,yn−1)>
<script id="MathJax-Element-18" type="math/tex; mode=display"><(x_0,y_0),(x_1,y_1),(x_2,y_2),.....(x_{n-1},y_{n-1})></script>
我们可以确定其系数形式;
我们称这种表达一个多项式的方法叫:点值表达
他有很多优点。比如可以 O(n) 时间内计算一个多项式乘法
再给出一个多项式:
B(x)=∑i=0n−1bixi
设 A(x) 与 B(x) 在取相同x的点值表达为:
<(x0,A0),(x1,A1),(x2,A2),.....(x2n−2,A2n−2)><(x0,B0),(x1,B1),(x2,B2),.....(x2n−2,B2n−2)>
<script id="MathJax-Element-23" type="math/tex; mode=display"><(x_0,A_0),(x_1,A_1),(x_2,A_2),.....(x_{2n-2},A_{2n-2})>\\ <(x_0,B_0),(x_1,B_1),(x_2,B_2),.....(x_{2n-2},B_{2n-2})></script>
这里 Ai=A(xi) , Bi=B(xi)
则 A(x)∗B(x) 为:
<(x0,A0B0),(x1,A1B1),(x2,A2B2),.....(x2n−2,A2n−2B2n−2)>
<script id="MathJax-Element-27" type="math/tex; mode=display"><(x_0,A_0B_0),(x_1,A_1B_1),(x_2,A_2B_2),.....(x_{2n-2},A_{2n-2}B_{2n-2})></script>3
非常方便。顺其自然;
之所以 A(x) 与 B(x) 多取一些点。是因为 A(x)∗B(x) 的次数界增加。
取点不足会导次数界达不到 2n−1
对于一个
n
次多项式。随机求其n个不同的点的朴素方法复杂度是O(n2)
假设 n为偶数。那么我们把 A(x) 。重组为两个多项式:
设其系数的下标为偶数组成的多项式为 A[0](x) :
A[0](x)=a0+a2x+a4x2+...+an−2xn2+1
设其系数的下标为奇数组成的多项式为 A[1](x) :
A[1](x)=a1+a3x+a5x2+...+an−1xn2+1
所以有:
A(xi)=A[0](x2i)+xiA[1](x2i)
4
好像并不会优化运算。
因为
<x20,x21,....x2n−1>
<script id="MathJax-Element-41" type="math/tex; mode=display">
</script>
依然有可能会组成n个不同的数字。
那么我们要计算的规模不会减半。
如果我们恰当选取一些x。是否会优化运算呢。
例如:
<x0,x1,....xn2−1,−x0,−x1,....−xn2−1>
<script id="MathJax-Element-42" type="math/tex; mode=display">
</script>
各个数字平方后。得到的集合大小减半:
<x20,x21,....x2n2−1>
<script id="MathJax-Element-43" type="math/tex; mode=display">
</script>
因为 (−x)2=x2
那么
A(−x)=A[0](x2)−xA[1](x2)A(x)=A[0](x2)+xA[1](x2)
但是这种关系不会传递下去。平方后得到的数字全是不小于0的数字。
再次递归又回到了原来的形式。很尴尬。
不过这启迪我们。取相反数。然后平放。可以把问题规模减半。但再一次递归就失效了。是否存在不失效的取值呢。
<script id="MathJax-Element-46" type="math/tex; mode=display"></script>
如果对于一个有偶数个元素的数字集合。每个元素平方后。得到的新集合。去除重复元素后。集合大小能够减半。并且得到的新集合如果为偶数。新集合依然满足上面的性质。我们称这个集合有 折半 的性质。
如果我们可以快速的找到一个满足折半性质的自变量
x
的取值集合.
分治就是可行的。
有一种东西。可以满足我们的需求。
那就是 —— n次单位复数根:
(使用复数确实让人有顾虑。这也是我学习FFT时最大 思想负担)
我们定义 i=−1‾‾‾√
那么形如
cosθ+i∗sinθ
的复数有诸多良好的性质。
例如:
(cosθ+i∗sinθ)k=cos kθ+i∗sin kθ
上面的性质可以用数学归纳法得到(较为简单。这里不做证明
我们记:
Wn=cos2πn+i∗sin2πn
我会告诉你 多项式在x取下面的值时。有助于我们进行DFT
<W0n,W1n,.....Wn−1n>
<script id="MathJax-Element-54" type="math/tex; mode=display">
</script>
上面的n个复数称为。n次单位复数根。(n次单位复数根是指这n个数)
如果我们x取 Wkn 时。根据刚才的分解:
A(Wkn)=A[0]( (Wkn)2)+WknA[1]((Wkn)2)
而
(Wkn)2=W2kn=cos 2k∗2πn+i∗sin2k∗2πn=cos 2πkn2+i∗sin2πkn2=Wkn2
因为三角函数的周期性:
我们有:
Wkn=cosk∗2πn+i∗sink∗2πn=cos(k mod n)∗2πn+i∗sin(k mod n)∗2πn=Wk mod nn
所以:
(Wkn)2=Wkn2=Wk mod n2n2
则:
<(W0n)2,(W1n)2,...(Wn−1n)2>等效于<W0n2,W1n2,...Wn2−1n2,W0n2,W1n2,...Wn2−1n2>
<script id="MathJax-Element-60" type="math/tex; mode=display"><(W_n^0)^2,(W_n^1)^2,...(W_n^{n-1})^2>\\等效于\\
</script>
惊不惊喜,意不意外。
x
的取值集合取单位复数根。不但满足折半的性质。而且还有一定的规律性。与原集合保持一致。
这意味着我们只需要计算取:
<W0n2,W1n2,...Wn2−1n2>
的 A[0],A[1]
问题范围减半。并且如果n为偶数。再次平方依然集合大小减半。
记:
A[0]k=A[0](Wkn2)=A[0]((Wkn)2)A[1]k=A[1](Wkn2)=A[1]((Wkn)2)
那么,当我们得到:
<A[0]0,A[0]1,..A[0]n2−1,A[1]0,A[1]1,...A[1]n2−1>
<script id="MathJax-Element-65" type="math/tex; mode=display">
</script>
这意味着我们得到了所有的 A[0]((Win)2),A[1]((Win)2)
则:
A( Wkn)=A[0]k modn2+WknA[1]k mod n2
即得到了:
<A0,A1,A2,...An−1>
<script id="MathJax-Element-68" type="math/tex; mode=display">
</script>
其中 Ak=A(Wkn)
当然。FFT的计算 n要取2的整数次幂 。这是因为每次减半。我们可以把不足的系数用0填充。
上面过程可以看作 ,取单位复数根计算目标y向量,如下图:
⎡⎣⎢⎢⎢⎢⎢(W0n)0(W1n)0⋮(Wn−1n)0(W0n)1(W1n)1⋮(Wn−1n)1⋯⋯⋱⋯(W0n)n−1(W1n)n−1⋮(Wn−1n)n−1⎤⎦⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢a0a1⋮an−1⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢A0A1⋮An−1⎤⎦⎥⎥⎥⎥
本着尽可能简单的原则。我们不在特别的说这个矩阵。
因为之前说关于系数的n个方程必然可解。所以上面的矩阵:
⎡⎣⎢⎢⎢⎢⎢(W0n)0(W1n)0⋮(Wn−1n)0(W0n)1(W1n)1⋮(Wn−1n)1⋯⋯⋱⋯(W0n)n−1(W1n)n−1⋮(Wn−1n)n−1⎤⎦⎥⎥⎥⎥⎥
必然存在逆矩阵5。(有线性代数基础的,或许不陌生)
可能你还是有点迷惑。没关系:
如果我们把FFT看作计算上面特别的矩阵相乘的算法呢。
利用FFT我们可以快速得到:
⎡⎣⎢⎢⎢⎢⎢(W0n)0(W1n)0⋮(Wn−1n)0(W0n)1(W1n)1⋮(Wn−1n)1⋯⋯⋱⋯(W0n)n−1(W1n)n−1⋮(Wn−1n)n−1⎤⎦⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢a0a1⋮an−1⎤⎦⎥⎥⎥⎥−(FFT)−>⎡⎣⎢⎢⎢⎢A0A1⋮An−1⎤⎦⎥⎥⎥⎥
同时 。我可以直接告诉你:
D=⎡⎣⎢⎢⎢⎢⎢(W0n)0(W1n)0⋮(Wn−1n)0(W0n)1(W1n)1⋮(Wn−1n)1⋯⋯⋱⋯(W0n)n−1(W1n)n−1⋮(Wn−1n)n−1⎤⎦⎥⎥⎥⎥⎥
的逆矩阵为:
E=⎡⎣⎢⎢⎢⎢⎢⎢⎢1n(W−0n)01n(W−1n)0⋮1n(W−(n−1)n)01n(W−0n)11n(W−1n)1⋮1n(W−(n−1)n)1⋯⋯⋱⋯1n(W−0n)n−11n(W−1n)n−1⋮1n(W−(n−1)n)n−1⎤⎦⎥⎥⎥⎥⎥⎥⎥
因为我直接给你了答案。。。所以怀疑的话。我们可以计算一下看看。
我们记Y=E*D
Y[t][j]=∑k=0n−1E[t][k]∗D[k][j], 其中t,j∈ [0,n−1]
Y[t][j]=∑k=0n−11nW−tkn∗Wkjn=1n∑k=0n−1(Wj−tn)k
上面可以看作等比数列求和。
当
t==j
时:
Y[t][j]=1n∑k=0n−11k=1
当 t≠j 时,令 u=j−t :
Y[t][j]=1n∑k=0n−1(Wun)k=1n∗(Wun)n−1Wun−1 (等比数列公式)
因为:
(Wun)n=Wunn=Wun mod nn=W0n=1
所以,
t≠j
时:
Y[t][j]=0
所以。Y为单位矩阵。E,D互为逆矩阵得证。
所以:
⎡⎣⎢⎢⎢⎢⎢(W0n)0(W−1n)0⋮(W−(n−1)n)0(W0n)1(W−1n)1⋮(W−(n−1)n)1⋯⋯⋱⋯(W0n)n−1(W−1n)n−1⋮(W−(n−1)n)n−1⎤⎦⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢A0A1⋮An−1⎤⎦⎥⎥⎥⎥−(FFT−)−>⎡⎣⎢⎢⎢⎢na0na1⋮nan−1⎤⎦⎥⎥⎥⎥
所以。 FFT− 其实就是 FFT 的 Wkn 变为 W−kn 得到。
并且所得目标向量每个元素除掉n就可以了6
FFT的高效实现:
通常。我们希望FFT的实现尽可能快。所以FFT算法也都使用迭代结构而不是递归结构。
下面介绍一种常用的去递归方法。
首先。上面介绍的FFT是只能处理2的整数次幂的次数界的多项式
所以不特别说明。都有 n=2k (系数不足的用0补足)
对于递归的第一步:
输入数组
<a0,a1,...an−1>
<script id="MathJax-Element-91" type="math/tex; mode=display">
</script>
被分解为,偶数和奇数两个部分数
<a0,a2,a4,...an−2> , <a1,a3,a5,...an−1>
<script id="MathJax-Element-92" type="math/tex; mode=display">
\ \ \ ,\ \ \
</script>
下标二进制形式右起第一个字符为 0 被分配到左集合
下标二进制形式右起第一个字符为 1 被分配到右集合
更一般的。对于第k次递归:
下标右起第k个字符为 0 被分配到它所在问题的左集合
下标右起第k个字符为 1 被分配到它所在问题的右集合
那么对于一个2进制形式的数字 B。将其表示为(注意: n=2k ):
B=∑t=0k−1bt∗2t
根据之前分析。这个数字递归后将会被分配到的位置为:
P(B)=∑t=0k−1bt∗2k−1−t, (因为每次都删去一半位置)
所以。对于一个数字的二进制形式的前k位对称过来。就是递归后的位置。
例如: k=4 , (0011)2−对称后−>(1100)2
我们调用FFT前对数组进行一次上述重拍。
便可 自底向上迭代 实现FFT
总结:
由于递归结构的FFT较好理解。而理解递归结构的FFT后。
不难写出非递归结构
所以在这里我们对递归结构的FFT基本框架做一个总结。
因为使用到了复数。所以不可避免的 。我们以下计算都在复数范围的
方便起见。我们用符号: Co(a,b) 来表示复数 a+bi
那么多项式:
A(x)=∑t=0n−1atxt=∑t=0n−1Co(at,0)xt
Wn=Co(cos2πn,sin2πn)
Wkn=Co(cos2kπn,sin2kπn)
令:
Y[t]=Co(at,0) , t∈[0,n−1]
当n=1时:
FFT(Y[],n)=Y[0]
n>1时:
令: Y[0][t]=Y[2t] , Y[1][t]=Y[2t+1]
计算 FFT(Y[0],n2) 与 FFT(Y[1],n2)
计算完成后,令:
Y[t]=Y[0][t]+WtnY[1][t]Y[t+n2]=Y[0][t]−WtnY[1][t] (其中t∈[0,n2−1])
返回 Y[];
————————————————————————
我对上述思想解释一下
A在 <W0n,W1n,...Wn−1n> <script id="MathJax-Element-111" type="math/tex"> </script>处的值保存在Y中。
并返回给上一层。所以对于当前调用:
A[0](Wtn2)=Y[0][t]A[1](Wtn2)=Y[1][t]
A(Wtn)=A[0]((Wtn)2)+WtnA[1]((Wtn)2)=A[0](Wtn2)+WtnA[1](Wtn2)
当 t∈[0,n2−1] 时:
A(Wtn)=Y[t]=A[0](Wtn2)+WtnA[1](Wtn2)=Y[0][t]+WtnY[1][t]
A(Wt+n2n)=Y[t+n2]=A[0](Wt+n2n2)+Wt+n2nA[1](Wt+n2n2)=A[0](Wtn2)−WtnA[1](Wtn2)=Y[0][t]−WtnY[1][t]
模版
给出一个迭代结构的FFT模版:
定义复数:
struct Complex
{
double x,y;
Complex(double x1=0.0 ,double y1=0.0)
{
x=x1;
y=y1;
}
Complex operator -(const Complex &b)const
{
return Complex(x-b.x,y-b.y);
}
Complex operator +(const Complex &b)const
{
return Complex(x+b.x,y+b.y);
}
Complex operator *(const Complex &b)const
{
return Complex(x*b.x-y*b.y,x*b.y+y*b.x);
}
};
去递归:
void change(Complex y[],int len)
{
int i,j,k;
for(i=1,j=len/2;i<len-1;i++)
{
if(i<j)swap(y[i],y[j]);
k=len/2;
while(j>=k)
{
j-=k;
k/=2;
}
j+=k;
}
}
迭代结构的FFT
on== 1 时: FFT
on==-1 时: FFT−
void FFT(Complex y[],int len,int on)
{
change(y,len);
for(int h=2;h<=len;h<<=1)
{
Complex wn(cos(on*2*pi/h),sin(on*2*pi/h));
for(int j=0;j<=len;j+=h)
{
Complex w(1,0);
for(int k=j;k<j+h/2;k++)
{
Complex u=y[k];
Complex t=w*y[k+h/2];
y[k]=u+t;
y[k+h/2]=u-t;
w=w*wn;
}
}
}
if(on==-1)
for(int i=0;i<len;i++)
y[i].x/=len;
}