「数学」FFT 基础学习笔记

快速傅里叶变换 - FFT 算法

多项式

多项式是由称为未知数的变量和称为系数的常数通过有限次加减法、乘法以及自然数幂次的乘方运算得到的代数表达式。多项式中单项式个数为多项式的项数。因为多项式由单项式组成,规定最高次项的次数叫做多项式的次数

为方便表示,我们记多项式 a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1} A(x) ,即

\begin{equation}A(x)=\sum\limits_{i=0}^{n-1}a_ix^i\end{equation}

我们规定:严格大于 n 的任意整数为这个多项式的次数界,即次数的范围。多项式中任意单项式次数均严格小于次数界。

(下面无特殊说明,多项式均用 A(x) 等类似形式表示

多项式表示法

系数表示法

通过系数来描述多项式叫做 系数表示

\begin{equation}A(x)=\sum\limits_{i=0}^{n}a_ix^i\end{equation}

点值表示

我们把多项式 A(x) 看做一个函数,选取不同的 n 个值 x_0,x_1,\cdots,x_{n-1} 代入 A(x) 中,可以得到 A(x) 的图像上 n 个不同的点,那么称点集 \alpha=\{ (x_i,A(x_i))\mid i\in[0,n-1],i\in \mathbb{Z}\} 为多项式 A(x) 点值表示

多项式加减法

两个多项式相加可以看作是对两组单项式的和进行重组与合并同类项。通过加法结合律,可以将同类项放在一起,合并之后就得到了两个多项式的和。下面是系数表示法的形式

\sum _ { i = 0 } ^ { n } a _ { i } X ^ { i } + \sum _ { i = 0 } ^ { n } b _ { i } X ^ { i } = \sum _ { i = 0 } ^ { n } \left( a _ { i } + b _ { i } \right) X ^ { i }

多项式乘法

计算两个多项式相乘时,首先使用乘法对加法的分配律将各项拆出,然后运用乘法结合律整合每一项,最后和加法一样整合同类项,就能得到乘积多项式。下面是系数表示法的形式

\begin{align}\left( \sum _ { i = 0 } ^ { n } a _ { i } X ^ { i } \right) \left( \sum _ { j = 0 } ^ { m } b _ { j } X ^ { j } \right) = \sum _ { k = 0 } ^ { m + n } \left( \sum _ { \mu + \nu = k } a _ { \mu } b _ { \nu } \right) X ^ { k }\end{align}

\sum _ { \mu + \nu = k } a _ { \mu } b _ { \nu } 也正是卷积

这样,朴素算法时间复杂度为 O(n^2) ,显然不是太好。

但是点值表达式的表示就很好了,直接乘起来就好了嘛!即:


A(x) 的点值表示为

\alpha=\{ (x_i,A(x_i))\mid i\in[0,n-1],i\in \mathbb{Z}\}

B(x) 的点值表示为

\beta=\{ (x_i,B(x_i))\mid i\in[0,n-1],i\in \mathbb{Z}\}

然后就能得到 C(x) 的点值表示为:

\gamma=\{ (x_i,A(x_i)B(x_i))\mid i\in[0,n-1],i\in \mathbb{Z}\}

但是这里只是两多项式的积。并不是真正的点值表示,不过也够用了

差值 - 点值转化

就要接近重点了同志们加油! 先定义两个运算:

  1. 定义一种运算,可以将系数表示转化为点值表示,记为 \text{DFT}(A(x))=\alpha

  2. 定义其逆运算,即从一个多项式的点值表示确定其系数也就是多项式插值。记为 \text{DFT}^{-1}(\alpha)=\text{IDFT}(\alpha)=A(x)

但是注意到 \text{DFT}(A(x)) 的复杂度是 O(n^2) 的(暴力计算),它的逆运算是 O(n^3) 的(高斯消元),现在就需要一个更高效的算法

复数

单位根

n 次单位根 n 次幂为 1 的复数。它们位于复平面的单位圆上,构成正 n 边形的顶点,其中一个顶点是 1

k 次单位根有 k 个,分别为:

e ^ { \frac { 2 \pi k i } { n } } \quad ( k = 0,1,2 , \cdots , n - 1 )

上面的式子是用欧拉公式证明的

欧拉公式

e^{i\theta}=\cos\theta+i\sin\theta

证明:

对于 e^{ix} ,有 e^{ix}\times e^{-ix}=e^0=1 e^{ix}\neq 0

f(x) = {\cos x+i\sin x\over e^{ix}}

求导,有

\begin{aligned} f ^ { \prime } ( x ) & = \frac { ( - \sin x + i \cos x ) \cdot e ^ { i x } - ( \cos x + i \sin x ) \cdot i \cdot e ^ { i x } } { \left( e ^ { i x } \right) ^ { 2 } } \\ & = \frac { - \sin x \cdot e ^ { i x } - i ^ { 2 } \sin x \cdot e ^ { i x } } { \left( e ^ { i x } \right) ^ { 2 } } \\ & = \frac { - \sin x \cdot e ^ { i x } + \sin x \cdot e ^ { i x } } { \left( e ^ { i x } \right) ^ { 2 } } \\ & = 0 \end{aligned}

证毕

几个引理

下面这些也可以通过复平面感性理解

消去引理

内容:当 d\not =0 时, \omega_{dn}^{dk}=\omega_n^k

证明:

\begin{equation}\begin{aligned}\omega_{dn}^{dk}&=e^{\frac{2dk\pi i}{dn}}\\&=e^{\frac{2k\pi i}{n}}=\omega_n^k\end{aligned}\end{equation}

折半引理

内容:如果 n>0 n 为偶数, (\omega_n^k)^2=\omega_{\frac{n}{2}}^k 证明:

\begin{equation}\begin{aligned}(\omega_{n}^{k})^2&=(e^{\frac{2k\pi i}{n}})^2\\&=e^{\frac{2k\pi i}{\frac{n}{2}}}=\omega_{\frac{n}{2}}^k\end{aligned}\end{equation}

更一般的结论:

\omega_n^{k+\frac{n}{2}}=\omega_n^{\frac{n}{2}}\omega_n^k=-\omega_n^k \omega_n^{\frac{n}{2}}=e^{\pi i}=-1

求和引理

内容:对于 \forall k\in(0,n) ,有

\sum\limits_{j=0}^{n-1}(\omega_n^k)^j=0

证明:

\begin{align}\sum _ { k = 0 } ^ { n - 1 } e ^ { \frac { 2 \pi k i } { n } } &= \frac { e ^ { \frac { 2 \pi k n i } { n } } - 1 } { e ^ { \frac { 2 \pi i } { n } } - 1 }\\& = \frac { 1 - 1 } { e ^ { \frac { 2 \pi i } { n } } - 1 } = 0\end{align}

离散傅里叶变换

DFT ,即 \text{Discrete Fourier Transform} ,离散傅里叶变换,就是我们之前定义的那个,现在有了复数,我们就可以在 O(n\log_2 n) 的时间内做完 DFT 运算。运用到了 \text{Cooley-Tukey} 方法

我们将多项式按照指数的奇偶分类,记原多项式为 A(x) ,那么构造两个多项式

\begin{equation}\begin{aligned}A_0(x)&=a_0+a_2x+a_4x^2+\cdots +a_{n-2}x^{\frac{n}{2}-1}\\&=\sum\limits_{i=0}^{\frac{n}{2}-1}a_{2i}x^i\\A_1(x)&=a_1+a_3x+a_5x^3+\cdots +a_{n-1}x^{\frac{n}{2}-1}\\&=\sum\limits_{i=0}^{\frac{n}{2}-1}a_{2i+1}x^i\\\end{aligned}\end{equation}

通过观察可知, A(x)=A_0(x^2)+xA_1(x^2) ,根据折半引理,我们考虑将 n 次单位根作为 x 代入计算。又因为折半引理,我们即可分治操作qwq

那么,对于 k<\frac{n}{2} ,有:

\begin{equation}\begin{aligned}A(\omega_n^k)&=A_0((\omega_n^k)^2)+\omega_n^kA_1((\omega_n^k)^2)\\&=A_0(\omega_{\frac{n}{2}}^k)+\omega_n^kA_1(\omega_{\frac{n}{2}}^k)\\&=\sum\limits_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^k\sum\limits_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}\\A(\omega_n^{k+\frac{n}{2}})&=A_0((\omega_n^k)^2)+\omega_n^{k+\frac{n}{2}}A_1((\omega_n^k)^2)\\&=A_0(\omega_{\frac{n}{2}}^k)-\omega_n^kA_1(\omega_{\frac{n}{2}}^k)\\&=\sum\limits_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}-\omega_n^k\sum\limits_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}\end{aligned}\end{equation}

