多项式与快速傅立叶变换

文章写的有点急。有错误的地方望指出
我学习 FFT 是一个比较慢的过程。 期间反反复复。 我写这篇博文只是一个非常浅显的理解。同时也可以帮助初学者在学习FFT的时候。有所偏重。避免太多思维上的负担。

直接正题吧:

首先: D F T DFT DFT本身并不负责多项式之间的乘法。
D F T DFT DFT只是一种变换。
F F T FFT FFT则是DFT的快速算法。(分治提高效率)

利用 F F T FFT FFT 。我们快速的将多项式变换为利于计算的形式
用这种方便计算的形式计算出来两个多项式的乘积。
这时候我们虽然已经得到目标多项式。但其形式并不是我们想要的
所以 之后利用 F F T FFT FFT的逆运算又快速的变换回去

我们记: F F T FFT FFT的逆运算为 F F T − FFT^- FFT

对于一个 n n n次多项A的表示。最常见的形式(系数表达):

A ( x ) = ∑ i = 0 n − 1 a i x i A(x)=\sum_{i=0}^{n-1}a_ix^i A(x)=i=0n1aixi

这 里 的 n 次 多 项 式 是 指 最 高 次 项 指 数 为 n − 1 的 多 项 式 这里的n次多项式是指最高次项指数为n-1的多项式 nn1
小学生手算 系数形式 的 多项式乘法 的 复杂度是 O ( n 2 ) O(n^2) O(n2)
如果我们知道了一个 n n n次多项式 的曲线上的 n n n个不同的的点。
我们是可以计算出来这个多项式的
w h y ? why? why
把系数看做未知数。列出来n个方程组。解出来这n个系数。
也就是说,给出曲线上的n个点:
&lt; ( x 0 , y 0 ) , ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . . . ( x n − 1 , y n − 1 ) &gt; &lt;(x_0,y_0),(x_1,y_1),(x_2,y_2),.....(x_{n-1},y_{n-1})&gt; <(x0,y0),(x1,y1),(x2,y2),.....(xn1,yn1)>
我们可以确定其系数形式;
我们称这种表达一个多项式的方法叫:点值表达
他有很多优点。比如可以 O ( n ) O(n) O(n)时间内计算一个多项式乘法
再给出一个多项式:

B ( x ) = ∑ i = 0 n − 1 b i x i B(x)=\sum_{i=0}^{n-1}b_ix^i B(x)=i=0n1bixi

A ( x ) A(x) A(x) B ( x ) B(x) B(x)在取相同x的点值表达为:
&lt; ( x 0 , A 0 ) , ( x 1 , A 1 ) , ( x 2 , A 2 ) , . . . . . ( x 2 n − 2 , A 2 n − 2 ) &gt; &lt;(x_0,A_0),(x_1,A_1),(x_2,A_2),.....(x_{2n-2},A_{2n-2})&gt; <(x0,A0),(x1,A1),(x2,A2),.....(x2n2,A2n2)>
&lt; ( x 0 , B 0 ) , ( x 1 , B 1 ) , ( x 2 , B 2 ) , . . . . . ( x 2 n − 2 , B 2 n − 2 ) &gt; &lt;(x_0,B_0),(x_1,B_1),(x_2,B_2),.....(x_{2n-2},B_{2n-2})&gt; <(x0,B0),(x1,B1),(x2,B2),.....(x2n2,B2n2)>
这里 A i = A ( x i ) A_i=A(x_i) Ai=A(xi) ,       B i = B ( x i ) \ \ \ \ \ B_i=B(x_i)      Bi=B(xi)
A ( x ) ∗ B ( x ) A(x)*B(x) A(x)B(x)为:
&lt; ( x 0 , A 0 B 0 ) , ( x 1 , A 1 B 1 ) , ( x 2 , A 2 B 2 ) , . . . . . ( x 2 n − 2 , A 2 n − 2 B 2 n − 2 ) &gt; &lt;(x_0,A_0B_0),(x_1,A_1B_1),(x_2,A_2B_2),.....(x_{2n-2},A_{2n-2}B_{2n-2})&gt; <(x0,A0B0),(x1,A1B1),(x2,A2B2),.....(x2n2,A2n2B2n2)>
非常方便。顺其自然;
之所以 A ( x ) A(x) A(x) B ( x ) B(x) B(x)多取一些点。是因为 A ( x ) ∗ B ( x ) A(x)*B(x) A(x)B(x)的次数界增加。
取点不足会导次数界达不到 2 n − 1 2n-1 2n1
对于一个 n n n次多项式。随机求其n个不同的点的朴素方法复杂度是 O ( n 2 ) O(n^2) O(n2)

