文章写的有点急。有错误的地方望指出
我学习 FFT 是一个比较慢的过程。 期间反反复复。 我写这篇博文只是一个非常浅显的理解。同时也可以帮助初学者在学习FFT的时候。有所偏重。避免太多思维上的负担。
直接正题吧:
首先:DFTDFTDFT本身并不负责多项式之间的乘法。
DFTDFTDFT只是一种变换。
FFTFFTFFT则是DFT的快速算法。(分治提高效率)
利用FFTFFTFFT 。我们快速的将多项式变换为利于计算的形式
用这种方便计算的形式计算出来两个多项式的乘积。
这时候我们虽然已经得到目标多项式。但其形式并不是我们想要的
所以 之后利用 FFTFFTFFT的逆运算又快速的变换回去
我们记:FFTFFTFFT的逆运算为 FFT−FFT^-FFT−
对于一个nnn次多项A的表示。最常见的形式(系数表达):
A(x)=∑i=0n−1aixiA(x)=\sum_{i=0}^{n-1}a_ix^iA(x)=i=0∑n−1aixi
这里的n次多项式是指最高次项指数为n−1的多项式这里的n次多项式是指最高次项指数为n-1的多项式这里的n次多项式是指最高次项指数为n−1的多项式
小学生手算 系数形式 的 多项式乘法 的 复杂度是 O(n2)O(n^2)O(n2)
如果我们知道了一个nnn次多项式 的曲线上的nnn个不同的的点。
我们是可以计算出来这个多项式的
why?why?why?
把系数看做未知数。列出来n个方程组。解出来这n个系数。
也就是说,给出曲线上的n个点:
<(x0,y0),(x1,y1),(x2,y2),.....(xn−1,yn−1)><(x_0,y_0),(x_1,y_1),(x_2,y_2),.....(x_{n-1},y_{n-1})><(x0,y0),(x1,y1),(x2,y2),.....(xn−1,yn−1)>
我们可以确定其系数形式;
我们称这种表达一个多项式的方法叫:点值表达
他有很多优点。比如可以O(n)O(n)O(n)时间内计算一个多项式乘法
再给出一个多项式:
B(x)=∑i=0n−1bixiB(x)=\sum_{i=0}^{n-1}b_ix^iB(x)=i=0∑n−1bixi
设A(x)A(x)A(x)与B(x)B(x)B(x)在取相同x的点值表达为:
<(x0,A0),(x1,A1),(x2,A2),.....(x2n−2,A2n−2)><(x_0,A_0),(x_1,A_1),(x_2,A_2),.....(x_{2n-2},A_{2n-2})><(x0,A0),(x1,A1),(x2,A2),.....(x2n−2,A2n−2)>
<(x0,B0),(x1,B1),(x2,B2),.....(x2n−2,B2n−2)><(x_0,B_0),(x_1,B_1),(x_2,B_2),.....(x_{2n-2},B_{2n-2})><(x0,B0),(x1,B1),(x2,B2),.....(x2n−2,B2n−2)>
这里Ai=A(xi)A_i=A(x_i)Ai=A(xi) , Bi=B(xi)\ \ \ \ \ B_i=B(x_i) Bi=B(xi)
则A(x)∗B(x)A(x)*B(x)A(x)∗B(x)为:
<(x0,A0B0),(x1,A1B1),(x2,A2B2),.....(x2n−2,A2n−2B2n−2)><(x_0,A_0B_0),(x_1,A_1B_1),(x_2,A_2B_2),.....(x_{2n-2},A_{2n-2}B_{2n-2})><(x0,A0B0),(x1,A1B1),(x2,A2B2),.....(x2n−2,A2n−2B2n−2)>
非常方便。顺其自然;
之所以A(x)A(x)A(x)与B(x)B(x)B(x)多取一些点。是因为A(x)∗B(x)A(x)*B(x)A(x)∗B(x)的次数界增加。
取点不足会导次数界达不到2n−12n-12n−1
对于一个nnn次多项式。随机求其n个不同的点的朴素方法复杂度是O(n2)O(n^2)O(n2)
假设 n为偶数。那么我们把A(x)A(x)A(x)。重组为两个多项式:
设其系数的下标为偶数组成的多项式为A[0](x)A^{[0]}(x)A[0](x):
A[0](x)=a0+a2x+a4x2+...+an−2xn2−1A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{\frac{n}{2}-1}A[0](x)=a0+a2x+a4x2+...+an−2x2n−1
设其系数的下标为奇数组成的多项式为A[1](x)A^{[1]}(x)A[1](x):
A[1](x)=a1+a3x+a5x2+...+an−1xn2−1A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{\frac{n}{2}-1}A[1](x)=a1+a3x+a5x2+...+an−1x2n−1
所以有:
A(xi)=A[0](xi2)+xiA[1](xi2)A(x_i)=A^{[0]}(x_i^2)+x_iA^{[1]}(x_i^2)A(xi)=A[0](xi2)+xiA[1](xi2)
好像并不会优化运算。
因为<x02,x12,....xn−12><x_0^2,x_1^2,....x_{n-1}^2><x02,x12,....xn−12>
依然有可能会组成n个不同的数字。
那么我们要计算的规模不会减半。
如果我们恰当选取一些x。是否会优化运算呢。
例如:
<x0,x1,....xn2−1,−x0,−x1,....−xn2−1><x_0,x_1,....x_{\frac{n}{2}-1},-x_0,-x_1,....-x_{\frac{n}{2}-1}><x0,x1,....x2n−1,−x0,−x1,....−x2n−1>
各个数字平方后。得到的集合大小减半:
<x02,x12,....xn2−12><x_0^2,x_1^2,....x_{\frac{n}{2}-1}^2><x02,x12,....x2n−12>
因为(−x)2=x2(-x)^2=x^2(−x)2=x2
那么A(−x)=A[0](x2)−xA[1](x2)A(x)=A[0](x2)+xA[1](x2)A(-x)=A^{[0]}(x^2)-xA^{[1]}(x^2)\\A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)A(−x)=A[0](x2)−xA[1](x2)A(x)=A[0](x2)+xA[1](x2)
但是这种关系不会传递下去。平方后得到的数字全是不小于0的数字。
再次递归又回到了原来的形式。很尴尬。
不过这启迪我们。取相反数。然后平放。可以把问题规模减半。但再一次递归就失效了。是否存在不失效的取值呢。
如果对于一个有偶数个元素的数字集合。每个元素平方后。得到的新集合。去除重复元素后。集合大小能够减半。并且得到的新集合如果为偶数。新集合依然满足上面的性质。我们称这个集合有折半的性质。
如果我们可以快速的找到一个满足折半性质的自变量xxx的取值集合.
分治就是可行的。
有一种东西。可以满足我们的需求。
那就是 。n次单位复数根:
(使用复数确实让人有顾虑。这也是我学习FFT时最大 思想负担)
我们定义i=−1i=\sqrt{-1}i=−1
那么形如cosθ+i∗sinθcos\theta +i*sin\thetacosθ+i∗sinθ
的复数有诸多良好的性质。
例如:(cosθ+i∗sinθ)k=cos kθ+i∗sin kθ(cos\theta +i*sin\theta)^k=cos\ k\theta +i*sin\ k\theta(cosθ+i∗sinθ)k=cos kθ+i∗sin kθ
上面的性质可以用数学归纳法得到(较为简单。这里不做证明
我们记:
Wn=cos2πn +i∗sin2πnW_n=cos\frac{2\pi}{n} +i*sin\frac{2\pi}{n}Wn=cosn2π +i∗sinn2π
我会告诉你 多项式在x取下面的值时。有助于我们进行DFT
<Wn0,Wn1,.....Wnn−1><W_n^0,W_n^1,.....W_n^{n-1}><Wn0,Wn1,.....Wnn−1>
上面的n个复数称为。n次单位复数根。(n次单位复数根是指这n个数)
如果我们x取WnkW_n^kWnk时。根据刚才的分解:
A(Wnk)=A[0]( (Wnk)2)+WnkA[1]((Wnk)2)A(W_n^k)=A^{[0]}(\ (W_n^k)^2)+W_n^kA^{[1]}((W_n^k)^2)A(Wnk)=A[0]( (Wnk)2)+WnkA[1]((Wnk)2)
而(Wnk)2=Wn2k=cos 2k∗2πn+i∗sin2k∗2πn=cos 2πkn2+i∗sin2πkn2=Wn2k(W_n^k)^2=W_n^{2k}=cos\ \frac{2k*2\pi}{n}+i*sin\frac{2k*2\pi}{n}\\=cos\ \frac{2\pi k}{\frac{n}{2}}+i*sin\frac{2\pi k}{\frac{n}{2}}=W_{\frac{n}{2}}^k(Wnk)2=Wn2k=cos n2k∗2π+i∗sinn2k∗2π=cos 2n2πk+i∗sin2n2πk=W2nk
因为三角函数的周期性:
我们有:
Wnk=cosk∗2πn+i∗sink∗2πn=cos(k mod n)∗2πn+i∗sin(k mod n)∗2πn=Wnk mod nW_n^k=cos\frac{k*2\pi}{n}+i*sin\frac{k*2\pi}{n}\\=cos\frac{(k\ mod\ n)*2\pi}{n}+i*sin\frac{(k\ mod\ n)*2\pi}{n}\\=W_n^{k\ mod \ n}Wnk=cosnk∗2π+i∗sinnk∗2π=cosn(k mod n)∗2π+i∗sinn(k mod n)∗2π=Wnk mod n
所以:(Wnk)2=Wn2k=Wn2k mod n2(W_n^k)^2=W_{\frac{n}{2}}^k=W_{\frac{n}{2}}^{k\ mod\ \frac{n}{2}}(Wnk)2=W2nk=W2nk mod 2n
则:<(Wn0)2,(Wn1)2,...(Wnn−1)2><(W_n^0)^2,(W_n^1)^2,...(W_n^{n-1})^2><(Wn0)2,(Wn1)2,...(Wnn−1)2>
等效于<Wn20,Wn21,...Wn2n2−1,Wn20,Wn21,...Wn2n2−1> <W_\frac{n}{2}^0,W_\frac{n}{2}^1,...W_\frac{n}{2}^{\frac{n}{2}-1},W_\frac{n}{2}^0,W_\frac{n}{2}^1,...W_\frac{n}{2}^{\frac{n}{2}-1}><W2n0,W2n1,...W2n2n−1,W2n0,W2n1,...W2n2n−1>
惊不惊喜,意不意外。
xxx的取值集合取单位复数根。不但满足折半的性质。而且还有一定的规律性。与原集合保持一致。
这意味着我们只需要计算取:
<Wn20,Wn21,...Wn2n2−1><W_\frac{n}{2}^0,W_\frac{n}{2}^1,...W_\frac{n}{2}^{\frac{n}{2}-1}><W2n0,W2n1,...W2n2n−1>
的A[0],A[1]A^{[0]},A^{[1]}A[0],A[1]
问题范围减半。并且如果n为偶数。再次平方依然集合大小减半。
记:Ak[0]=A[0](Wn2k)=A[0]((Wnk)2)Ak[1]=A[1](Wn2k)=A[1]((Wnk)2)A_{k}^{[0]}=A^{[0]}(W_\frac{n}{2}^k)=A^{[0]}((W_n^k)^2)\\A_{k}^{[1]}=A^{[1]}(W_\frac{n}{2}^k)=A^{[1]}((W_n^k)^2)Ak[0]=A[0](W2nk)=A[0]((Wnk)2)Ak[1]=A[1](W2nk)=A[1]((Wnk)2)
那么,当我们得到:
<A0[0],A1[0],..An2−1[0],A0[1],A1[1],...An2−1[1]><A_{0}^{[0]},A_{1}^{[0]},..A_{\frac{n}{2}-1}^{[0]},A_{0}^{[1]},A_{1}^{[1]},...A_{\frac{n}{2}-1}^{[1]}><A0[0],A1[0],..A2n−1[0],A0[1],A1[1],...A2n−1[1]>
这意味着我们得到了所有的A[0]((Wni)2),A[1]((Wni)2)A^{[0]}((W_n^i)^2),A^{[1]}((W_n^i)^2)A[0]((Wni)2),A[1]((Wni)2)
则:
A( Wnk)=Ak modn2[0]+WnkAk mod n2[1]A(\ W_n^k)=A_{k\ mod \frac{n}{2}}^{[0]}+W_n^kA_{k\ mod\ \frac{n}{2}}^{[1]}A( Wnk)=Ak mod2n[0]+WnkAk mod 2n[1]
即得到了:<A0,A1,A2,...An−1><A_0,A_1,A_2,...A_{n-1}><A0,A1,A2,...An−1>
其中Ak=A(Wnk)A_k=A(W_n^k)Ak=A(Wnk)
当然。FFT的计算n要取2的整数次幂。这是因为每次减半。我们可以把不足的系数用0填充。
上面过程可以看作 ,取单位复数根计算目标y向量,如下图:
[(Wn0)0(Wn0)1⋯(Wn0)n−1(Wn1)0(Wn1)1⋯(Wn1)n−1⋮⋮⋱⋮(Wnn−1)0(Wnn−1)1⋯(Wnn−1)n−1][a0a1⋮an−1]=[y0y1⋮yn−1]\left[\begin{matrix} (W_n^0)^0 & (W_n^0)^1 & \cdots & (W_n^0)^{n-1} \\(W_n^1)^0 & (W_n^1)^1 & \cdots & (W_n^1)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (W_n^{n-1})^0 & (W_n^{n-1})^1 & \cdots & (W_n^{n-1})^{n-1} \\\end{matrix}\right]\left[\begin{matrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \\\end{matrix}\right]=\left[\begin{matrix} y_0 \\ y_1 \\ \vdots \\ y_{n-1} \\\end{matrix}\right]⎣⎢⎢⎢⎡(Wn0)0(Wn1)0⋮(Wnn−1)0(Wn0)1(Wn1)1⋮(Wnn−1)1⋯⋯⋱⋯(Wn0)n−1(Wn1)n−1⋮(Wnn−1)n−1⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡a0a1⋮an−1⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡y0y1⋮yn−1⎦⎥⎥⎥⎤
本着尽可能简单的原则。我们不在特别的说这个矩阵。
因为之前说关于系数的n个方程必然可解。所以上面的矩阵:
[(Wn0)0(Wn0)1⋯(Wn0)n−1(Wn1)0(Wn1)1⋯(Wn1)n−1⋮⋮⋱⋮(Wnn−1)0(Wnn−1)1⋯(Wnn−1)n−1]\left[\begin{matrix} (W_n^0)^0 & (W_n^0)^1 & \cdots & (W_n^0)^{n-1} \\(W_n^1)^0 & (W_n^1)^1 & \cdots & (W_n^1)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (W_n^{n-1})^0 & (W_n^{n-1})^1 & \cdots & (W_n^{n-1})^{n-1} \\\end{matrix}\right]⎣⎢⎢⎢⎡(Wn0)0(Wn1)0⋮(Wnn−1)0(Wn0)1(Wn1)1⋮(Wnn−1)1⋯⋯⋱⋯(Wn0)n−1(Wn1)n−1⋮(Wnn−1)n−1⎦⎥⎥⎥⎤
必然存在逆矩阵。(有线性代数基础的,或许不陌生)
可能你还是有点迷惑。没关系:
如果我们把FFT看作计算上面特别的矩阵相乘的算法呢。
利用FFT我们可以快速得到:
[(Wn0)0(Wn0)1⋯(Wn0)n−1(Wn1)0(Wn1)1⋯(Wn1)n−1⋮⋮⋱⋮(Wnn−1)0(Wnn−1)1⋯(Wnn−1)n−1][a0a1⋮an−1]−(FFT)−>[A0A1⋮An−1]\left[\begin{matrix} (W_n^0)^0 & (W_n^0)^1 & \cdots & (W_n^0)^{n-1} \\(W_n^1)^0 & (W_n^1)^1 & \cdots & (W_n^1)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (W_n^{n-1})^0 & (W_n^{n-1})^1 & \cdots & (W_n^{n-1})^{n-1} \\\end{matrix}\right]\left[\begin{matrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \\\end{matrix}\right]-(FFT)->\left[\begin{matrix} A_0 \\ A_1 \\ \vdots \\ A_{n-1} \\\end{matrix}\right]⎣⎢⎢⎢⎡(Wn0)0(Wn1)0⋮(Wnn−1)0(Wn0)1(Wn1)1⋮(Wnn−1)1⋯⋯⋱⋯(Wn0)n−1(Wn1)n−1⋮(Wnn−1)n−1⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡a0a1⋮an−1⎦⎥⎥⎥⎤−(FFT)−>⎣⎢⎢⎢⎡A0A1⋮An−1⎦⎥⎥⎥⎤
同时 。我可以直接告诉你:D=[(Wn0)0(Wn0)1⋯(Wn0)n−1(Wn1)0(Wn1)1⋯(Wn1)n−1⋮⋮⋱⋮(Wnn−1)0(Wnn−1)1⋯(Wnn−1)n−1]D=\left[\begin{matrix} (W_n^0)^0 & (W_n^0)^1 & \cdots & (W_n^0)^{n-1} \\(W_n^1)^0 & (W_n^1)^1 & \cdots & (W_n^1)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (W_n^{n-1})^0 & (W_n^{n-1})^1 & \cdots & (W_n^{n-1})^{n-1} \\\end{matrix}\right]D=⎣⎢⎢⎢⎡(Wn0)0(Wn1)0⋮(Wnn−1)0(Wn0)1(Wn1)1⋮(Wnn−1)1⋯⋯⋱⋯(Wn0)n−1(Wn1)n−1⋮(Wnn−1)n−1⎦⎥⎥⎥⎤
的逆矩阵为:
E=[1n(Wn−0)01n(Wn−0)1⋯1n(Wn−0)n−11n(Wn−1)01n(Wn−1)1⋯1n(Wn−1)n−1⋮⋮⋱⋮1n(Wn−(n−1))01n(Wn−(n−1))1⋯1n(Wn−(n−1))n−1]E=\left[\begin{matrix} \frac{1}{n}(W_n^{-0})^0 & \frac{1}{n}(W_n^{-0})^1 & \cdots & \frac{1}{n}(W_n^{-0})^{n-1} \\\frac{1}{n}(W_n^{-1})^0 &\frac{1}{n} (W_n^{-1})^1 & \cdots &\frac{1}{n} (W_n^{-1})^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{1}{n}(W_n^{-(n-1)})^0 &\frac{1}{n} (W_n^{-(n-1)})^1 & \cdots & \frac{1}{n}(W_n^{-(n-1)})^{n-1} \\\end{matrix}\right]E=⎣⎢⎢⎢⎡n1(Wn−0)0n1(Wn−1)0⋮n1(Wn−(n−1))0n1(Wn−0)1n1(Wn−1)1⋮n1(Wn−(n−1))1⋯⋯⋱⋯n1(Wn−0)n−1n1(Wn−1)n−1⋮n1(Wn−(n−1))n−1⎦⎥⎥⎥⎤
因为我直接给你了答案。。。所以怀疑的话。我们可以计算一下看看。
我们记Y=E*D
Y[t][j]=∑k=0n−1E[t][k]∗D[k][j], 其中t,j∈ [0,n−1]Y[t][j]=\sum_{k=0}^{n-1}E[t][k]*D[k][j],\ \ \ \ 其中t,j\in\ [0,n-1]Y[t][j]=k=0∑n−1E[t][k]∗D[k][j], 其中t,j∈ [0,n−1]
Y[t][j]=∑k=0n−11nWn−tk∗Wnkj=1n∑k=0n−1(Wnj−t)kY[t][j]=\sum_{k=0}^{n-1}\frac{1}{n}W_n^{-tk}*W_n^{kj}=\frac{1}{n}\sum_{k=0}^{n-1}(W_n^{j-t})^kY[t][j]=k=0∑n−1n1Wn−tk∗Wnkj=n1k=0∑n−1(Wnj−t)k
上面可以看作等比数列求和。
当t==jt==jt==j时:Y[t][j]=1n∑k=0n−11k=1Y[t][j]=\frac{1}{n}\sum_{k=0}^{n-1}1^k=1Y[t][j]=n1k=0∑n−11k=1
当t≠jt\ne jt̸=j时,令u=j−tu=j-tu=j−t:
Y[t][j]=1n∑k=0n−1(Wnu)k=1n∗(Wnu)n−1Wnu−1 (等比数列公式)Y[t][j]=\frac{1}{n}\sum_{k=0}^{n-1}(W_n^{u})^k\\=\frac{1}{n}*\frac{(W_n^u)^n-1}{W_n^u-1}\ \ \ \ \ \ (等比数列公式) Y[t][j]=n1k=0∑n−1(Wnu)k=n1∗Wnu−1(Wnu)n−1 (等比数列公式)
因为:(Wnu)n=Wnun=Wnun mod n=Wn0=1(W_n^u)^n=W_n^{un}=W_n^{un\ mod\ n}=W_n^0=1(Wnu)n=Wnun=Wnun mod n=Wn0=1
所以,t≠jt\ne jt̸=j时:Y[t][j]=0Y[t][j]=0Y[t][j]=0
所以。Y为单位矩阵。E,D互为逆矩阵得证。
所以:
[(Wn0)0(Wn0)1⋯(Wn0)n−1(Wn−1)0(Wn−1)1⋯(Wn−1)n−1⋮⋮⋱⋮(Wn−(n−1))0(Wn−(n−1))1⋯(Wn−(n−1))n−1][A0A1⋮An−1]−(FFT)−>[na0na1⋮nan−1]\left[\begin{matrix} (W_n^0)^0 & (W_n^0)^1 & \cdots & (W_n^0)^{n-1} \\(W_n^{-1})^0 & (W_n^{-1})^1 & \cdots & (W_n^{-1})^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (W_n^{-(n-1)})^0 & (W_n^{-(n-1)})^1 & \cdots & (W_n^{-(n-1)})^{n-1} \\\end{matrix}\right]\left[\begin{matrix} A_0 \\ A_1 \\ \vdots \\ A_{n-1} \\\end{matrix}\right]-{(FFT)}->\left[\begin{matrix} na_0 \\ na_1 \\ \vdots \\ na_{n-1} \\\end{matrix}\right]⎣⎢⎢⎢⎡(Wn0)0(Wn−1)0⋮(Wn−(n−1))0(Wn0)1(Wn−1)1⋮(Wn−(n−1))1⋯⋯⋱⋯(Wn0)n−1(Wn−1)n−1⋮(Wn−(n−1))n−1⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡A0A1⋮An−1⎦⎥⎥⎥⎤−(FFT)−>⎣⎢⎢⎢⎡na0na1⋮nan−1⎦⎥⎥⎥⎤
所以。FFT−FFT^-FFT−其实就是FFTFFTFFT的WnkW_n^kWnk变为Wn−kW_n^{-k}Wn−k得到。
并且所得目标向量每个元素除掉n就可以了
FFT的高效实现:
通常。我们希望FFT的实现尽可能快。所以FFT算法也都使用迭代结构而不是递归结构。
下面介绍一种常用的去递归方法。
首先。上面介绍的FFT是只能处理2的整数次幂的次数界的多项式
所以不特别说明。都有n=2kn=2^kn=2k(系数不足的用0补足)
对于递归的第一步:
输入数组<a0,a1,...an−1><a_0,a_1,...a_{n-1}><a0,a1,...an−1>
被分解为,偶数和奇数两个部分数
<a0,a2,a4,...an−2> , <a1,a3,a5,...an−1><a_0,a_2,a_4,...a_{n-2}>\ \ \ ,\ \ \ <a_1,a_3,a_5,...a_{n-1}><a0,a2,a4,...an−2> , <a1,a3,a5,...an−1>
下标二进制形式右起第一个字符为 0 被分配到左集合
下标二进制形式右起第一个字符为 1 被分配到右集合
更一般的。对于第k次递归:
下标右起第k个字符为 0 被分配到它所在问题的左集合
下标右起第k个字符为 1 被分配到它所在问题的右集合
那么对于一个2进制形式的数字 B。将其表示为(注意:n=2kn=2^kn=2k):
B=∑t=0k−1bt∗2tB=\sum_{t=0}^{k-1}b_t*2^tB=t=0∑k−1bt∗2t
根据之前分析。这个数字递归后将会被分配到的位置为:
P(B)=∑t=0k−1bt∗2k−1−t, (因为每次都删去一半位置)P(B)=\sum_{t=0}^{k-1}b_t*2^{k-1-t},\ \ \ \ \ (因为每次都删去一半位置)P(B)=t=0∑k−1bt∗2k−1−t, (因为每次都删去一半位置)
所以。对于一个数字的二进制形式的前k位对称过来。就是递归后的位置。
例如:k=4k=4k=4, (0011)2−对称后−>(1100)2(0011)_2-对称后->(1100)_2(0011)2−对称后−>(1100)2
我们调用FFT前对数组进行一次上述重拍。
便可自底向上迭代实现FFT
总结:
由于递归结构的FFT较好理解。而理解递归结构的FFT后。
不难写出非递归结构
所以在这里我们对递归结构的FFT基本框架做一个总结。
因为使用到了复数。所以不可避免的 。我们以下计算都在复数范围的
方便起见。我们用符号:Co(a,b)Co(a,b)Co(a,b)来表示复数a+bia+bia+bi
那么多项式:
A(x)=∑t=0n−1atxt=∑t=0n−1Co(at,0)xtA(x)=\sum_{t=0}^{n-1}a_tx^t=\sum_{t=0}^{n-1}Co(a_t,0)x^tA(x)=t=0∑n−1atxt=t=0∑n−1Co(at,0)xt
Wn=Co(cos2πn,sin2πn)W_n=Co(cos\frac{2\pi}{n},sin\frac{2\pi}{n})Wn=Co(cosn2π,sinn2π)
Wnk=Co(cos2kπn,sin2kπn)W_n^k=Co(cos\frac{2k\pi}{n},sin\frac{2k\pi}{n})Wnk=Co(cosn2kπ,sinn2kπ)
令:Y[t]=Co(at,0) , t∈[0,n−1]Y[t]=Co(a_t,0)\ \ ,\ t\in[0,n-1]Y[t]=Co(at,0) , t∈[0,n−1]
当n=1时:
FFT(Y[],n)=Y[0]FFT(Y[],n)=Y[0]FFT(Y[],n)=Y[0]
n>1时:
令:Y[0][t]=Y[2t]Y^{[0]}[t]=Y[2t]Y[0][t]=Y[2t] , Y[1][t]=Y[2t+1]Y^{[1]}[t]=Y[2t+1]Y[1][t]=Y[2t+1]
计算FFT(Y[0],n2)FFT(Y^{[0]},\frac{n}{2})FFT(Y[0],2n)与FFT(Y[1],n2)FFT(Y^{[1]},\frac{n}{2})FFT(Y[1],2n)
计算完成后,令:
Y[t]=Y[0][t]+WntY[1][t]Y[t+n2]=Y[0][t]−WntY[1][t] (其中t∈[0,n2−1])Y[t]=Y^{[0]}[t]+W_n^tY^{[1]}[t]\\Y[t+\frac{n}{2}]=Y^{[0]}[t]-W_n^tY^{[1]}[t]\ \ \ \ \ \ \\(其中t\in [0,\frac{n}{2}-1])Y[t]=Y[0][t]+WntY[1][t]Y[t+2n]=Y[0][t]−WntY[1][t] (其中t∈[0,2n−1])
返回Y[];Y[];Y[];
------------------------------------------------------------------------
我对上述思想解释一下
A在<Wn0,Wn1,...Wnn−1><W_n^0,W_n^1,...W_{n}^{n-1}><Wn0,Wn1,...Wnn−1>处的值保存在Y中。
并返回给上一层。所以对于当前调用:
A[0](Wn2t)=Y[0][t]A[1](Wn2t)=Y[1][t]A^{[0]}(W_{\frac{n}{2}}^t)=Y^{[0]}[t]\\A^{[1]}(W_{\frac{n}{2}}^t)=Y^{[1]}[t]A[0](W2nt)=Y[0][t]A[1](W2nt)=Y[1][t]
A(Wnt)=A[0]((Wnt)2)+WntA[1]((Wnt)2)=A[0](Wn2t)+WntA[1](Wn2t)A(W_n^t)=A^{[0]}((W_n^t)^2)+W_n^tA^{[1]}((W_n^t)^2)\\=A^{[0]}(W_{\frac{n}{2}}^t)+W_n^tA^{[1]}(W_{\frac{n}{2}}^t)A(Wnt)=A[0]((Wnt)2)+WntA[1]((Wnt)2)=A[0](W2nt)+WntA[1](W2nt)
当t∈[0,n2−1]t\in[0,\frac{n}{2}-1]t∈[0,2n−1] 时:
A(Wnt)=Y[t]=A[0](Wn2t)+WntA[1](Wn2t)=Y[0][t]+WntY[1][t]A(W_n^t)=Y[t]=A^{[0]}(W_{\frac{n}{2}}^t)+W_n^tA^{[1]}(W_{\frac{n}{2}}^t)=Y^{[0]}[t]+W_n^tY^{[1]}[t]A(Wnt)=Y[t]=A[0](W2nt)+WntA[1](W2nt)=Y[0][t]+WntY[1][t]
A(Wnt+n2)=Y[t+n2]=A[0](Wn2t+n2)+Wnt+n2A[1](Wn2t+n2)=A[0](Wn2t)−WntA[1](Wn2t)=Y[0][t]−WntY[1][t]A(W_n^{t+\frac{n}{2}})=Y[t+\frac{n}{2}]=A^{[0]}(W_{\frac{n}{2}}^{t+\frac{n}{2}})+W_n^{t+\frac{n}{2}}A^{[1]}(W_{\frac{n}{2}}^{t+\frac{n}{2}})\\=A^{[0]}(W_{\frac{n}{2}}^t)-W_n^tA^{[1]}(W_{\frac{n}{2}}^t)\\=Y^{[0]}[t]-W_n^tY^{[1]}[t]A(Wnt+2n)=Y[t+2n]=A[0](W2nt+2n)+Wnt+2nA[1](W2nt+2n)=A[0](W2nt)−WntA[1](W2nt)=Y[0][t]−WntY[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,s;
int b=len>>1;
for(i=1,j=len>>1,s=len-1;i<s;i++)
{
if(i<j)swap(y[i],y[j]);
k=b;
while(j>=k)
{
j-=k;
k>>=1;
}
j+=k;
}
}
迭代结构的FFT(
on== 1 时: FFTFFTFFT
on==-1时: FFT−FFT^-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,s=h>>1;j<len;j+=h,s+=h)
{
Complex w(1,0);
for(int k=j,b=s;k<s;k++,b++)
{
Complex u=y[k];
Complex t=w*y[b];
y[k]=u+t;
y[b]=u-t;
w=w*wn;
}
}
}
if(on==-1)
for(int i=0;i<len;i++)
y[i].x/=len;
}