这样,由奇偶性使得需要代入的值范围减半,递归做下去就行了qwq

我们来分析一下时间复杂度qwq

\begin{equation}T(n)=2T\left(\frac{n}{2}\right)+O(n)=O(n\log_2 n)\end{equation}

因为永远是奇偶合并,如果想要更方便地合并必然 n 2 的整数次幂,否则合并到一定时候左右不一定都是相同次数(可以考虑一个完全二叉树)。所以做的时候取 2 的整数次幂作为 n

由于这个算法是实质是将 DFT 分解为较小长度的多个 DFT ,因此它可以同任一种其他的 DFT 算法联合使用

逆离散傅里叶变换

IDFT ,即 \text{Inverse Discrete Fourier Transform} ,逆离散傅里叶变换。我们在 O(n\log_2 n) 的时间里搞定了 DFT ,可是 IDFT 还是个 O(n^3) 的,所以继续优化算法

我们说过,IDFT 相当于求系数:

\begin{equation}\begin{bmatrix}(\omega_n^0)^0& (\omega_n^0)^1&(\omega_n^0)^2 &\cdots &(\omega_n^0)^{n-1}\\(\omega_n^1)^0& (\omega_n^1)^1&(\omega_n^1)^2&\cdots &(\omega_n^1)^{n-1}\\\vdots&\vdots &\vdots &\ddots &\vdots \\(\omega_n^{n-1})^0 & (\omega_n^{n-1})^1& (\omega_n^{n-1})^2&\cdots & (\omega_n^{n-1})^{n-1}\end{bmatrix}\times\begin{bmatrix}a_0\\a_1\\ \vdots \\a_{n-1}\end{bmatrix}=\begin{bmatrix}A(\omega_n^0)\\A(\omega_n^1)\\A(\omega_n^2)\\\vdots \\A(\omega_n^{n-1})\\\end{bmatrix}\end{equation}

