声明

首先，我需要声明，本文是在转载的基础上稍微修饰的，经过原创作者 ZLH_HHHH（佐理慧学姐）的许可方才转载并修饰的，由于我就是初学者，并且是数学渣滓，所以我学姐建议我写一下残疾人手册，我当然是欣然接受！！！

正文：

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

直接正题吧：

首先： DFT 本身并不负责多项式之间的乘法。

DFT 只是一种变换。

FFT 则是DFT的快速算法。（分治提高效率）

利用 FFT 。我们快速的将多项式变换为利于计算的形式

用这种方便计算的形式计算出来两个多项式的乘积。

这时候我们虽然已经得到目标多项式。但其形式并不是我们想要的

所以之后利用 FFT 的逆运算又快速的变换回去

我们记： FFT 的逆运算为 FFT− 1

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

$A (x) = \sum i = 0 n - 1 a i x i$

$这里的 n 次多项式是指最高次项指数为 n - 1 的多项式$ 2

小学生手算系数形式的多项式乘法的复杂度是 O(n2)

如果我们知道了一个 n 次多项式的曲线上的 n 个不同的的点。

我们是可以计算出来这个多项式的

why？

把系数看做未知数。列出来n个方程组。解出来这n个系数。

也就是说，给出曲线上的n个点：

$< (x 0, y 0), (x 1, y 1), (x 2, y 2), . . . . . (x n - 1, y n - 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) = \sum i = 0 n - 1 b i x i$

设 A(x) 与 B(x) 在取相同x的点值表达为：

$< (x 0, A 0), (x 1, A 1), (x 2, A 2), . . . . . (x 2 n - 2, A 2 n - 2) > < (x 0, B 0), (x 1, B 1), (x 2, B 2), . . . . . (x 2 n - 2, B 2 n - 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) 为：

$< (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) >$ <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) = a 0 + a 2 x + a 4 x 2 + . . . + a n - 2 x n 2 + 1$

设其系数的下标为奇数组成的多项式为 A[1](x) ：

$A [1] (x) = a 1 + a 3 x + a 5 x 2 + . . . + a n - 1 x n 2 + 1$

所以有：

$A (x i) = A [0] (x 2 i) + x i A [1] (x 2 i)$ 4

好像并不会优化运算。

因为 $< x 20, x 21, . . . . x 2 n - 1 >$ <script id="MathJax-Element-41" type="math/tex; mode=display"> </script>

依然有可能会组成n个不同的数字。

那么我们要计算的规模不会减半。

如果我们恰当选取一些x。是否会优化运算呢。

例如：

$< x 0, x 1, . . . . x n 2 - 1, - x 0, - x 1, . . . . - x n 2 - 1 >$ <script id="MathJax-Element-42" type="math/tex; mode=display"> </script>

各个数字平方后。得到的集合大小减半：

$< x 20, x 21, . . . . x 2 n 2 - 1 >$ <script id="MathJax-Element-43" type="math/tex; mode=display"> </script>

因为 (−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)$

但是这种关系不会传递下去。平方后得到的数字全是不小于0的数字。

再次递归又回到了原来的形式。很尴尬。

不过这启迪我们。取相反数。然后平放。可以把问题规模减半。但再一次递归就失效了。是否存在不失效的取值呢。

<script id="MathJax-Element-46" type="math/tex; mode=display"></script>

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

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

分治就是可行的。

有一种东西。可以满足我们的需求。

那就是 —— n次单位复数根：

（使用复数确实让人有顾虑。这也是我学习FFT时最大思想负担）

我们定义 i=−1‾‾‾√

那么形如 $c o s θ + i * s i n θ$

的复数有诸多良好的性质。

例如： $(c o s θ + i * s i n θ) k = c o s k θ + i * s i n k θ$

上面的性质可以用数学归纳法得到（较为简单。这里不做证明

我们记：

$W n = c o s 2 π n + i * s i n 2 π n$

我会告诉你多项式在x取下面的值时。有助于我们进行DFT

$< W 0 n, W 1 n, . . . . . W n - 1 n >$ <script id="MathJax-Element-54" type="math/tex; mode=display"> </script>

上面的n个复数称为。n次单位复数根。（n次单位复数根是指这n个数）

如果我们x取 Wkn 时。根据刚才的分解：

$A (W k n) = A [0] ((W k n) 2) + W k n A [1] ((W k n) 2)$

而 $(W k n) 2 = W 2 k n = 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 k n 2$

因为三角函数的周期性：

我们有：

$W k n = 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 k m o d n n$

所以： $(W k n) 2 = W k n 2 = W k m o d n 2 n 2$