假设 n为偶数。那么我们把 A ( x ) A(x) A(x)。重组为两个多项式:
设其系数的下标为偶数组成的多项式为 A [ 0 ] ( x ) A^{[0]}(x) A[0](x)
A [ 0 ] ( x ) = a 0 + a 2 x + a 4 x 2 + . . . + a n − 2 x n 2 − 1 A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{\frac{n}{2}-1} A[0](x)=a0+a2x+a4x2+...+an2x2n1
设其系数的下标为奇数组成的多项式为 A [ 1 ] ( x ) A^{[1]}(x) A[1](x)
A [ 1 ] ( x ) = a 1 + a 3 x + a 5 x 2 + . . . + a n − 1 x n 2 − 1 A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{\frac{n}{2}-1} A[1](x)=a1+a3x+a5x2+...+an1x2n1
所以有:
A ( x i ) = A [ 0 ] ( x i 2 ) + x i A [ 1 ] ( x i 2 ) A(x_i)=A^{[0]}(x_i^2)+x_iA^{[1]}(x_i^2) A(xi)=A[0](xi2)+xiA[1](xi2)
好像并不会优化运算。
因为 &lt; x 0 2 , x 1 2 , . . . . x n − 1 2 &gt; &lt;x_0^2,x_1^2,....x_{n-1}^2&gt; <x02,x12,....xn12>
依然有可能会组成n个不同的数字。
那么我们要计算的规模不会减半。
如果我们恰当选取一些x。是否会优化运算呢。
例如:
&lt; x 0 , x 1 , . . . . x n 2 − 1 , − x 0 , − x 1 , . . . . − x n 2 − 1 &gt; &lt;x_0,x_1,....x_{\frac{n}{2}-1},-x_0,-x_1,....-x_{\frac{n}{2}-1}&gt; <x0,x1,....x2n1,x0,x1,....x2n1>
各个数字平方后。得到的集合大小减半:

&lt; x 0 2 , x 1 2 , . . . . x n 2 − 1 2 &gt; &lt;x_0^2,x_1^2,....x_{\frac{n}{2}-1}^2&gt; <x02,x12,....x2n12>

因为 ( − x ) 2 = x 2 (-x)^2=x^2 (x)2=x2
那么 A ( − x ) = A [ 0 ] ( x 2 ) − x A [ 1 ] ( x 2 ) A ( x ) = A [ 0 ] ( x 2 ) + x A [ 1 ] ( x 2 ) 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的数字。
再次递归又回到了原来的形式。很尴尬。
不过这启迪我们。取相反数。然后平放。可以把问题规模减半。但再一次递归就失效了。是否存在不失效的取值呢。

如果对于一个有偶数个元素的数字集合。每个元素平方后。得到的新集合。去除重复元素后。集合大小能够减半。并且得到的新集合如果为偶数。新集合依然满足上面的性质。我们称这个集合有折半的性质。

如果我们可以快速的找到一个满足折半性质的自变量 x x x的取值集合.

分治就是可行的。

有一种东西。可以满足我们的需求。
那就是 。n次单位复数根:
(使用复数确实让人有顾虑。这也是我学习FFT时最大 思想负担)

