线性积分与傅里叶变换

文章被题目大致分为了两个部分：线性积分和 傅里叶变换。其实这么分实际上是不甚妥当的，因为傅里叶变换就是线性积分变换中的一种。如果要结合题目细讲，写起来怕是一个浩大的工程了，因此讲解的内容大多就仅止步于在信息学奥赛中的应用。当然，博主也会尽力深挖，争取让文章不只是初步。

写在前面

很多人学快速傅里叶变换，学它的微分关系，学它的卷积特性，但是始终都没有思考过这个玩意是干嘛用的。毕竟本篇文章实际上想要讲的就是傅里叶变换，我们也首先要了解到傅里叶变换的作用是什么。

从数学上来说，傅里叶变换是用于解决两个多项式的 卷积，简单来说就是两个多项式相乘的次数，如果直接暴力计算，那么时间复杂度应该是$O(N^2)$，而快速傅里叶变换可以将时间复杂度降为$O(NlogN)$。这在时间维度上就体现了它宝贵的价值，因为多项式相乘的普遍性，复杂的快速傅里叶变换也渐渐走进人们的视野。

而从物理学或者工程学中，傅里叶变换的常见用途是 信号处理，而所属的信息与系统又是大学中大部分工科的基础，足以见其重要性。具体的讲就是 将给定信号把时间映射到振幅。而关于时间和频率的有关内容还请参考Heinrich关于傅里叶变换的教程：傅里叶分析之掐死教程 ，例如其在讯号处理的经典应用就是将讯号分解为振幅分量和频率分量。

线性积分变换

积分变换是数学中作用于函数的算子，用以处理微分方程等问题。

一个十分生涩的定义，我们假设当前有一个函数$f(x)$，将操作$T$表示为一个积分转换，而$f$经过此积分变换后的函数表示为$Tf(y)$：$Tf(y) = \int_{x_1}^{x_2}K(x, y)f(x) dx$。

其中$K()$是一个确定的含有两个参数的函数，称为此次积分变换的核函数（简称核）。而核使我们自己选择的，当我们选择不同的核函数$K()$或者积分域$(x_1,x_2)$就得到了不同的积分变换，这应该是很显然的。而相对应的，积分变换也有相对的反积分变换，也就是：$f(x) = \int_{y_1}^{y_2}K^{-1}(y, x)Tf(y)dy$

其中$K^{-1}(y, x)$被称为反核。

而积分变换与反积分变换构成了傅里叶变换的框架，在傅里叶变换中，核函数为$\frac{e^{ixy}}{\sqrt{2 \pi}}$，其积分域为$-\infty, \infty$。

下面我们会从基础开始一步步地将傅里叶变换拆成我们好理解的东西。

傅里叶变换

多项式

我想来到这里的人肯定不需要我再讲这东西了，如果果真有人不会，还请自行维基。为了讲解需要，在此强调几个多项式的概念性内容：

多项式的框架式形式：$\sum_{i = 0}^{N - 1}a_ix_i$

若其最高次的非零系数为$a_k$，则该多项式的次数为k。

任何一个严格大于其次数的整数都是这个多项式的次数界。

对于多项式的计算，加法显然很简单。$C(x) = A(x) + B(x) = \sum_{i = 0}^{N - 1}a_ix^i + \sum_{i = 0}^{N - 1}b_ix^i = \sum_{i = 1}^{N - 1}c_ix^i$

只需要将对应次数的系数相加即可，我们可以在$O(N)$的时间复杂度内得到答案。

而对于多项式乘法，方法是将$A(x)$中的每一项与$B(x)$中的每一项相乘，之后合并同类项即可。$C(x) = \sum_{i = 0}^{2N - 2}c_ix_i$其中$c_i = \sum_{j = 0}^ia_jb_i-j$称为两个多项式的卷积，表示为$a \otimes b$

对于一个多项式来说我们有两种表示方法：系数表示法和点值表示法，其实两种表示方法应该是等价的，但是在计算上来说，对于点值表示的多项式，求解的时间复杂度是$O(N)$，这已经非常优秀，但是系数表示的多项式的暴力求解却是$O(N^2)$，而快速傅里叶变换就可以做到将系数表示的多项式在$O(NlogN)$的时间内求解。