则： $< (W 0 n) 2, (W 1 n) 2, . . . (W n - 1 n) 2 > 等效于 < W 0 n 2, W 1 n 2, . . . W n 2 - 1 n 2, W 0 n 2, W 1 n 2, . . . W n 2 - 1 n 2 >$ <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 的取值集合取单位复数根。不但满足折半的性质。而且还有一定的规律性。与原集合保持一致。

这意味着我们只需要计算取：

$< W 0 n 2, W 1 n 2, . . . W n 2 - 1 n 2 >$

的 A[0],A[1]

问题范围减半。并且如果n为偶数。再次平方依然集合大小减半。

记： $A [0] k = A [0] (W k n 2) = A [0] ((W k n) 2) A [1] k = A [1] (W k n 2) = A [1] ((W k n) 2)$

那么，当我们得到：

$< A [0] 0, A [0] 1, . . A [0] n 2 - 1, A [1] 0, A [1] 1, . . . A [1] n 2 - 1 >$ <script id="MathJax-Element-65" type="math/tex; mode=display"> </script>

这意味着我们得到了所有的 A[0]((Win)2),A[1]((Win)2)

则：

$A (W k n) = A [0] k m o d n 2 + W k n A [1] k m o d n 2$

即得到了： $< A 0, A 1, A 2, . . . A n - 1 >$ <script id="MathJax-Element-68" type="math/tex; mode=display"> </script>

其中 Ak=A(Wkn)

当然。FFT的计算 n要取2的整数次幂。这是因为每次减半。我们可以把不足的系数用0填充。

上面过程可以看作，取单位复数根计算目标y向量，如下图：

$⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ (W 0 n) 0 (W 1 n) 0 ⋮ (W n - 1 n) 0 (W 0 n) 1 (W 1 n) 1 ⋮ (W n - 1 n) 1 \dots \dots ⋱ \dots (W 0 n) n - 1 (W 1 n) n - 1 ⋮ (W n - 1 n) n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ a 0 a 1 ⋮ a n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ A 0 A 1 ⋮ A n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥$

本着尽可能简单的原则。我们不在特别的说这个矩阵。

因为之前说关于系数的n个方程必然可解。所以上面的矩阵：

$⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ (W 0 n) 0 (W 1 n) 0 ⋮ (W n - 1 n) 0 (W 0 n) 1 (W 1 n) 1 ⋮ (W n - 1 n) 1 \dots \dots ⋱ \dots (W 0 n) n - 1 (W 1 n) n - 1 ⋮ (W n - 1 n) n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥$

必然存在逆矩阵5。（有线性代数基础的，或许不陌生）

可能你还是有点迷惑。没关系：

如果我们把FFT看作计算上面特别的矩阵相乘的算法呢。

利用FFT我们可以快速得到：

$⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ (W 0 n) 0 (W 1 n) 0 ⋮ (W n - 1 n) 0 (W 0 n) 1 (W 1 n) 1 ⋮ (W n - 1 n) 1 \dots \dots ⋱ \dots (W 0 n) n - 1 (W 1 n) n - 1 ⋮ (W n - 1 n) n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ a 0 a 1 ⋮ a n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ - (F F T) - > ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ A 0 A 1 ⋮ A n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥$

同时。我可以直接告诉你： $D = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ (W 0 n) 0 (W 1 n) 0 ⋮ (W n - 1 n) 0 (W 0 n) 1 (W 1 n) 1 ⋮ (W n - 1 n) 1 \dots \dots ⋱ \dots (W 0 n) n - 1 (W 1 n) n - 1 ⋮ (W n - 1 n) n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥$

的逆矩阵为：

$E = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 1 n (W - 0 n) 0 1 n (W - 1 n) 0 ⋮ 1 n (W - (n - 1) n) 0 1 n (W - 0 n) 1 1 n (W - 1 n) 1 ⋮ 1 n (W - (n - 1) n) 1 \dots \dots ⋱ \dots 1 n (W - 0 n) n - 1 1 n (W - 1 n) n - 1 ⋮ 1 n (W - (n - 1) n) n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥$

因为我直接给你了答案。。。所以怀疑的话。我们可以计算一下看看。

我们记Y=E*D

$Y [t] [j] = \sum k = 0 n - 1 E [t] [k] * D [k] [j], 其中 t, j \in [0, n - 1]$

$Y [t] [j] = \sum k = 0 n - 1 1 n W - t k n * W k j n = 1 n \sum k = 0 n - 1 (W j - t n) k$

上面可以看作等比数列求和。