这同时也是做 DFT 时我们乘的矩阵。只不过这次 a_i 变成未知数了

记上面的系数矩阵为 D ,再考虑以下矩阵:

\begin{bmatrix}(\omega_n^0)^0& (\omega_n^0)^1&(\omega_n^0)^2 &\cdots &(\omega_n^0)^{n-1}\\(\omega_n^{-1})^0& (\omega_n^{-1})^1&(\omega_n^{-1})^2&\cdots &(\omega_n^{-1})^{n-1}\\\vdots&\vdots &\vdots &\ddots &\vdots \\(\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1& (\omega_n^{-(n-1)})^2&\cdots & (\omega_n^{-(n-1)})^{n-1}\end{bmatrix}

记此矩阵为 W ,则记 E=D\times W 。可知:

\begin{equation}\begin{aligned}E_{ij}&=\sum\limits_{k=0}^{n-1}D_{ik}W_{kj}\\&=\sum\limits_{k=0}^{n-1}\omega_n^{-ik}\omega_n^{kj}\\&=\sum\limits_{k=0}^{n-1}\omega_n^{k(j-i)}\end{aligned}\end{equation}

i=j 时:

\begin{equation}\begin{aligned}E_{ij}&=\sum\limits_{k=0}^{n-1}\omega_n^{k(j-i)}\\&=\sum\limits_{k=0}^{n-1}\omega_n^0=n\\\end{aligned}\end{equation}

i\not =j 时,有

\begin{equation}\begin{bmatrix}a_0\\a_1\\ \vdots \\a_{n-1}\end{bmatrix}=\begin{bmatrix}(\omega_n^0)^0& (\omega_n^0)^1&(\omega_n^0)^2 &\cdots &(\omega_n^0)^{n-1}\\(\omega_n^{-1})^0& (\omega_n^{-1})^1&(\omega_n^{-1})^2&\cdots &(\omega_n^{-1})^{n-1}\\\vdots&\vdots &\vdots &\ddots &\vdots \\(\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1& (\omega_n^{-(n-1)})^2&\cdots & (\omega_n^{-(n-1)})^{n-1}\end{bmatrix}\times\begin{bmatrix}A(\omega_n^0)\\A(\omega_n^1)\\A(\omega_n^2)\\\vdots \\A(\omega_n^{n-1})\\\end{bmatrix}\end{equation}

这和 DFT 的操作很像,所以用 DFT 的过程实现 IDFT 即可

若有负号,就把 \omega_n^1 改成 \omega_n^{-1} 即可;

现在 IDFT 也是 O(n\log_2 n) 的了

至此,所有的时间都是 O(n\log_2 n) 的了!即 FFT 的复杂度为 O(n\log_2 n)

事实上,可以通过迭代的方式降低常数

比如对于系数为 \{a_{0\to 7}\} 按照奇偶性分类的时候,有

1
2
3
#0 : {0, 1, 2, 3, 4, 5, 6, 7}
#1 : {0, 2, 4, 6}, {1, 3, 5, 7}
#2 : {0, 4}, {2, 6}, {1, 5}, {3, 7}

对于 #1 ,就相当于对这些数的二进制的第一位进行分类。 #2 ,就相当于对这些数的二进制的第二位进行分类。通过简单的计算就可以得到,位置 i 最后的数就是 i 的二进制的逆序上最后的数

后记

快速傅立叶变换具有一些实现上的缺点,举例来说,原向量必须乘上复数系数的矩阵加以处理,而且每个复数系数的实部和虚部是一个正弦及余弦函数,因此大部分的系数都是浮点数,也就是说,我们必须做复数而且是浮点数的运算,因此计算量会比较大,而且浮点数运算产生的误差会比较大

奇偶分组这个玩意就很坑,容易发现,设 r_i 为将 i 的高位翻转到低位得到的数字(如 0011 就变成 1100 ),那么初始位置为 i 的数就会被移动到 r_i ,所以我们可以预先完成移动,然后合并计算。时间就会加快很多了

不过 FFT 本身常数巨大,主要因为 double 类型的计算时间,还有使用 \sin , \cos 函数带来的计算时间问题,所以常数问题一定要注意

我们看到, DFT 做完后,得到的是将 n 次单位根代入后表达式的点值表示, IDFT​ 做完后,得到的是两多项式相乘的系数表示

由于 FFT 涉及复数和除法,容易出现精度问题,所以考虑一种整数在模的前提下的 FFT ,就出现了 FNT ,也叫 NTT ,即快速数论变换Fast Number-theory TransformNumber Theory Transform)...

如果我们计算卷积的时候 a , b 的下标不是相加为定值,而是一或二元逻辑运算(与,或,异或等等)后为定值,则可应用快速沃尔什变换(Fast Walsh Hadamard Transform).....

FFT 真是大毒瘤……

FFT FWT 都在图像处理,音频处理方面有较大应用

下面给出代码

代码

(没有迭代的算法)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
// Code by Menci
const double PI = acos(-1);

bool inversed = false;

inline std::complex<double> omega(const int n, const int k) {
if (!inversed) return std::complex<double>(cos(2 * PI / n * k), sin(2 * PI / n * k));
else return std::conj(std::complex<double>(cos(2 * PI / n * k), sin(2 * PI / n * k)));
}

inline void transform(std::complex<double> *a, const int n) {
if (n == 1) return;

static std::complex<double> buf[MAXN];
const int m = n / 2;
// 按照系数奇偶划分为两半
for (int i = 0; i < m; i++) {
buf[i] = a[i * 2];
buf[i + m] = a[i * 2 + 1];
}
std::copy(buf, buf + n, a);

// 分治
std::complex<double> *a1 = a, *a2 = a + m;
fft(a1, m);
fft(a2, m);

// 合并两个子问题
for (int i = 0; i < m; i++) {
std::complex<double> x = omega(n, i);
buf[i] = a1[i] + x * a2[i];
buf[i + m] = a1[i] - x * a2[i];
}

std::copy(buf, buf + n, a);
}
坚持原创技术分享,您的支持将鼓励我继续创作!