非正式引出

此部分参考这个视频
此部分可能不够严谨，目的主要是引出FFT的分治思想，请读者见谅

多项式的乘积问题

首先给出一个问题：
::: align-center
给定两个多项式，计算二者的乘积

:::
设计一个算法来解决这个问题

暴力

首先，最暴力的想法是采用乘法分配律
e.g.
::: align-center
\begin{aligned} &(x^2+3x+2)(2x^2+1)\\ &=x^2(2x^2+1)+3x(2x^2+1)+2(2x^2+1)\\ &=... \end{aligned}

:::
如此，我们可以考虑将多项式的前系数按照次数由高到低（或由低到高，下面采用由低到高）储存在两个数组之中

这个算法的复杂度为\mathcal{O}(n^2)

新的表示法

现在我们考虑优化上面的算法，寻求更低的复杂度

我们需要用到下面这个定理
::: align-center
n次多项式可以由n+1个两两不同的点唯一确定

:::
证明为：
给定d次多项式P(x)=p_0(x)+p_1(x)x+\cdots+p_dx^d，以及d+1个点\{(x_0,P(x_0)),\cdots,(x_d,P(x_d))\}
将其写为
::: align-center
\begin{bmatrix} P(x_0)\\ P(x_1)\\ \vdots\\ P(x_d) \end{bmatrix}=\begin{bmatrix} 1&x_0&x_0^2&\cdots&x_0^d\\ 1&x_1&x_1^2&\cdots&x_1^d\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&x_d&x_d^2&\cdots&x_d^d \end{bmatrix} \begin{bmatrix} p_0\\p_1\\\vdots\\p_d \end{bmatrix}

:::
容易知道等号右边的矩阵是一个范德蒙德矩阵，在点两两不相同时它是可逆矩阵，因此此矩阵满秩序，任何点向量施加此矩阵的解向量是唯一的，故此多项式唯一确定

有了这个结论，我们就可以在保持x不变的条件下取多项式上两两不同的点。对于上一节的两个多项式，因为它们的乘积是一个四次多项式，所以我们只需要在两个乘积前的多项式各取5个点，并且计算这两对点的“乘积”得到新的5个点，之后再由这5个点确定唯一的多项式，那么这个新的多项式，应该就是我们需要求的

如此，我们的算法复杂度看起来缩减到了\mathcal{O}(n)

但是如此同时，我们也需要考虑
::: align-center
如何取点？

:::
显然，对于n阶多项式，我们需要取n+1个点，但是若我们采用前面的数组+存储系数，我们取点的复杂度将高达\mathcal{O}(n^2)（n+1个要取的点及内部遍历数组）

但是，我们可以考虑改造多项式，采用分治策略，提出其中奇函数与偶函数的部分，对于其中奇函数的部分，再提出一个x使其变为偶函数
例如
::: align-center
\begin{aligned} &P(x)=3x^5+2x^4+x^3+7x^2+5x+1\\ &=(2x^4+7x^2+1)+x(3x^4+x^2+5)\\ &=P_e[x^2]+xP_o[x^2] \end{aligned}

:::

其中
::: align-center
\begin{aligned} &P_e[x^2]=2(x^2)^2+7(x^2)+1\\ &P_o[x^2]=3(x^2)^2+(x^2)+5 \end{aligned}

:::

如此便可利用函数奇偶性
::: align-center
\begin{aligned} &P(x_i)=P_e[x_i^2]+x_iP_o(x_i^2)\\ &P(-x_i)=P_e[x_i^2]-x_iP_o(x_i^2) \end{aligned}

:::
要取两个点(\pm x_i,P(\pm x_i))只需要计算一遍P_e[x_i^2],P_o[x_i^2]，之后使用上式即可得到函数值

同时，针对P_e,P_o的取值问题又可以被以类似的方法拆分为使用奇偶性解决的问题，如果我们的取值问题就变成了一个递归问题，我们使用
::: align-center
P(x_i)=P_e[x_i^2]+x_iP_o(x_i^2)