我们定义 i = − 1 i=\sqrt{-1} i=1
那么形如 c o s θ + i ∗ s i n θ cos\theta +i*sin\theta cosθ+isinθ
的复数有诸多良好的性质。
例如: ( c o s θ + i ∗ s i n θ ) k = c o s   k θ + i ∗ s i n   k θ (cos\theta +i*sin\theta)^k=cos\ k\theta +i*sin\ k\theta (cosθ+isinθ)k=cos kθ+isin kθ
上面的性质可以用数学归纳法得到(较为简单。这里不做证明
我们记:

W n = c o s 2 π n   + i ∗ s i n 2 π n W_n=cos\frac{2\pi}{n}  +i*sin\frac{2\pi}{n} Wn=cosn2π +isinn2π

我会告诉你 多项式在x取下面的值时。有助于我们进行DFT
&lt; W n 0 , W n 1 , . . . . . W n n − 1 &gt; &lt;W_n^0,W_n^1,.....W_n^{n-1}&gt; <Wn0,Wn1,.....Wnn1>
上面的n个复数称为。n次单位复数根。(n次单位复数根是指这n个数)
如果我们x取 W n k W_n^k Wnk时。根据刚才的分解:
A ( W n k ) = A [ 0 ] (   ( W n k ) 2 ) + W n k A [ 1 ] ( ( W n k ) 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)

( W n k ) 2 = W n 2 k = c o s   2 k ∗ 2 π n + i ∗ s i n 2 k ∗ 2 π n = c o s   2 π k n 2 + i ∗ s i n 2 π k n 2 = W n 2 k (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 n2k2π+isinn2k2π=cos 2n2πk+isin2n2πk=W2nk

因为三角函数的周期性:
我们有:

W n k = c o s k ∗ 2 π n + i ∗ s i n k ∗ 2 π n = c o s ( k   m o d   n ) ∗ 2 π n + i ∗ s i n ( k   m o d   n ) ∗ 2 π n = W n k   m o d   n W_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=cosnk2π+isinnk2π=cosn(k mod n)2π+isinn(k mod n)2π=Wnk mod n

所以: ( W n k ) 2 = W n 2 k = W n 2 k   m o d   n 2 (W_n^k)^2=W_{\frac{n}{2}}^k=W_{\frac{n}{2}}^{k\ mod\ \frac{n}{2}} (Wnk)2=W2nk=W2nk mod 2n

则: &lt; ( W n 0 ) 2 , ( W n 1 ) 2 , . . . ( W n n − 1 ) 2 &gt; &lt;(W_n^0)^2,(W_n^1)^2,...(W_n^{n-1})^2&gt; <(Wn0)2,(Wn1)2,...(Wnn1)2>

等效于 &lt; W n 2 0 , W n 2 1 , . . . W n 2 n 2 − 1 , W n 2 0 , W n 2 1 , . . . W n 2 n 2 − 1 &gt; &lt;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}&gt; <W2n0,W2n1,...W2n2n1,W2n0,W2n1,...W2n2n1>

惊不惊喜,意不意外。
x x x的取值集合取单位复数根。不但满足折半的性质。而且还有一定的规律性。与原集合保持一致。
这意味着我们只需要计算取:

&lt; W n 2 0 , W n 2 1 , . . . W n 2 n 2 − 1 &gt; &lt;W_\frac{n}{2}^0,W_\frac{n}{2}^1,...W_\frac{n}{2}^{\frac{n}{2}-1}&gt; <W2n0,W2n1,...W2n2n1>

A [ 0 ] , A [ 1 ] A^{[0]},A^{[1]} A[0],A[1]
问题范围减半。并且如果n为偶数。再次平方依然集合大小减半。
记: A k [ 0 ] = A [ 0 ] ( W n 2 k ) = A [ 0 ] ( ( W n k ) 2 ) A k [ 1 ] = A [ 1 ] ( W n 2 k ) = A [ 1 ] ( ( W n k ) 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)
那么,当我们得到:

&lt; A 0 [ 0 ] , A 1 [ 0 ] , . . A n 2 − 1 [ 0 ] , A 0 [ 1 ] , A 1 [ 1 ] , . . . A n 2 − 1 [ 1 ] &gt; &lt;A_{0}^{[0]},A_{1}^{[0]},..A_{\frac{n}{2}-1}^{[0]},A_{0}^{[1]},A_{1}^{[1]},...A_{\frac{n}{2}-1}^{[1]}&gt; <A0[0],A1[0],..A2n1[0],A0[1],A1[1],...A2n1[1]>

这意味着我们得到了所有的 A [ 0 ] ( ( W n i ) 2 ) , A [ 1 ] ( ( W n i ) 2 ) A^{[0]}((W_n^i)^2),A^{[1]}((W_n^i)^2) A[0]((Wni)2),A[1]((Wni)2)
则:
A (   W n k ) = A k   m o d n 2 [ 0 ] + W n k A k   m o d   n 2 [ 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]
即得到了: &lt; A 0 , A 1 , A 2 , . . . A n − 1 &gt; &lt;A_0,A_1,A_2,...A_{n-1}&gt; <A0,A1,A2,...An1>
其中 A k = A ( W n k ) A_k=A(W_n^k) Ak=A(Wnk)
当然。FFT的计算n要取2的整数次幂。这是因为每次减半。我们可以把不足的系数用0填充。
上面过程可以看作 ,取单位复数根计算目标y向量,如下图:
[ ( W n 0 ) 0 ( W n 0 ) 1 ⋯ ( W n 0 ) n − 1 ( W n 1 ) 0 ( W n 1 ) 1 ⋯ ( W n 1 ) n − 1 ⋮ ⋮ ⋱ ⋮ ( W n n − 1 ) 0 ( W n n − 1 ) 1 ⋯ ( W n n − 1 ) n − 1 ] [ a 0 a 1 ⋮ a n − 1 ] = [ y 0 y 1 ⋮ y n − 1 ] \left[\begin{matrix} (W_n^0)^0 &amp; (W_n^0)^1 &amp; \cdots &amp; (W_n^0)^{n-1} \\(W_n^1)^0 &amp; (W_n^1)^1 &amp; \cdots &amp; (W_n^1)^{n-1} \\ \vdots &amp; \vdots &amp; \ddots &amp; \vdots \\ (W_n^{n-1})^0 &amp; (W_n^{n-1})^1 &amp; \cdots &amp; (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(Wnn1)0(Wn0)1(Wn1)1(Wnn1)1(Wn0)n1(Wn1)n1(Wnn1)n1a0a1an1=y0y1yn1
本着尽可能简单的原则。我们不在特别的说这个矩阵。
因为之前说关于系数的n个方程必然可解。所以上面的矩阵:
[ ( W n 0 ) 0 ( W n 0 ) 1 ⋯ ( W n 0 ) n − 1 ( W n 1 ) 0 ( W n 1 ) 1 ⋯ ( W n 1 ) n − 1 ⋮ ⋮ ⋱ ⋮ ( W n n − 1 ) 0 ( W n n − 1 ) 1 ⋯ ( W n n − 1 ) n − 1 ] \left[\begin{matrix} (W_n^0)^0 &amp; (W_n^0)^1 &amp; \cdots &amp; (W_n^0)^{n-1} \\(W_n^1)^0 &amp; (W_n^1)^1 &amp; \cdots &amp; (W_n^1)^{n-1} \\ \vdots &amp; \vdots &amp; \ddots &amp; \vdots \\ (W_n^{n-1})^0 &amp; (W_n^{n-1})^1 &amp; \cdots &amp; (W_n^{n-1})^{n-1} \\\end{matrix}\right] (Wn0)0(Wn1)0(Wnn1)0(Wn0)1(Wn1)1(Wnn1)1(Wn0)n1(Wn1)n1(Wnn1)n1
必然存在逆矩阵。(有线性代数基础的,或许不陌生)
可能你还是有点迷惑。没关系:
如果我们把FFT看作计算上面特别的矩阵相乘的算法呢。
利用FFT我们可以快速得到:
[ ( W n 0 ) 0 ( W n 0 ) 1 ⋯ ( W n 0 ) n − 1 ( W n 1 ) 0 ( W n 1 ) 1 ⋯ ( W n 1 ) n − 1 ⋮ ⋮ ⋱ ⋮ ( W n n − 1 ) 0 ( W n n − 1 ) 1 ⋯ ( W n n − 1 ) n − 1 ] [ a 0 a 1 ⋮ a n − 1 ] − ( F F T ) − &gt; [ A 0 A 1 ⋮ A n − 1 ] \left[\begin{matrix} (W_n^0)^0 &amp; (W_n^0)^1 &amp; \cdots &amp; (W_n^0)^{n-1} \\(W_n^1)^0 &amp; (W_n^1)^1 &amp; \cdots &amp; (W_n^1)^{n-1} \\ \vdots &amp; \vdots &amp; \ddots &amp; \vdots \\ (W_n^{n-1})^0 &amp; (W_n^{n-1})^1 &amp; \cdots &amp; (W_n^{n-1})^{n-1} \\\end{matrix}\right]\left[\begin{matrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \\\end{matrix}\right]-(FFT)-&gt;\left[\begin{matrix} A_0 \\ A_1 \\ \vdots \\ A_{n-1} \\\end{matrix}\right] (Wn0)0(Wn1)0(Wnn1)0(Wn0)1(Wn1)1(Wnn1)1(Wn0)n1(Wn1)n1(Wnn1)n1a0a1an1(FFT)>A0A1An1
同时 。我可以直接告诉你: D = [ ( W n 0 ) 0 ( W n 0 ) 1 ⋯ ( W n 0 ) n − 1 ( W n 1 ) 0 ( W n 1 ) 1 ⋯ ( W n 1 ) n − 1 ⋮ ⋮ ⋱ ⋮ ( W n n − 1 ) 0 ( W n n − 1 ) 1 ⋯ ( W n n − 1 ) n − 1 ] D=\left[\begin{matrix} (W_n^0)^0 &amp; (W_n^0)^1 &amp; \cdots &amp; (W_n^0)^{n-1} \\(W_n^1)^0 &amp; (W_n^1)^1 &amp; \cdots &amp; (W_n^1)^{n-1} \\ \vdots &amp; \vdots &amp; \ddots &amp; \vdots \\ (W_n^{n-1})^0 &amp; (W_n^{n-1})^1 &amp; \cdots &amp; (W_n^{n-1})^{n-1} \\\end{matrix}\right] D=(Wn0)0(Wn1)0(Wnn1)0(Wn0)1(Wn1)1(Wnn1)1(Wn0)n1(Wn1)n1(Wnn1)n1
的逆矩阵为:
E = [ 1 n ( W n − 0 ) 0 1 n ( W n − 0 ) 1 ⋯ 1 n ( W n − 0 ) n − 1 1 n ( W n − 1 ) 0 1 n ( W n − 1 ) 1 ⋯ 1 n ( W n − 1 ) n − 1 ⋮ ⋮ ⋱ ⋮ 1 n ( W n − ( n − 1 ) ) 0 1 n ( W n − ( n − 1 ) ) 1 ⋯ 1 n ( W n − ( n − 1 ) ) n − 1 ] E=\left[\begin{matrix} \frac{1}{n}(W_n^{-0})^0 &amp; \frac{1}{n}(W_n^{-0})^1 &amp; \cdots &amp; \frac{1}{n}(W_n^{-0})^{n-1} \\\frac{1}{n}(W_n^{-1})^0 &amp;\frac{1}{n} (W_n^{-1})^1 &amp; \cdots &amp;\frac{1}{n} (W_n^{-1})^{n-1} \\ \vdots &amp; \vdots &amp; \ddots &amp; \vdots \\ \frac{1}{n}(W_n^{-(n-1)})^0 &amp;\frac{1}{n} (W_n^{-(n-1)})^1 &amp; \cdots &amp; \frac{1}{n}(W_n^{-(n-1)})^{n-1} \\\end{matrix}\right] E=n1(Wn0)0n1(Wn1)0n1(Wn(n1))0n1(Wn0)1n1(Wn1)1n1(Wn(n1))1n1(Wn0)n1n1(Wn1)n1n1(Wn(n1))n1
因为我直接给你了答案。。。所以怀疑的话。我们可以计算一下看看。
我们记Y=E*D
Y [ t ] [ j ] = ∑ k = 0 n − 1 E [ 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=0n1E[t][k]D[k][j],    t,j [0,n1]
Y [ t ] [ j ] = ∑ k = 0 n − 1 1 n W n − t k ∗ W n k j = 1 n ∑ k = 0 n − 1 ( W n j − t ) k Y[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})^k Y[t][j]=k=0n1n1WntkWnkj=n1k=0n1(Wnjt)k
上面可以看作等比数列求和。
t = = j t==j t==j时: Y [ t ] [ j ] = 1 n ∑ k = 0 n − 1 1 k = 1 Y[t][j]=\frac{1}{n}\sum_{k=0}^{n-1}1^k=1 Y[t][j]=n1k=0n11k=1
t ≠ j t\ne j t̸=j时,令 u = j − t u=j-t u=jt
Y [ t ] [ j ] = 1 n ∑ k = 0 n − 1 ( W n u ) k = 1 n ∗ ( W n u ) n − 1 W n u − 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=0n1(Wnu)k=n1Wnu1(Wnu)n1      
因为: ( W n u ) n = W n u n = W n u n   m o d   n = W n 0 = 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 ≠ j t\ne j t̸=j时: Y [ t ] [ j ] = 0 Y[t][j]=0 Y[t][j]=0
所以。Y为单位矩阵。E,D互为逆矩阵得证。
所以:
[ ( W n 0 ) 0 ( W n 0 ) 1 ⋯ ( W n 0 ) n − 1 ( W n − 1 ) 0 ( W n − 1 ) 1 ⋯ ( W n − 1 ) n − 1 ⋮ ⋮ ⋱ ⋮ ( W n − ( n − 1 ) ) 0 ( W n − ( n − 1 ) ) 1 ⋯ ( W n − ( n − 1 ) ) n − 1 ] [ A 0 A 1 ⋮ A n − 1 ] − ( F F T ) − &gt; [ n a 0 n a 1 ⋮ n a n − 1 ] \left[\begin{matrix} (W_n^0)^0 &amp; (W_n^0)^1 &amp; \cdots &amp; (W_n^0)^{n-1} \\(W_n^{-1})^0 &amp; (W_n^{-1})^1 &amp; \cdots &amp; (W_n^{-1})^{n-1} \\ \vdots &amp; \vdots &amp; \ddots &amp; \vdots \\ (W_n^{-(n-1)})^0 &amp; (W_n^{-(n-1)})^1 &amp; \cdots &amp; (W_n^{-(n-1)})^{n-1} \\\end{matrix}\right]\left[\begin{matrix} A_0 \\ A_1 \\ \vdots \\ A_{n-1} \\\end{matrix}\right]-{(FFT)}-&gt;\left[\begin{matrix} na_0 \\ na_1 \\ \vdots \\ na_{n-1} \\\end{matrix}\right] (Wn0)0(Wn1)0(Wn(n1))0(Wn0)1(Wn1)1(Wn(n1))1(Wn0)n1(Wn1)n1(Wn(n1))n1A0A1An1(FFT)>na0na1nan1
所以。 F F T − FFT^- FFT其实就是 F F T FFT FFT W n k W_n^k Wnk变为 W n − k W_n^{-k} Wnk得到。
并且所得目标向量每个元素除掉n就可以了

FFT的高效实现:

通常。我们希望FFT的实现尽可能快。所以FFT算法也都使用迭代结构而不是递归结构。
下面介绍一种常用的去递归方法。
首先。上面介绍的FFT是只能处理2的整数次幂的次数界的多项式
所以不特别说明。都有 n = 2 k n=2^k n=2k(系数不足的用0补足)
对于递归的第一步:
输入数组 &lt; a 0 , a 1 , . . . a n − 1 &gt; &lt;a_0,a_1,...a_{n-1}&gt; <a0,a1,...an1>
被分解为,偶数和奇数两个部分数
&lt; a 0 , a 2 , a 4 , . . . a n − 2 &gt;     ,     &lt; a 1 , a 3 , a 5 , . . . a n − 1 &gt; &lt;a_0,a_2,a_4,...a_{n-2}&gt;\ \ \ ,\ \ \ &lt;a_1,a_3,a_5,...a_{n-1}&gt; <a0,a2,a4,...an2>   ,   <a1,a3,a5,...an1>
下标二进制形式右起第一个字符为 0 被分配到左集合
下标二进制形式右起第一个字符为 1 被分配到右集合
更一般的。对于第k次递归:
下标右起第k个字符为 0 被分配到它所在问题的左集合
下标右起第k个字符为 1 被分配到它所在问题的右集合
那么对于一个2进制形式的数字 B。将其表示为(注意: n = 2 k n=2^k n=2k):

B = ∑ t = 0 k − 1 b t ∗ 2 t B=\sum_{t=0}^{k-1}b_t*2^t B=t=0k1bt2t

根据之前分析。这个数字递归后将会被分配到的位置为:

P ( B ) = ∑ t = 0 k − 1 b t ∗ 2 k − 1 − t ,       ( 因 为 每 次 都 删 去 一 半 位 置 ) P(B)=\sum_{t=0}^{k-1}b_t*2^{k-1-t},\ \ \ \ \ (因为每次都删去一半位置) P(B)=t=0k1bt2k1t,     ()

所以。对于一个数字的二进制形式的前k位对称过来。就是递归后的位置。
例如: k = 4 k=4 k=4, ( 0011 ) 2 − 对 称 后 − &gt; ( 1100 ) 2 (0011)_2-对称后-&gt;(1100)_2 (0011)2>(1100)2
我们调用FFT前对数组进行一次上述重拍。
便可自底向上迭代实现FFT

总结:

由于递归结构的FFT较好理解。而理解递归结构的FFT后。

不难写出非递归结构

所以在这里我们对递归结构的FFT基本框架做一个总结。

因为使用到了复数。所以不可避免的 。我们以下计算都在复数范围的

方便起见。我们用符号: C o ( a , b ) Co(a,b) Co(a,b)来表示复数 a + b i a+bi a+bi
那么多项式:

A ( x ) = ∑ t = 0 n − 1 a t x t = ∑ t = 0 n − 1 C o ( a t , 0 ) x t A(x)=\sum_{t=0}^{n-1}a_tx^t=\sum_{t=0}^{n-1}Co(a_t,0)x^t A(x)=t=0n1atxt=t=0n1Co(at,0)xt

W n = C o ( c o s 2 π n , s i n 2 π n ) W_n=Co(cos\frac{2\pi}{n},sin\frac{2\pi}{n}) Wn=Co(cosn2π,sinn2π)

W n k = C o ( c o s 2 k π n , s i n 2 k π n ) W_n^k=Co(cos\frac{2k\pi}{n},sin\frac{2k\pi}{n}) Wnk=Co(cosn2kπ,sinn2kπ)

令: Y [ t ] = C o ( a t , 0 )    ,   t ∈ [ 0 , n − 1 ] Y[t]=Co(a_t,0)\ \ ,\ t\in[0,n-1] Y[t]=Co(at,0)  , t[0,n1]

当n=1时:

F F T ( Y [ ] , n ) = Y [ 0 ] FFT(Y[],n)=Y[0] FFT(Y[],n)=Y[0]

n>1时:

令: Y [ 0 ] [ t ] = Y [ 2 t ] Y^{[0]}[t]=Y[2t] Y[0][t]=Y[2t] , Y [ 1 ] [ t ] = Y [ 2 t + 1 ] Y^{[1]}[t]=Y[2t+1] Y[1][t]=Y[2t+1]

计算 F F T ( Y [ 0 ] , n 2 ) FFT(Y^{[0]},\frac{n}{2}) FFT(Y[0],2n) F F T ( Y [ 1 ] , n 2 ) FFT(Y^{[1]},\frac{n}{2}) FFT(Y[1],2n)

计算完成后,令:

Y [ t ] = Y [ 0 ] [ t ] + W n t Y [ 1 ] [ t ] Y [ t + n 2 ] = Y [ 0 ] [ t ] − W n t Y [ 1 ] [ t ]        ( 其 中 t ∈ [ 0 , n 2 − 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,2n1]

返回 Y [ ] ; Y[]; Y[];

------------------------------------------------------------------------

我对上述思想解释一下

A在 &lt; W n 0 , W n 1 , . . . W n n − 1 &gt; &lt;W_n^0,W_n^1,...W_{n}^{n-1}&gt; <Wn0,Wn1,...Wnn1>处的值保存在Y中。

并返回给上一层。所以对于当前调用:

A [ 0 ] ( W n 2 t ) = Y [ 0 ] [ t ] A [ 1 ] ( W n 2 t ) = 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 ( W n t ) = A [ 0 ] ( ( W n t ) 2 ) + W n t A [ 1 ] ( ( W n t ) 2 ) = A [ 0 ] ( W n 2 t ) + W n t A [ 1 ] ( W n 2 t ) 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 , n 2 − 1 ] t\in[0,\frac{n}{2}-1] t[0,2n1] 时:

A ( W n t ) = Y [ t ] = A [ 0 ] ( W n 2 t ) + W n t A [ 1 ] ( W n 2 t ) = Y [ 0 ] [ t ] + W n t Y [ 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 ( W n t + n 2 ) = Y [ t + n 2 ] = A [ 0 ] ( W n 2 t + n 2 ) + W n t + n 2 A [ 1 ] ( W n 2 t + n 2 ) = A [ 0 ] ( W n 2 t ) − W n t A [ 1 ] ( W n 2 t ) = Y [ 0 ] [ t ] − W n t Y [ 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 时: F F T FFT FFT

on==-1时: F F T − 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;
}
评论 18
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值