对于点值表示的多项式求值运算的逆运算（已知点值表示的多项式求系数表示）称为插值，关于插值内容还请参照博主的另一篇博客：对插值法及拉格朗日插值多项式的初步理解运用。

单位复数根

N次单位复数根$\omega$满足$\omega ^N = 1$

定义很好理解，但是这一部分作为前置知识的内容实际上算是难点了。下图表示了$n$个单位复数根均匀的分布在以复平面的原点为圆心的单位半径的圆上。

而这究竟是怎么理解呢？Heinrich的教程很形象的做了解释：

我们假设有一条长度为1的线段在一条数轴上，那么乘以-1之后的线段就与其相反，而$i = \sqrt{-1}$，那么我们将其乘以i的线段就在一个垂直的虚数轴上，于是我们就得到了一个由一个实数轴做x轴，一个虚数轴做y轴的复平面。

然后这里就可以引入另一个重要的公式：欧拉公式：$e^{ix} = cosx + i sinx$。运用这个式子我们就得到了上图的8次单位复数根，并且也知道n次单位复数根一共有n个。而这些单位根就是$e^{\frac{2 \pi i k}{n}}$我们把n次单位根的第m个写作$\omega_n^m$。

下面是单位复数根的几个性质或推论：

1. $\omega _{d \times n}^{d \times m} = \omega _n^m$，其中$d > 0$
2. 对于任何偶数$n>0$有$\omega _n^{n/2} = \omega_2 = -1$
3. 对于偶数$n>0$，n个n次单位复数根的平方的集合就是$\frac n2$个$\frac n2$次单位复数根的集合。
4. 对于任意$n \geq1$ 和不能被n整除的非负整数k，有$\sum_{i = 0}^{N - 1}(\omega_n^m)^i = 0$

至于具体的证明就不再说了，都是很简单的推论，博友们大可自行脑补，权当是练习了。

离散傅里叶变换

好了铺垫完所有的前置知识之后我们终于开始了傅里叶变换的内容。对于一个次数界为n的多项式$F(x) = \sum_{i = 1} ^{n - 1}a_ix^i$，我们在插值方面的知识可以知道n个点就可以确定这个多项式，那么我们同样也可以代入n次的n个单位复数根，来确定这个多项式，可以 求出其结果$k_k = F(\omega_n^k) = \sum_{i = 0}^{n - 1}a_i \omega _n^{k_i}$。

而$y​$就被称为是系数向量a的离散傅里叶变换。

这个东西计算的时间复杂度依然是$O(N^2)$，那它到底有什么用处呢？这在算法上或许无法体现，但是它是傅里叶变换在时域和频域上都呈现离散的形式，具体依然可以参考Heinrich的讲解，使其在频谱分析和数据压缩等领域发挥了巨大优势。

快速傅里叶变换

也就是常说的FFT，可以说是对于DFT（离散傅里叶变换）在时间复杂度上的一大改进。

设多项式$F(x)$的系数向量$a(a_0, a_1, a_2…a_{n - 1})$，对于多项式$F(x) = \sum_{i = 0} ^ {n - 1}a_ix_i$来说，可以按照下标的奇偶性分为两个部分：

$$F(x) = (a_0 + a_2x^2 +...+ a_{n - 2}x^{n - 2}) + (a_1 + a_3x^3...a_{n - 1}x^{n - 1})$$

我们设

$$F_1(x) = a_0 + a_2x + a_4x^+...+ a_{n - 2}x^{\frac n2 - 1}$$ $$F_2(x) = a_1 + a_3x + a_5x^+...+ a_{n - 1}x^{\frac n2 - 1}$$

则

$$F(x) = F_1(x^2) + xF_2(x^2)$$

和离散傅里叶变换一样，我们把$\omega _n^k$代入得到

$$F(\omega _n^k) = F_1(\omega _n^{2k}) + \omega _n^kF_2(\omega _n^{2k})$$

然后同理，将$\omega _n^{k + \frac n2}$代入可以得到：