:::
来计算x_i处的函数值，对于P_e[x_i^2],P_o[x_i^2]，我们分别令它们为新的P(x_i)，每次计算出在x_i处的函数值后我们都可以直接利用上面奇偶性的式子获得在-x_i处的函数值，第一次递归后获得在两个点x=\pm x_i处的两个值，分别令其为x_1,x_2，它们互为相反数；第二次再次递归获得在x_1=\pm x_{i1},x_2=\pm x_{i2}四个值，这是从上到下的解释。

或者我们也可以从下到上的解释，假设我们需要取在四个点x=x_1,-x_1,x_2,-x_2处的值，那么我们只需要计算P_e[x_1^2],P_o[x_1^2]以及P_e[x_2^2],P_o[x_2^2]，因为互为相反数的两个数平方后是相等的，所以原来的4次遍历取值运算现在变为了2次对P_e[x^2],P_o[x^2]的取值运算。此时我们再次考虑x_1^2,x_2^2也是一对相反数，令其为x_i，那么我们最后只需要计算P_e[x_i^2],P_o[x_i^2]，也就是1次运算，我们的运算量在以指数级别的速度下降！

但是，我们需要讨论一个问题
::: align-center
具体取哪些值？或者说，满足上面要求的x_i存在吗？

:::

因为要实现上面的这个算法，我们需要一个x_i，它开2次根后的两个解是互为相反数的，而每个解再开二次根后的两个解依然是相反数，如果我们能找到这样一个x_i，那么上面的减少P_e[x^2],P_o[x^2]运算的算法就是成立的，对于上面的问题，我们只需要得到x_i的两个二次方根解和四个四次方根解，之后用两个二次方根解计算P_e[x^2],P_o[x^2]，然后使用四个四次方根解和前面的结果，代入

::: align-center
\begin{aligned} &P(x_i)=P_e[x_i^2]+\mathbf{x_i}P_o(x_i^2)\\ &P(-x_i)=P_e[x_i^2]-\mathbf{x_i}P_o(x_i^2) \end{aligned}

:::
下面的两个式子，最后就可以实现相较于原算法而言指数级别的规模下降，我们的算法复杂度将下降到\mathcal{O}(n\log n)

所以问题来了
::: align-center
这样的一个“魔数”存在吗？

:::

“神奇”的1的n次方根

首先，我们需要明白，在前一节中提到的“不停开二次根”的操作在实数范围内已经不成立了，因为对于一个实数而言它的二次方根在实数范围内可能是没有解的。所以要解决这个问题，首先需要做的就是将范围扩充到复数域。此时若我们考虑1的n次方根，它在复数域下的性质，完美地满足我们的要求。

首先，1的任意偶次方根的解一定是两两互为相反数的，考虑其n次方根在欧拉公式下的表示
::: align-center
w_N^j=e^{\frac{2\pi i}{n} j},j=0,1,2\cdots,n-1

:::
容易证明
::: align-center
w_N^{j+\frac{n}{2}}=e^{\pi i}w_N^j=-w_N^j

:::
也就是说，将1的n次方根在n为偶数时一定是两两相反配对的，而且，将这些n次方根平方以后，它们依然两两配对。我们的多项式取值的问题，居然隐藏在“1”之中

库利-图基算法的引入

到了这里我们不再讨论如何计算多项式乘积的问题，我们只需要多项式“取值”分治策略的思想

我们知道，在DFT中（请见这篇文章）我们需要利用到w^j进行变换计算，这里我们不妨考虑，所谓的DFT，其实就是一个多项式乘积的问题，我们的离散信号x_j,j=1,2,\cdots,n=1就是所谓的“多项式的系数”（对完整信号的离散采样），它们代表着一个唯一确定的“多项式”（信号），对于某一个给定的数k，现在要对x_jw_N^{k\cdot 0},x_jw_N^{k\cdot 1},\cdots,x_jw_N^{k\cdot(n-1)}依次进行“取值”操作，相当于有j=1,2,\cdots,n-1共n个信号，所以取1的“n次方根”（这里的想法已经很接近笔者在DFT文章中的想法了）