当 t==j 时： $Y [t] [j] = 1 n \sum k = 0 n - 1 1 k = 1$

当 t≠j 时,令 u=j−t ：

$Y [t] [j] = 1 n \sum k = 0 n - 1 (W u n) k = 1 n * ( W u n ) n - 1 W u n - 1 （等比数列公式）$

因为： $(W u n) n = W u n n = W u n m o d n n = W 0 n = 1$

所以, t≠j 时： $Y [t] [j] = 0$

所以。Y为单位矩阵。E,D互为逆矩阵得证。

所以：

$⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ (W 0 n) 0 (W - 1 n) 0 ⋮ (W - (n - 1) n) 0 (W 0 n) 1 (W - 1 n) 1 ⋮ (W - (n - 1) n) 1 \dots \dots ⋱ \dots (W 0 n) n - 1 (W - 1 n) n - 1 ⋮ (W - (n - 1) n) n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ A 0 A 1 ⋮ A n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ - (F F T -) - > ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ n a 0 n a 1 ⋮ n a n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥$

所以。 FFT− 其实就是 FFT 的 Wkn 变为 W−kn 得到。

并且所得目标向量每个元素除掉n就可以了6

FFT的高效实现：

通常。我们希望FFT的实现尽可能快。所以FFT算法也都使用迭代结构而不是递归结构。

下面介绍一种常用的去递归方法。

首先。上面介绍的FFT是只能处理2的整数次幂的次数界的多项式

所以不特别说明。都有 n=2k （系数不足的用0补足）

对于递归的第一步：

输入数组 $< a 0, a 1, . . . a n - 1 >$ <script id="MathJax-Element-91" type="math/tex; mode=display"> </script>

被分解为，偶数和奇数两个部分数

$< a 0, a 2, a 4, . . . a n - 2 >, < a 1, a 3, a 5, . . . a n - 1 >$ <script id="MathJax-Element-92" type="math/tex; mode=display"> \ \ \ ,\ \ \ </script>

下标二进制形式右起第一个字符为 0 被分配到左集合

下标二进制形式右起第一个字符为 1 被分配到右集合

更一般的。对于第k次递归：

下标右起第k个字符为 0 被分配到它所在问题的左集合

下标右起第k个字符为 1 被分配到它所在问题的右集合

那么对于一个2进制形式的数字 B。将其表示为(注意： n=2k )：

$B = \sum t = 0 k - 1 b t * 2 t$

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

$P (B) = \sum t = 0 k - 1 b t * 2 k - 1 - t, (因为每次都删去一半位置)$

所以。对于一个数字的二进制形式的前k位对称过来。就是递归后的位置。

例如： k=4 , (0011)2−对称后−>(1100)2

我们调用FFT前对数组进行一次上述重拍。

便可自底向上迭代实现FFT

总结：

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

不难写出非递归结构

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

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

方便起见。我们用符号： Co(a,b) 来表示复数 a+bi

那么多项式：

$A (x) = \sum t = 0 n - 1 a t x t = \sum t = 0 n - 1 C o (a t, 0) x t$

$W n = C o (c o s 2 π n, s i n 2 π n)$

$W k n = C o (c o s 2 k π n, s i n 2 k π n)$

令: $Y [t] = C o (a t, 0), t \in [0, n - 1]$

当n=1时：

$F F T (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] + W t n Y [1] [t] Y [t + n 2] = Y [0] [t] - W t n Y [1] [t] （其中 t \in [0, n 2 - 1] ）$

返回 Y[];

————————————————————————

我对上述思想解释一下

A在 <W0n,W1n,...Wn−1n> <script id="MathJax-Element-111" type="math/tex"> </script>处的值保存在Y中。

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

$A [0] (W t n 2) = Y [0] [t] A [1] (W t n 2) = 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;
}

残疾人手册

《FFT》模版
 《多项式与快速傅立叶变换》原文

命名部分不用太纠结，另一种比较通用的称谓是将 FFT 的两个过程正逆分别称为 DFT 与 IDFT 。 ↩
这里沿用《算法导论》左闭右开的原则，所以重新定义一下。 ↩
C(x)=A(x)∗B(x)=AiBi ↩
因为在拆分的时候降次了，所以我们需要代入的是 x2i ，并且在奇数的情况前加上 xi 。 ↩
逆矩阵是 FFT− 的核心，如果不存在逆矩阵将无法逆转运算。 ↩
在逆运算过程中因为为了方便运算两边都乘以了 n ，所以最后的结果是我们所需要的结果的 n 倍，除以即可。 ↩

FFT（最详细最通俗的入门手册）

声明

正文：