$$F(\omega _n^{k + \frac n2}) = F_1(\omega _n^{2k + n}) - F(\omega _n^kF_2(F(\omega _n^{2k} \times \omega_n^n)$$

最后得到

$$F(F(\omega _n^{k + \frac n2}) = F_1(\omega _n^{2k}) - \omega _n^kF_2(\omega _n^{2k})$$

如此我们发现，两个式子化简后，只有一个常数项相反，其他都相同，于是再计算出第一个式子的时候我们可以$O(1)$求出第二个式子，并且发现两个式子将最初的范围缩小了一半，于是原问题缩小了一般，然后我们发现子问题也满足原问题性质，$k$与$k +\frac n2$同时取遍了(0, n - 1)，于是时间复杂度缩小到$O(NlogN)$。

这个算法在FFT的计算中最为常见，是在1965年由J.W.Cooley和J.W.Tuky提出的，因此也被称为Cooley-Tukey算法，实际上是基于分治的思想实现的，

快速傅里叶逆变换

也称为傅里叶反变换，上文我们在线性积分变换中提到过反积分变换，而傅里叶反变换就是傅里叶变换的反积分变换，在数学中的意义是点值表示的多项式转化为系数表示，也就是说我们要从点值向量$(a_0, a_1, a_2…a_{n - 1})$得到系数向量$(y_0, y_1, y_2…y_{n - 1})$。则我们可以设$(y_0, y_1, y_2…y_{n - 1})$是$(a_0, a_1, a_2…a_{n - 1})$进行傅里叶逆变换后的结果，并
设有多项式$F(x) = \sum_{i = 0} ^{n - 1}y_0x^i$，
假设有向量$(c_0, c_1, c_2…c_{n - 1})$表示多项式在$\omega_n^0, \omega_n^{-1}…\omega_n^{1 - n}$上的点值表示，则向量c满足

$$c_k = \sum_{i = 0} ^{n - 1}y_i(\omega_n^{-k})^i$$

然后我们考虑将式子展开。

$$c_k = \sum_{i = 0} ^ {n - 1}y_i(\omega _n^{-k}) ^ i$$ $$= \sum_{i = 0} ^ {n - 1}(\sum_{j = 0} ^ {n - 1}a_j(\omega _n^{j} k) ^ j)(\omega_n^{-k})^i$$ $$= \sum_{i = 0}^{n - 1}(\sum_{j = 0} ^{n - 1}a_j(\omega_n^j)^i)(\omega_n^{-k})^i ~=~ \sum_{i = 0}^{n - 1}(\sum_{j = 0} ^{n - 1}a_j(\omega_n^j)^i(\omega_n^{-k})^i)$$ $$= \sum_{i = 0}^{n - 1}\sum_{j = 0}^{n - 1}a_j(\omega_n^{j - k})^i ~=~ \sum_{j = 0}^{n - 1}a_j(\sum_{i = 0}^{n - 1}(\omega_n^{j - k})^i)$$

设一个式子

$$S(\omega_n^k) = \sum_{i = 0} ^{n - 1}(\omega_n^k)^i$$

在两边同时乘一个$\omega_n^k$得：

$$\omega_n^kS(\omega_n^k) = \sum_{i = 0} ^{n - 1}(\omega_n^k)^{(i+1)}$$

将前后两个式子相减，得到

$$\omega _n^kS(\omega_n^k) - S(\omega _n^k) = (\omega_n^k)^n - 1$$ $$S(\omega _n^k) = \frac {(\omega_n^k)^n - 1}{\omega_n^k - 1}= 0$$

则当$k = 0$的时候，其值为0

而原来的式子

$$c_k = \sum_{j = 0}^{n - 1}a_j(\sum_{i = 0} ^ {j - k})^i = \sum_{j = 0} ^ {n-1}a_jS(\omega_n^{j - k})$$

当$j = k$的时候$S = n$，否则$S = 0$，所以$a_i = \frac 1n c_i$

这样我们得到了点值和系数的关系式，所以，从结论来说，将单位根幂上-1，然后做一次快速傅里叶变换，将结果的数除以n就是傅里叶逆变换的结果。

代码实现

一般的分治都会采用递归操作实现，但是在FFT中递归实现的效率比较低，所以大多数情况下采用迭代操作实现。我们观察二进制下的序列规律

发现分治到边界条件的时候每个数的下标等于原下标的二进制翻转，于是我们省去了分奇偶的操作，直接枚举每一个二进制位，然后向上合并就可以了。下面是迭代操作实现的代码。

int l = 0 ; while ((1 << l) < N) l ++ ;
for (int i = 0 ; i < N ; i ++) {
	 int T = 0 ;
	 for (int j = 0 ; j < l ; j ++)
	 	if (i & (1 << j)) T /= (1 << (k - j - 1)) ;
	 if (i < t) swap(a[i], a[T]) ;
}

蝴蝶操作

我们把$\omega_n^k$称为旋转因子并保留下来，进入两个向量$a_k^0$和$a_k^1$时，旋转因子$\omega_n^k$乘以$a_k^1$，输出与$a_k^0$的和与差，这一个操作被称为蝴蝶操作。那这个东西究竟有什么作用呢？

重新考虑我们向上合并两个子问题时，假设有$A_1(\omega_{n/2}^k)$储存在$a(k)$中，$A_2(\omega _{n/2}^k)$储存在$a(k +\frac n2)$中，并且这两个值将要被储存在$b(k)$和$b(k + \frac n2)$中，则合并的操作可以如下表示：

$$b(k) ← a(k) + a(k +\frac n2) \times\omega_n^k$$ $$b(k +\frac n2) ← a(k) - a(k +\frac n2) \times\omega_n^k$$

于是我们考虑将两个值存放在原来的a中，取消b数组的存在，但是需要覆盖原来的值，所以就需要一个临时变量T。

$$T ←a(k +\frac n2) \times\omega_n^k$$ $$a(k +\frac n2) ← a(k) - T$$ $$a(k) ← a(k) - T$$

名字很好听，其实操作很简单，仅仅有一行而已。

模板代码

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std ;
typedef long long LL ;
const int MAXN = 10000010 ;

inline int Read() {
	int X = 0, F = 1 ; char ch = getchar() ;
	while (ch > '9' || ch < '0') F = (ch == '-' ? - 1 : 1), ch = getchar() ;
	while (ch >= '0' && ch <= '9') X=(X<<1)+(X<<3)+(ch^48), ch = getchar() ;
	return X * F ;
}

int N, M, l, r[MAXN], L = 1 ;
const double Pi = acos(-1.0);
struct complex {
    double x, y ;
    complex (double xx = 0, double yy = 0) {x = xx, y = yy ;}
} a[MAXN], b[MAXN] ;

complex operator + (complex a, complex b) {  
	return complex(a.x + b.x , a.y + b.y) ;
}
complex operator - (complex a, complex b) { 
	return complex(a.x - b.x , a.y - b.y) ;
}
complex operator * (complex a, complex b) { 
return complex(a.x * b.x - a.y * b.y , a.x * b.y + a.y * b.x) ;	
}
inline void FFT(complex *A, int type) {
	for (int i = 0 ; i < L ; i ++) if (i < r[i]) swap(A[i], A[r[i]]) ;
    for (int Mid = 1 ; Mid < L ; Mid <<= 1) {
    	complex W(cos(Pi / Mid), type * sin(Pi / Mid)) ;
    	for (int R = Mid << 1, j = 0 ; j < L ; j += R) {
    		complex w(1, 0) ;
    		for (int k = 0 ; k < Mid ;k ++, w = w * W) {
    			complex x = A[j + k], y = w * A[j + Mid + k] ;
    			A[j + k] = x + y ; A[j + Mid + k] = x - y ;
    		}
    	}
    }

}
int main() {
	int N = Read(), M = Read() ;
	for (int i = 0 ; i <= N ; i ++) a[i].x = Read() ;
	for (int i = 0 ; i <= M ; i ++) b[i].x = Read() ;
	while (L <= N + M) L <<= 1, l ++ ;
	for (int i = 0 ; i < L ; i ++)
		r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)) ;
	FFT (a, 1) ; FFT (b, 1) ;
	for (int i = 0 ; i <= L ; i ++) a[i] = a[i] * b[i] ; 
    FFT(a, - 1) ;
    for (int i = 0 ; i <= N + M ; i ++)
    	printf("%d ", (int)(a[i].x/ L + 0.5)) ; 
}

参考文献

同济高等数学，

Thomas H.Cormen算法导论，

Miskcoo 从多项式乘法到快速傅里叶变换，

郭晓旭 FFT讲义，

维基百科 傅里叶变换，

白霂凡 Fast Fourier Transform，

Menci FFT学习笔记，

Heinrich 傅里叶分析之掐死教程，

attack 快速傅里叶变换(FFT)详解