需要注意的是常数k其实不会影响结果，因为w_N^{k\cdot(n-1)}由于复数的性质，等同于w_N^{k\cdot(n-1)\space mod \space N}，即与此项等价的复数幂次永远不会超过N-1，我们最多只会有N-1个不同的w_N^{k\cdot(n-1)}

考虑将信号索引j类比为前面“多项式的幂次”，将w_N^{k\cdot(n-1)}按照j的奇偶分为两个部分,此时每个部分均有\frac{n}{2}个信号，对应的，我们取1的\frac{n}{2}次方根

这里我们依然类比假设信号有4个采样点，那么我们按奇偶，将其分为x_0,x_2和x_1,x_3，此时变换需要用到w_4^j，即1的4次方根，也即“取四个值”

所以对于我们的“多项式”来说，需要
::: align-center
\begin{aligned} &P(x_i)=P_e[x_i^2]+x_iP_o[x_i^2]\\ &P(-x_i)=P_e[x_i^2]-x_iP_o[x_i^2] \end{aligned}

:::
其中“取值点”x_i就是w_4^j
故
::: align-center
\begin{aligned} &P(w_4^j)=P_e[w_4^{2j}]+x_iP_o[w_4^{2j}]\\ &P(-w_4^j)=P_e[w_4^{2j}]-x_iP_o[w_4^{2j}] \end{aligned}

:::
进一步利用上面的性质
::: align-center
\begin{aligned} &P(w_4^j)=P_e[w_4^{2j}]+w_4^jP_o[w_4^{2j}]\\ &P(w_4^{j+\frac{n}{2}})=P_e[w_4^{2j}]-w_4^jP_o[w_4^{2j}] \end{aligned}

:::
这里的P已经不再是函数的对应关系了，而是FFT的对应关系
进一步
::: align-center
\begin{aligned} &FFT[w_4^j]=FFT_e[w_4^{2j}]+w_4^jFFT_o[w_4^{2j}]\\ &FFT[w_4^{j+\frac{n}{2}}]=FFT_e[w_4^{2j}]-w_4^jFFT_o[w_4^{2j}] \end{aligned}

:::
其中FFT是伪算子，其后跟的是此变换需要对信号施加的旋转因子w_N^j

我们将信号不断拆分为奇偶部分，对其进行FFT，之后使用上面的对称性性质的式子（又乘蝶形算法公式）对结果进行合并，如此我们的计算的规模会指数级的下降

而逆FFT的做法，也是类似的，只不过旋转因子变成了w_n^{-j}，而且加上了\frac{1}{n}修正了频域到时域的能量

这就是Cooley-Tukey FFT（库利-图基FFT）算法的精髓所在，它利用了DFT中1的n次方根的对称性，采用分治策略，不断缩减问题的规模，使得复杂度从DFT的\mathcal{O}(n^2)缩减到了\mathcal{O}(n\log n)

库利-图基算法的正式介绍

另一种解释

这里我们希望使用另一种更为“数学”的解释

我们把N次DFT分解为N_1N_2次DFT，即N=N_1\times N_2次DFT，读者可以将其理解为N_1步每次处理N_2点的小DFT,亦或者是N_2步每次处理N_1点的小DFT

我们的目标是寻找这诸如N_2步里面的每个N_1点小DFT之间存不存在一些“关系”，使我们可以节省运算。在数学上，这相当于我们在寻找这种“新的表示方法”下的DFT计算式中有无诸如相互抵消的项之类的现象。

所以不妨先将要DFT的信号分解为N_2个N_1点信号
::: align-center
\begin{aligned} &x_0,x_1\cdots x_{N_1-1}\\ &x_{N_1},x_{N_1+1},\cdots x_{2N_1-1}\\ &\vdots\\ &x_{(N_2-1)N_1},x_{(N_2-1)N_1+1},\cdots,x_{N_2N_1-1} \end{aligned}

:::

在这种“按行分块”的表示法下，信号的“列数”是固定的，我们通过调整行号a(0\le a\le N_2-1)和由每行的第b(0\le b\le N_1-1)个来索引信号点，如此信号索引j的新表示方法为
::: align-center
j(a,b)=aN_1+b,0\le a\le N_2-1,0\le b\le N_1-1

:::

现在我们来考虑旋转因子w_N，对于每个N_1点的小DFT，每个小DFT中的第k个小信号要施加的旋转因子w_N是相同的！

这也就意味着，对于上面的信号排布来说，每一列的信号点要施加的旋转因子是相同的，这也就意味着这些相同列上的信号经DFT变换后在一个频谱上，所以我们可以在旋转因子的索引上“按列分块”，类似的可得DFT频谱索引k的新表示方法为
::: align-center
k(c,d)=cN_2+d,0\le c\le N_1-1,0\le d\le N_2-1

:::

索引编制完成后，我们开始正式探究DFT式子中有无简化计算的方法，将原DFT式子中的j,k修改为新索引
::: align-center
\begin{aligned} &\hat{x}_k=\sum_{j=0}^{N-1}x_jw_N^{jk}\\ &=\sum_{a,b}x_{a,b}w_N^{(aN_1+b)(cN_2+d)}\\ &=\sum_{b=0}^{N_1-1}w_N^{b(cN_2+d)}\sum_{a=0}^{N_2-1}x_{a,b}w_N^{aN_1(cN_2+d)} \end{aligned}

:::
这里我们展开w_N^{aN_1(cN_2+d)}
::: align-center
\begin{aligned} &w_N^{aN_1(cN_2+d)}=w_N^{acN_1N_2}w_N^{adN_1}\\ &=w_{N_2}^{ad} \end{aligned}

:::
前一项为1，因为N=N_1N_2，取模N简化计算，得w_N^0=1
后一项展开后可以由\frac{N_1}{N}=\frac{1}{N_2}化简

如此式子变为
::: align-center
\begin{aligned} \sum_{b=0}^{N_1-1}w_N^{b(cN_2+d)}\sum_{a=0}^{N_2-1}x_{a,b}w_{N_2}^{ad} \end{aligned}

:::

此时很容易知道我们的式子消去了一些项，因此运算复杂度是下降的，分析其复杂度

• 内部求和，考虑变量a,b,d，复杂度\mathcal{O}(N_2^2N_1)
• 外部求和，考虑变量b,c,d，复杂度\mathcal{O}(N_1^2N_2)
• 共计复杂度\mathcal{O}(N_1N_2(N_1+N_2))

通常采用二基底FFT，即将信号分为两部分（通常是奇偶）N_1=\frac{N}{2},N_2=2，此时单次FFT复杂度已经已经降为\mathcal{O}(\frac{N^2}{2})，即原DFT算法的一半复杂度，若继续递归二分，最终复杂度将收敛于\mathcal{O}(n\log n)

表示

工程：蝶形算法

可以使用蝶形算法方便的表示，所谓的“蝶形”值得就是刚才我们所说的蝶形运算：
::: align-center
\begin{aligned} &X_1(k)=x_1(k)+w_N^{k}x_2(k)\\ &X_2(k)=x_1(k)-w_N^kx_2(k) \end{aligned}

:::
它可以使用下面的图表示：

这里读者可以方便的看出两个运算使用了同一个旋转因子w_N^k

我们以8点基2FFT为例，根据前面的推导，单次N点的DFT经由FFT可以表示为对2个\frac{N}{2}小DFT（内层求和）的蝶形运算（外层求和）

进一步分解两个\frac{N}{2}DFT为FFT的蝶形运算

进一步分解下去，整个FFT就可以只使用蝶形运算图表示，最后得到下面8点基2FFT的蝶形运算图：

蝶形运算的表示方法使得FFT在工程设计上变得非常直观

计算机科学：伪代码

# P为信号点集
function FFT(P):
  n=lenth(P)
  if (n==1):
    return P
  w = exp(-2 * pi * i / n)
  # 奇数索引部分
  p_even = [p0,p2,...,pn-2]
  # 偶数索引部分
  p_odd = [p1,p3,...,pn-1]
  # 递归
  fft_p_even = FFT(p_even)
  fft_p_odd = FFT(p_odd)
  # 初始化最后的结果
  y = [0] * n
  #蝶形运算
  for j in range(n/2):
    y[j] = fft_p_even[j] + (w**j) * fft_p_odd[j]
    y[j+n/2] = fft_p_even[j] - (w**j) * fft_p_odd[j]

  return y

应用-多项式乘法

FFT可用于加速多项式乘法，这主要利用了圆周卷积的性质

多项式的处理

首先，多项式乘法是一个线性卷积运算
考虑下面的多项式乘积，其中A(x)为n次多项式，B(x)为m次多项式
::: align-center
\begin{aligned} &C(x)=A(x)B(x)\\ &=(\sum_{j=0}^n a_jx^j)(\sum_{k=0}^mb_kx^k)\\ &=\sum_{j=0}^n\sum_{k=0}^ma_jb_kx^{j+k}\\ &=\sum_{l=0}^{n+m}(\sum_{l=j+k}a_jb_k)\\ \end{aligned}

:::
接下来，考虑内层求和的上下限，已知多项式A(x)系数索引的范围为0\le j\le n，B(x)系数索引k的范围为0\le k\le m，令新索引r=j\Rightarrow k=l-r，考虑

• 考虑索引r=j，它应在[0,n]内
• 考虑索引k，它在[0,m]内，所以0\le l-r\le m\Rightarrow l-m\le r\le l
• 综合以上结果，上限为\min(l,n)，下限为\max(0,l-m)

所以内层的系数计算，实际上为一个线性卷积
::: align-center
\sum_{r=\max(0,l-m)}^{\min(l,n)}a_rb_{l-r}

:::

要使用FFT加速这个线性卷积，首先要考虑将线性卷积表示为圆周卷积

对于多项式的系数向量x,y，两个多项式乘积多项式的系数
若两个向量等长，直接计算圆周卷积，
::: align-center
x\circledast y=\sum_{l=0}^{N-1}x_ly_{j-l\space mod \space M}

:::
此时考虑y的索引j-l，当j-l<0，即j时，圆周卷积的取模运算会使索引映射到N+j-l，这会导致原来线性卷积的后半部分叠加到前半部分

为了正确的使得圆周卷积的结果和线性卷积的结果相同，需要将两个多项式的系数（x为N维，y为M维）向量补0增广到N+M-1维以上
::: align-center
\begin{aligned} &x=(x_0,x_1,\cdots,x_{N-1},0,\cdots,0)\\ &y=(y_0,y_1,\cdots,y_{M-1},0,\cdots,0) \end{aligned}

:::

当系数长度扩充到L\ge N+M-1时，对任意卷积后的系数索引0\le j\le N+M-2，有

• x的索引l\le N-1时：依然有效
• y的索引(j-l)\le M-1时：依然有效
  实际上，可以由不等式的性质得到j-l的范围为-(N-1)\le j-l\le N+M-2，在圆周卷积中，索引溢出的部分-(N-1)\le j-l<0取模L会得到L+(j-l)
  但因为L+(j-l)\ge L-(N-1)\ge(N+M-1)-(N-1)=M
  又因为y的M到L的部分都是补充的0，所以不会出现原来直接计算圆周卷积时溢出索引取模后与原来的部分重叠的问题
• 其他情况：均为0，对结果无贡献

一般地来说，设多项式系数向量x,y的维数分别为n,m维，在这一步我们需要将两个系数向量补零增广到n+m-1维

FFT加速

将多项式乘法表示为圆周卷积后，因为DFT有卷积定理（对于FFT也成立）：时域圆周卷积等于频域乘积
::: align-center
\mathcal{F}(x\circledast y)=\mathcal{F}(x)\times\mathcal{F}(y)

:::
所以我们可以将两个多项式系数分别作FFT后直接逐点相乘，之后对结果进行逆FFT变换，取前N+M-1个点，就是多项式系数乘积的结果

在原来计算线性卷积时，我们需要遍历j,l，复杂度为\mathcal{O}(n^2)，而使用FFT+逐点乘积+IFFT加速后我们的复杂度只有\mathcal{O}(n\log n)，相比原来有质的提升