傅里叶变换

本文转载自查看原文 2021-09-23 10:18 733

1. 连续傅立叶变换（Continuous Fourier Transform）

对于时域连续函数 $x(t)$ ，它的傅立叶正变换（FT）定义为

$\mathcal F[ x(t)] = X(\omega) = \int_{-\infty}^\infty x(t) e^{-i\omega t}dt$ （用角频率 $\omega$ 表示) 或者 $\mathcal F[ x(t)] = X(f) = \int_{-\infty}^\infty x(t) e^{-i2\pi f t}dt$ （用频率 $f$ 表示, $\omega = 2\pi f$ )

傅立叶逆变换（inverse FT）定义为

$\mathcal F ^{-1}[ X(\omega)] = x(t) = \frac{1}{2\pi} \int_{-\infty}^\infty X(\omega) e^{+i\omega t}d\omega = \frac{1}{2\pi} \int_{-\infty}^\infty X(f) e^{+i2\pi f t}df$

2. 离散傅立叶变换（Discrete Fourier Transform, DFT）

如果有等时间间距 $\Delta t$ 的时域信号N个： $x[n],\text{where } n=0,\ldots,N-1$ ,它的离散傅立叶变换（DFT）定义为

$\mathcal F[x[n]]=X(e^{i\omega}) = X[k] = \sum_{n=0}^{N-1}x[n] e^{-i\frac{2\pi}{N}nk}, \;\;(k=0,\ldots,N-1)\\$

这里 $\omega = \frac{2\pi}{N} k$ . 相应的逆变换定义为

$\mathcal F^{-1}[X[k]]=X(e^{i\omega}) = x[n] = \frac{1}{N}\sum_{k=0}^{N-1}X[k] e^{i\frac{2\pi}{N}nk}, \;\;(n=0,\ldots,N-1)\\$

3. 快速傅立叶变换（Fast Fourier Transform, FFT）

快速傅立叶变换FFT是一种更加高效的计算离散傅立叶变换的算法，有着 $O(N log_2 N)$ 的计算复杂度，比原始的DFT $O(N^2)$ 计算复杂度有更好的可扩展性。Cooley-Tuckey算法是一种常见的FFT实现方式，但是需要限制输入数据点数 $N=2^m$ （2的乘方）形式。

如果你的数据数量不是2的几次幂，可以通过在后面加0（zero padding）的方式凑成下一个2的乘方个数据，添加的0不会影响频域谱，只会增加表观数据密度，即多出来的频域点其实是有效数据的线性内插值 Linear interpolation，如下图所示：

对于zero padding证明如下：

如果记原始时域信号的离散傅立叶变换为 $F(e^{i\omega}) = \sum_{n=0}^{M-1} f[n]e^{-1\omega n}$ ，其中时域数据数为 $M$ . 我们添加0到后面，记zero padding的信号为 $g[n]=f[n]+zeros, (M\le n \le N-1)$ , N 是2的乘方。所以 $g[n], (n=0,\ldots,N-1)$ 的离散傅立叶变换为 $G(e^{i\omega}) = \sum_{n=0}^{N-1} g[n]e^{-i\omega n}= \sum_{n=0}^{M-1} f[n]e^{-i\omega n}+ \sum_{n=M}^{N-1} 0\times e^{-i\omega n}=F(e^{i\omega}) \qquad \square\\$

离散傅立叶变换得到的频域信号是以 $\Delta f=\frac{1}{N\Delta t}$ 为间距的，最大频率称做Nyquist频率，即为 $f_{Nyquist}=f_{max} = \frac{1}{2\Delta t}$ ，它只跟时域信号的步长有关。经过FFT变换的频域信号对应的下列频率

$\begin{align} f=\frac{\omega}{2\pi} = & 0, \Delta f, 2\Delta f, \ldots, (N/2-1)\Delta f,\\ & \pm (N/2)\Delta f, -(N/2-1)\Delta f, -(N/2-2)\Delta f, \ldots, -\Delta f. \end{align}\\$

也就是说，前一半是正的频率，后一半是负的频率，中间点是正负最大频率。一般画图的话，需要把后一半挪到前面。

连续傅里叶变换

一般情况下，若“傅里叶变换”一词的前面未加任何限定语，则指的是“ 连续傅里叶变换”。“连续傅里叶变换”将平方可积的函数

表示成复指数函数的积分形式：

上式其实表示的是连续傅里叶变换的逆变换，即将时间域的函数表示为频率域的函数

的积分。反过来，其正变换恰好是将频率域的函数

表示为时间域的函数

的积分形式。一般可称函数

为原函数，而称函数

为傅里叶变换的像函数，原函数和像函数构成一个傅里叶变换对（transform pair）。

当

为奇函数（或偶函数）时，其余弦（或正弦）分量为零，而可以称这时的变换为余弦变换（或正弦变换）。

离散傅里叶变换

为了在科学计算和数字信号处理等领域使用计算机进行傅里叶变换，必须将函数定义在离散点上而非连续域内，且须满足有限性或周期性条件。这种情况下，序列

的离散傅里叶变换（discrete Fourier transform, DFT）为。

其逆变换为

直接使用DFT的定义计算的计算复杂度为

，而快速傅里叶变换（fast Fourier transform, FFT）可以将复杂度改进为

。计算复杂度的降低以及数字电路计算能力的发展使得DFT成为在信号处理领域十分实用且重要的方法。

如果写成一般形式，用 $x[n]$ 替换原信号，用 $N$ 替换40，得到

$X_k=\sum_{n=0}^{N-1}x[n]cos(2\pi\frac{kn}{N})$

问题来了，虽然貌似联系很紧密，但这怎么跟DFT的公式长得不一样。。。DFT的公式应该是这样的：

$X_k=\sum_{n=0}^{N-1}x[n]e^{-2\pi j\frac{kn}{N}}$

就算用欧拉公式展开，我们得到的是

$X_k=\sum_{n=0}^{N-1}x[n]cos(2\pi\frac{kn}{N})-j\sum_{n=0}^{N-1}x[n]sin(2\pi\frac{kn}{N})$

作者：肖畅
链接：https://www.zhihu.com/question/21314374/answer/542909849
来源：知乎
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

觉得现有答案都是从连续傅里叶变换引出离散的，各位大神功底很扎实啊。但问题是很多同学就是因为对那些繁琐的 $\int_{-\infty}^{\infty}$ 和 $e^{jwt}$ 感到晕厥，才转到没有积分，只有 $\sum$ ，看着比较简单的DFT来试试。所以DFT这样明显更简单的一个东西为什么一定要从连续傅里叶变换出发来讲呢？换个思路不行吗？

本文直接从DFT出发，不需要连续傅里叶变换的知识。

DFT做了什么事？

对一段有限长的离散信号，找出它含有的各个频率的正弦波分量。（这是句大家都知道的废话）

但如果换个单独针对DFT的解释就会好理解很多：

在长度为N的离散信号中，针对k=(0,1,2...)，分别找出在长度N内振动k个周期的三角波分量的权值

怎么做呢？举个例子，大家可以打开MATLAB一起做

抛开时间抛开采样频率，只看点的个数，我们对某个余弦信号在两个周期内采样了40次：

$x[n]=cos(2\pi \frac{n}{20}), n=0,1,\cdots,39$

n = 0:39; y = cos(2*pi*(2*n)/40); stem(n,y);

现在我们想通过DFT知道它在40次采样时间内振动了几个周期，算法是怎么处理的呢？其实很暴力：

事先选取40个长度为40个点的基信号，它们分别长这样：

第一个，40次采样内振动0个周期： $cos(2\pi\frac{0n}{40})$ ，即常值

40次采样内振动0个周期

第二个，40次采样内振动1个周期： $cos(2\pi \frac{n}{40})$

40次采样内振动1个周期

第三个，40次采样内振动2个周期： $cos(2\pi \frac{2n}{40})$

40次采样内振动2个周期

第四个，40次采样内振动3个周期： $cos(2\pi \frac{3n}{40})$

40次采样内振动3个周期

第五个，40次采样内振动4个周期： $cos(2\pi\frac{4n}{40})$

40次采样内振动4个周期

……

以此类推一直到40个采样内振动39个周期。

接下来，对于上述40个基信号，判断它们跟原信号的相关程度。代数里判断两个长度相同的离散函数x[k],y[k]的相关程度，是用把它们在同一点的函数值相乘，再对得到的积求和的方法（就是向量求内积的方法），即下面的公式：

$correlation(x,y) = \sum_k{x[k]y[k]}$

这个值越大，x[k]和y[k]长得越像。于是DFT把

$cos(2\pi\frac{0n}{40})$ 到 $cos(2\pi\frac{39n}{40})$ 这40个基函数都跟当前函数 $cos(2\pi\frac{2n}{40})$ 比较了一下，发现

$cos(2\pi \frac{2n}{40})$ 和 $cos(2\pi \frac{38n}{40})$ 跟它长得最像！

这很显然，因为 $cos(2\pi\frac{38n}{40})=cos(2\pi n-2\pi\frac{38n}{40})=cos(2\pi\frac{2n}{40})$

于是 $X_2=X_{38}=\sum_{n=0}^{39}cos(2\pi\frac{2n}{40})cos(2\pi\frac{2n}{40})=20$

下面，如果我们把这40次每次比较的correlation值记下来，就得到了原信号在每个频率上的分量大小，比如，我要看它跟40次采样中振动8次的那个基信号相似程度，就是：

把这40个correlation值记下来，就得到了原信号的频域X：

stem(n, abs(fft(y)));

可以看到 $X_2 = X_{38} = 20$ ，其他值都是0

如果写成一般形式，用 $x[n]$ 替换原信号，用 $N$ 替换40，得到

$X_k=\sum_{n=0}^{N-1}x[n]cos(2\pi\frac{kn}{N})$

问题来了，虽然貌似联系很紧密，但这怎么跟DFT的公式长得不一样。。。DFT的公式应该是这样的：

$X_k=\sum_{n=0}^{N-1}x[n]e^{-2\pi j\frac{kn}{N}}$

就算用欧拉公式展开，我们得到的是

$X_k=\sum_{n=0}^{N-1}x[n]cos(2\pi\frac{kn}{N})-j\sum_{n=0}^{N-1}x[n]sin(2\pi\frac{kn}{N})$

这又是为什么？

这是因为，对于一个信号，如果只跟余弦函数比较，会损失一些信息，比如相位。

如果我们的原信号有一些相位偏移， $x=cos(2\pi \frac{2n}{40} + \frac{\pi}{4})$

如果对这个函数同样按照上面的方法计算频域，结果会有一些不一样：

n = 0:39; x = cos(2 * pi * (2 * n)/40 + pi/4); stem(n, real(fft(x)));

虽然对应频率还是筛选出来了，但是权值大小不一样了，如果计算一下得知

$X_2=X_{38}=\sum_{n=0}^{39}cos(2\pi\frac{2n}{40}+\frac{\pi}{4})cos(2\pi\frac{2n}{40})=10\sqrt{2}$

在相位偏移前，这个值是20，是原来的 $\frac{\sqrt{2}}{2}$ 倍

注意如果我们再找一个信号 $y$ ，没有相位偏移，但是直接把幅值砍到 $\frac{\sqrt{2}}{2}$ ，即：

$y=\frac{\sqrt{2}}{2}cos(2\pi \frac{2n}{40})$

那么这个信号的

$Y_2=Y_{38}=\sum_{n=0}^{39}\frac{\sqrt{2}}{2}\cos(2\pi\frac{2n}{40})\cos(2\pi\frac{2n}{40})=10\sqrt{2}$

跟 $x$ 信号的结果一模一样，这样我们就无法通过频域恢复信号了，或者说我们损失了信息

解决方法是另选一组以正弦函数（实际上选了负正弦）为基准的“基信号”，即 $-\sin(2\pi \frac{0n}{40})$ 到 $-\sin(2\pi \frac{39n}{40})$ ，计算另一组原信号与正弦基的相关系数，这两组系数一起作为DFT的最终结果。而复数只是一个工具，用来方便地同时存储两组计算结果。当然它还有一个好处就是能够比较直观地表现出模和相位。

选负正弦还是正弦做基信号其实无所谓，只是最后的结果算出来的相位反一下而已，幅值是一样的。如果一个信号跟某频率余弦和负正弦的相关系数分别为a,b，那么代表这个信号差不多型如

$a\cos(\frac{kn}{N})-b\sin(\frac{kn}{N})$

根据高中数学可以求得其模为 $\sqrt{a^2+b^2}$ ，（相对余弦的）相位为 $\arctan(\frac{b}{a})$ ，这与复数 $a+bi$ 的模和相位是相同的，因此DFT的公式

$X_k=\sum_{n=0}^{N-1}x[n]e^{-2\pi j\frac{kn}{N}}$

相当于同时把x[n]做了跟N个余弦基、N个负正弦基的比对，结果用N个复数存储。如果想要看某个频率的相位和模，就看对应复数的相位和模。

我们再来看看上面有相位偏移的那个例子：

原信号： $x[n]=cos(2\pi\frac{2n}{40}+\frac{\pi}{4})$

与余弦比对： $X_2 =\sum_{n=0}^{39} cos(2\pi\frac{2n}{40}+\frac{\pi}{4})cos(2\pi\frac{2n}{40})=10\sqrt{2}$

与负正弦比对： $X_2 =\sum_{n=0}^{39} cos(2\pi\frac{2n}{40}+\frac{\pi}{4})sin(-2\pi\frac{2n}{40})=10\sqrt{2}$

在40个点内振动两个周期这个频率上，其DFT的结果为 $10\sqrt{2}+10\sqrt{2}j$

其模为20，与相位偏移前相同。

其相位为 $\frac{\pi}{4}$ ，也没有问题。

如果将原信号变为 $x[n]=sin(2\pi\frac{2n}{40}+\frac{\pi}{4})$ ，会求得该频率DFT结果为 $10\sqrt{2}-10\sqrt{2}j$ ，求得其相位为 $-\frac{\pi}{4}$ 。因此，根据DFT结果求得的相位是相对余弦信号的相位。

对于实际信号

上面做的分析没有考虑时间和采样频率，但只要给定一个采样频率，时域的横轴就可以换算成时间，频域的横轴就可以换算成频率。

DFT所做的只是选取了在给定长度内振动了整数个周期的正弦和余弦波作为基。

上面的信号是40个采样，如果给一个采样频率是100Hz，那么信号长度就是0.4s，原信号在40个采样内振动了两个周期，可以算出其频率为5Hz

对于频域，每个“在40个采样内振动了k个周期”的基信号的实际频率为 $\frac{100}{40}k=2.5k$ ，也就是说频域图中的每一个点代表2.5Hz，这个系统的频域分辨率是2.5Hz，所以我们从 $X_2=20$ 得知原信号的实际频率为 $2\times2.5=5Hz$

FFT(Fast Fourier Transform)是什么

FFT是 Cooley & Tuket 两人1965年提出的快速计算DFT的算法。这背后还有个故事，美国和苏联1963年签了个核试验禁令，互相约定大家都不搞核试验了。但是美国不放心啊，怕毛子说一套做一套，肯尼迪就请了一堆科学家开会，说想搞一套不用去苏联检查就能探测到核试验的设备。美国在苏联周围放了一圈声波探测仪，但是遇到个问题——DFT算得太慢了。正巧 Tuket 和IBM的一个叫 Richard Garwin 的出席了那次会议。会上Garwin就想起来Tuket和Cooley好像合作过一个类似的算法，就让Tuket去找Cooley，但是隐藏了真实目的，告诉Cooley这个算法是为了探测氦3晶体的自旋周期……（也是醉了orz）。两人的算法最后于1965年发表，极大提升了DFT的计算速度，甚至被后世誉为“20世纪最伟大的工程算法”。

从DFT的公式可以看出其算法复杂度是 $O(n^2)$ 的，两人的FFT算法表明，在N为合数的情况下 $N=n_1n_2$ ，原长度为N的DFT可以分解为两个长度分别为 $n_1$ 和 $n_2$ 的DFT！

当然算法不停止于此，如果 $n_1$ 和 $n_2$ 还可以继续分解，就可以分解为更多小DFT来计算。

使用了算法中的分治（divide and conquer）策略。

而由于DFT的长度其实是可以任意选择的，因此通常在使用FFT时，都选择2的整数次幂作为FFT长度，这样可以一直分解到N/2个长度为2的DFT，复杂度直接降到 $O(nlogn)$

更详细的解答，请移步我做的视频：https://www.bilibili.com/video/av44600709

想从理论层面更深入了解DFT的，可以看我另外一个视频：https://www.bilibili.com/video/BV1a

0x00 写在前面

为了让更多人能够看到这个教程，希望大家收藏之前，也要点赞哦！！！蟹蟹大家的认可和鼓励。

傅里叶变换

快速傅里叶变换（Fast Fourier Transform，FFT）是一种可在 $O(nlogn)$ 时间内完成的离散傅里叶变换（Discrete Fourier transform，DFT）算法。

在算法竞赛中的运用主要是用来加速多项式的乘法。

考虑到两个多项式 $A(x),B(x)$ 的乘积 $C(x)$ ，假设 $A(x)$ 的项数为 $n$ ，其系数构成的 $n$ 维向量为 $(a_0,a_1,a_2,...,a_{n-1})$ ， $B(x)$ 的项数为 $m$ ，其系数构成的 $m$ 维向量为 $(b_0,b_1,b_2,...,b_{m-1})$ 。

我们要求 $C(x)$ 的系数构成的 $n+m-1$ 维的向量，先考虑朴素做法。

可以用这段代码表示：

for ( int i = 0 ; i < n ; ++ i ) for ( int j = 0 ; j < m ; ++ j ) { c [i + j] += a [i] * b [j] ; }

思路非常清晰，其时间复杂度是 $O(n^2)$ 的。

所以我们来学习快速傅里叶变换。

0x01 关于多项式

多项式有两种表示方法，系数表达法与点值表达法

多项式的系数表示法

设多项式 $A(x)$ 为一个 $n-1$ 次的多项式，显然，所有项的系数组成的系数向量 $(a_0,a_1,a_2,...,a_{n-1})$ 唯一确定了这个多项式。

$A(x)=\sum_{i=0}^{n-1}a_i\cdot x^i$

多项式的点值表示法

将一组互不相同的 $(x_0,x_1,x_2,...,x_n)$ （叫插值节点）分别带入 $A(x)$ ，得到 $n$ 个取值 $(y_0,y_1,y_2,...,y_n)$ .

其中

$y_i=\sum_{j=0}^{n-1}a_j\cdot x_i^j$

定理:

一个 $n-1$ 次多项式在 $n$ 个不同点的取值唯一确定了该多项式。

证明：

假设命题不成立，存在两个不同的 $n-1$ 次多项式 $A(x),B(x)$ ，满足对于任何 $i \in [0,\ n - 1]$ ，有 $A(x_i) = B(x_i)$ 。
令 $C(x) = A(x) - B(x)$ ，则 $C(x)$ 也是一个 $n-1$ 次多项式。对于任何 $i \in [0,\ n - 1]$ ，都有 $C(x_i) = 0$ 。
即 $C(x)$ 有 $n$ 个根，这与代数基本定理（一个 $n-1$ 次多项式在复数域上有且仅有 $n-1$ 个根）相矛盾，故 $C(x)$ 并不是一个 $n-1$ 次多项式，推到矛盾。
原命题成立，证毕。

如果我们按照定义求一个多项式的点值表示，时间复杂度为 $O(n^2)$

已知多项式的点值表示，求其系数表示，可以使用插值。朴素的插值算法时间复杂度为 $O(n^2)$ 。

关于多项式的乘法

已知在一组插值节点 $(x_0,x_1,x_2,...,x_n)$ 中 $A(x),B(x)$ (假设个多项式的项数相同，没有的视为 $0$ )的点值向量分别为 $(y_{a0},y_{a1},y_{a_2},...,y_{an}),(y_{b0},y_{b1},y_{b_2},...,y_{bn})$ ，那么 $C(x)=A(x)\cdot B(x)$ 的点值表达式可以在 $O(n)$ 的时间内求出，为 $(y_{a0}\cdot y_{b0},y_{a1}\cdot y_{b1},y_{a_2}\cdot y_{b2},...,y_{an}\cdot y_{bn})$ 。

因为 $C(x)$ 的项数为 $A(x),B(x)$ 的项数之和。

设 $A(x),B(x)$ 分别有 $n,m$ 项所以我们带入的插值节点有至少有 $n+m$ 个。

如果我们能快速通过点值表式求出系数表示，那么就搭起了它们之间的一座桥了。

这也是快速傅里叶变换的基本思路，由系数表达式到点值表达式到结果的点值表达式再到结果的系数表达式。

0x02 关于复数的基本了解

我们把形如 $a+bi$ 这样的数叫做复数，复数集合用 $C$ 来表示。其中 $a$ 称为实部 $(real\;part)$ ， $b$ 称为虚部 $(imaginary\;part)$ , $i$ 为虚数单位，指满足 $x^2=-1$ 的一个解 $\sqrt{-1}$ ；此外，对于这样对复数开偶次幂的数叫做虚数 $( imaginary\;number)$ .

每一个复数 $a+bi$ 都对应了一个平面上的向量 $(a,b)$ 我们把这样的平面称为复平面 $(complex\;plane)$ ，它是由水平的实轴与垂直的虚轴建立起来的复数的几何表示。

故每一个复数唯一对应了一个复平面上的向量，每一个复平面上的向量也唯一对应了一个复数。其中 $0$ 既被认为是实数，也被认为是虚数。

其中复数 $z=a+bi$ 的模长 $\left| z \right|$ 定义为 $z$ 在复平面的距离到原点的距离， $\left| z \right|=\sqrt{a^2+b^2}$ 。幅角 $\theta$ 为实轴的正半轴正方向（逆时针）旋转到 $z$ 的有向角度。

由于虚数无法比较大小。复数之间的大小关系只存在等于与不等于两种关系，两个复数相等当且仅当实部虚部对应相等。对于虚部为 $0$ 的复数之间是可以比较大小的，相当于实数之间的比较。

复数之间的运算满足结合律，交换律和分配律。

由此定义复数之间的运算法则：

$(a+bi)+(c+di)=(a+c)+(b+d)i$

$(a+bi)-(c+di)=(a-c)+(b-d)i$ $(a+bi)\cdot(c+di)=ac+adi+bci+bdi^2=(ac-bd)+(ad+cb)i$

$\frac{a+bi}{c+di}=\frac{(a+bi)\cdot(c-di)}{(c+di)\cdot(c-di)}=\frac{(ac+bd)+(bc-ad)i}{c^2+d^2}=\frac{(ac+bd)}{c^2+d^2}+\frac{(bc-ad)i}{c^2+d^2}$

复数运算的加法满足平行四边形法则，乘法满足幅角相加，模长相乘。

对于一个复数 $z=a+bi$ ,它的共轭复数是 $z^{'}=a-bi$ , $z^{'}$ 称为 $z$ 的复共轭 $(complex\;conjugate)$ .

共轭复数有一些性质

$z\cdot z^{'}=a^2+b^2$

$\left| z \right| = \left | z^{'} \right|$

0x03 复数中的单位根

复平面中的单位圆

其中 $\bar{OA}$ 单位根，表示为 $e^{i\theta}$ ，可知 $e^{i\theta}=cos\theta+i\cdot sin\theta$

（顺便一提著名的欧拉幅角公式 $e^{i\pi}+1=0$ 其实是由定义来的...）

将单位圆等分成 $n$ 个部分（以单位圆与实轴正半轴的交点一个等分点），以原点为起点，圆的这 $n$ 个 $n$ 等分点为终点，作出 $n$ 个向量。

其中幅角为正且最小的向量称为 $n$ 次单位向量，记为 $\omega_{n}^{1}$ 。

（有没有大佬帮我补张图啊，画不来）

其余的 $n-1$ 个向量分别为 $\omega_{n}^{2},\omega_{n}^{3},......,\omega_{n}^{n}$ ，它们可以由复数之间的乘法得来 $w_{n}^{k}=w_{n}^{k-1}\cdot w_{n}^{1}\ (2 \leq k \leq n)$ 。

容易看出 $w_{n}^{n}=w_{n}^{0}=1$ 。

对于 $w_{n}^{k}$ ,它事实上就是 $e^{2\pi\cdot \frac{k}{n}i}$ 。

所以 $\omega_{n}^{k}=cos(2\pi\cdot\frac{k}{n})+i\cdot sin(2\pi\cdot\frac{k}{n})$

关于单位根有两个性质

性质一（又称为折半引理）：

$\omega_{2n}^{2k}=\omega_{n}^{k}$

证明一：

由几何意义，这两者表示的向量终点是一样的。

证明二：

由计算的公式：

$\omega_{2n}^{2k}=cos(2\pi\frac{2k}{2n})+i\cdot sin(2\pi\frac{2k}{2n})=cos(2\pi\frac{k}{n})+i\cdot sin(2\pi\frac{k}{n})=\omega_{n}^{k}$

其实由此我们可以引申出

$\omega_{mn}^{mk}=\omega_{n}^{k}$

性质二（又称为消去引理）

$\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}$

证明一：

由几何意义，这两者表示的向量终点是相反的，左边较右边在单位圆上多转了半圈。

证明二：

由计算的公式：

$\omega_{n}^{k+\frac{n}{2}}=cos(2\pi\frac{k+\frac{n}{2}}{n})+i\cdot sin(2\pi\frac{k+\frac{n}{2}}{n})=cos(2\pi\frac{k}{n}+\pi)+i\cdot sin(2\pi\frac{k}{n}+\pi)=-cos(2\pi\frac{k}{n})-i\cdot sin(2\pi\frac{k}{n})=-\omega_{n}^{k}$

最后一步由三角恒等变换得到。

0x04 离散傅里叶变换（Discrete Fourier Transform）

首先我们单独考虑一个 $n$ 项（ $n=2^x$ ）的多项式 $A(x)$ ，其系数向量为 $(a_0,a_1,a_2,...,a_{n-1})$ 。我们将 $n$ 次单位根的 $0$ ~ $n-1$ 次幂分别带入 $A(x)$ 得到其点值向量 $(A(w_n^{0}),A(w_n^{1}),A(w_n^{2}),...,A(w_n^{n-1}))$ 。

这个过程称为离散傅里叶变换（Discrete Fourier Transform）。

如果朴素带入，时间复杂度也是 $O(n^2)$ 的。

所以我们必须要利用到单位根 $\omega$ 的特殊性质。

对于 $A(x)=a_0+a_1\cdot x^1+a_2\cdot x^2+a_3\cdot x^3+...+a_{n-1}\cdot x^{n-1}$

考虑将其按照奇偶分组

$A(x)=(a_0+a_2\cdot x^2+a_{4}\cdot x^{4}...+a_{n-2}\cdot x^{n-2})+(a_1\cdot x^1+a_3\cdot x^3+a_{5}\cdot x^{5}+...+a_{n-1}\cdot x^{n-1})$

$A(x)=(a_0+a_2\cdot x^2+a_{4}\cdot x^{4}...+a_{n-2}\cdot x^{n-2})+x\cdot(a_1+a_3\cdot x^2+a_{5}\cdot x^{4}+...+a_{n-1}\cdot x^{n-2})$

令

$A1(x)=(a_0+a_2\cdot x+a_{4}\cdot x^{2}...+a_{n-2}\cdot x^{\frac{n-2}{2}})$

$A2(x)=(a_1+a_3\cdot x+a_{5}\cdot x^{2}...+a_{n-1}\cdot x^{\frac{n-2}{2}})$

则可得到

$A(x)=A1(x^2)+x\cdot A2(x^2)$

分类讨论

设 $0\leq k\leq \frac{n}{2}-1$ ， $k\in Z$

$A(\omega_{n}^{k})=A1(\omega_{n}^{2k})+\omega_{n}^{k}\cdot A2(\omega_{n}^{2k})$

由上文提到的折半引理

$A(\omega_{n}^{k})=A1(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k}\cdot A2(\omega_{\frac{n}{2}}^{k})$

对于 $\frac{n}{2}\leq k+\frac{n}{2}\leq n-1$

$A(\omega_{n}^{k+\frac{n}{2}})=A1(\omega_{n}^{2k+n})+\omega_{n}^{k+\frac{n}{2}}\cdot A2(\omega_{n}^{2k+n})$

其中 $\omega_{n}^{2k+n}=\omega_{n}^{2k}\cdot \omega_{n}^{n}=\omega_{n}^{2k}=\omega_{\frac{n}{2}}^{k}$

由消去引理 $\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}$

故 $A(\omega_{n}^{k+\frac{n}{2}})=A1(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}\cdot A2(\omega_{\frac{n}{2}}^{k})$

注意， $k$ 与 $k+\frac{n}{2}$ 取遍了 $[0,n-1]$ 中的 $n$ 个整数，保证了可以由这 $n$ 个点值反推解出系数（上文已证明）。

于是我们可以知道

如果已知了 $A1(x),A2(x)$ 分别在 $\omega_{\frac{n}{2}}^{0},\omega_{\frac{n}{2}}^{1},...,\omega_{\frac{n}{2}}^{\frac{n}{2}-1},$ 的取值，可以在 $O(n)$ 的时间内求出 $A(x)$ 的取值。

而 $A1(x),A2(x)$ 都是 $A(x)$ 一半的规模，显然可以转化为子问题递归求解。

时间复杂度：

$T(n)=2T(\frac{n}{2})+O(n)=O(nlogn)$

0x05 离散傅里叶反变换（Inverse Discrete Fourier Transform）

使用快速傅里叶变换将点值表示的多项式转化为系数表示，这个过程叫做离散傅里叶反变换（Inverse Discrete Fourier Transform）。

即由 $n$ 维点值向量 $(A(x_0),A(x_1),...,A(x_{n-1}))$ 推出 $n$ 维系数向量 $(a_0,a_1,...,a_{n-1})$ 。

设 $(d_0,d_1,...,d_{n-1})$ 为 $(a_0,a_1,...,a_{n-1})$ 得到的离散傅里叶变换的结果。

我们构造一个多项式 $F(x)=d_0+d_1\cdot x+d_2\cdot x^2+...+d_{n-1}\cdot x^{n-1}$

设向量 $(c_0,c_1,...,c_{n-1})$ 中

$c_k$ 为 $F(x)$ 在 $x=\omega_{n}^{-k}$ 的点值表示

即 $c_k=\sum_{i=0}^{n-1}d_i\cdot(\omega_{n}^{-k})^i$ ，

我们考虑对 $d_i$ 进行还原

于是

$c_k=\sum_{i=0}^{n-1}[\sum_{j=0}^{n-1}a_j\cdot(\omega_{n}^{i})^j]\cdot(\omega_{n}^{-k})^i$

由和式的性质

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

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

令 $S(j,k)=\sum_{i=0}^{n-1}(\omega_{n}^{i})^{j-k}$

对其进行化简

设 $j-k=\delta$

则 $S(j,k)=\omega_{n}^{0}+\omega_{n}^{\delta}+\omega_{n}^{2\delta}+...+\omega_{n}^{(n-1)\delta}$

其公比为 $\omega_{n}^{\delta}$

当 $\omega_{n}^{\delta}=1$ 即 $\delta=0$ 时

$S(j,k)=n$ 此时 $\delta=0\Rightarrow j-k=0 \Rightarrow j=k$

当 $\omega_{n}^{\delta}\ne1$ 即 $\delta\ne0$ 时

由等比数列求和公式

$S(j,k)=\frac{\omega_{n}^{0}[(\omega_{n}^{\delta})^{n}-1]}{\omega_{n}^{\delta}-1}=\frac{1 [(\omega_{n}^{n})^{\delta}-1]}{\omega_{n}^{\delta}-1}=\frac{[(1)^{\delta}-1]}{\omega_{n}^{\delta}-1}=\frac{0}{\omega_{n}^{\delta}-1}=0$ ,此时 $j\ne k$ .

所以

$S(j,k)=[j=k]\cdot n$

将 $S(j,k)$ 带入原式

$c_k=\sum_{j=0}^{n-1}a_j\cdot S(j,k)=\sum_{j=0}^{n-1}a_j\cdot [j=k]\cdot n=a_k\cdot n$

所以 $a_k=\frac{c_k}{n}$ .

其中 $a_k$ 为原多项式 $A(x)$ 的系数向量 $(a_0,a_1,...,a_n)$ 中的 $a_k$ .

由此得到：

对于多项式 $A(x)$ 由插值节点 $(\omega_{n}^{0},\omega_{n}^{1},\omega_{n}^{2},...,\omega_{n}^{n-1})$ 做离散傅里叶变换得到的点值向量 $(d_0,d_1,...,d_{n-1})$ 。我们将 $(\omega_{n}^{0},\omega_{n}^{-1},\omega_{n}^{-2},...,\omega_{n}^{-(n-1)})$ 作为插值节点， $(d_0,d_1,...,d_{n-1})$ 作为系数向量，做一次离散傅里叶变换得到的向量每一项都除以 $n$ 之后得到的 $(\frac{c_0}{n},\frac{c_1}{n},...,\frac{c_{n-1}}{n})$ 就是多项式的系数向量 $(a_0,a_1,...,a_{n-1})$ 。

注意到 $\omega_{n}^{-k}$ 是 $\omega_{n}^{k}$ 的共轭复数。

这个过程称为离散傅里叶反变换。

0x06 关于FFT在C++的实现

首先要解决复数运算的问题，我们可以使用C++STL自带的 $std :: complex < T >$ 依照精度要求 $T$ 一般为 $double,long\;double$ 。

也可以自己封装，下面是我封装的复数类。

struct Complex { double r, i ; Complex ( ) { } Complex ( double r, double i ) : r ( r ), i ( i ) { } inline void real ( const double& x ) { r = x ; } inline double real ( ) { return r ; } inline Complex operator + ( const Complex& rhs ) const { return Complex ( r + rhs.r, i + rhs.i ) ; } inline Complex operator - ( const Complex& rhs ) const { return Complex ( r - rhs.r, i - rhs.i ) ; } inline Complex operator * ( const Complex& rhs ) const { return Complex ( r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r ) ; } inline void operator /= ( const double& x ) { r /= x, i /= x ; } inline void operator *= ( const Complex& rhs ) { *this = Complex ( r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r ) ; } inline void operator += ( const Complex& rhs ) { r += rhs.r, i += rhs.i ; } inline Complex conj ( ) { return Complex ( r, -i ) ; } } ;

我们由上面的分析可以得到这个递归的写法。

bool inverse = false ; inline Complex omega ( const int& n, const int& k ) { if ( ! inverse ) return Complex ( cos ( 2 * PI / n * k ), sin ( 2 * PI / n * k ) ) ; return Complex ( cos ( 2 * PI / n * k ), sin ( 2 * PI / n * k ) ).conj ( ) ; } inline void fft ( Complex *a, const int& n ) { if ( n == 1 ) return ; static Complex buf [N] ; const int m = n >> 1 ; for ( int i = 0 ; i < m ; ++ i ) { buf [i] = a [i << 1] ; buf [i + m] = a [i << 1 | 1] ; } memcpy ( a, buf, sizeof ( Complex ) * ( n + 1 ) ) ; Complex *a1 = a, *a2 = a + m; fft ( a1, m ) ; fft ( a2, m ) ; for ( int i = 0 ; i < m ; ++ i ) { Complex t = omega ( n, i ) ; buf [i] = a1 [i] + t * a2 [i] ; buf [i + m] = a1 [i] - t * a2 [i] ; } memcpy ( a, buf, sizeof ( Complex ) * ( n + 1 ) ) ; }

但是这样的 $FFT$ 要用到辅助数组，并且常数比较大。

能不能优化呢？

我们把每一次分组的情况推演出来

递归分类的每一层

观察到每一个位置的数其实都是原来位置上的数的二进制后 $log_{2}n$ 位 $reverse$ 了一下。

于是我们可以想，先将原数组调整成最底层的位置（很好调整吧）。

然后从倒数第二层由底向上计算。

这就是我们一般用来实现 $FFT$ 的 $Cooley-Tukey$ 算法。

考虑怎么合并？

在 $Cooley-Tukey$ 算法中，合并操作被称作是蝴蝶操作。

虑合并两个子问题的过程，这一层有 $n$ 项需要处理。假设 $A_1(\omega_{ \frac{n}{2} } ^ k)$ 和 $A_2(\omega_{ \frac{n}{2} } ^ k)$ 分别存在 $a(k)$ 和 $a(\frac{n}{2} + k)$ 中， $A(\omega_n ^ {k})$ 和 $A(\omega_n ^ {k + \frac{n}{2} })$ 将要被存放在 $buf(k)$ 和 $buf(\frac{n}{2} + k)$ 中，合并的单位操作可表示为

$buf(k):=a(k)+\omega_{n}^{k}\cdot a(k+\frac{n}{2})$

$buf(k+\frac{n}{2}):=a(k)-\omega_{n}^{k}\cdot a(k+\frac{n}{2})$

只要将合并顺序换一下，再加入一个临时变量，合并过程就可以在原数组中进行。

令 $t:=\omega_{n}^{k}\cdot a(k+\frac{n}{2})$

合并过程如下：

$a(k+\frac{n}{2}):=a(k)-t$

$a(k):=a(k)+t$ 。

至此，我们可以给出 $Cooley-Tukey$ 算法的实现。

struct FastFourierTransform { Complex omega [N], omegaInverse [N] ; void init ( const int& n ) { for ( int i = 0 ; i < n ; ++ i ) { omega [i] = Complex ( cos ( 2 * PI / n * i), sin ( 2 * PI / n * i ) ) ; omegaInverse [i] = omega [i].conj ( ) ; } } void transform ( Complex *a, const int& n, const Complex* omega ) { for ( int i = 0, j = 0 ; i < n ; ++ i ) { if ( i > j ) std :: swap ( a [i], a [j] ) ; for( int l = n >> 1 ; ( j ^= l ) < l ; l >>= 1 ) ; } for ( int l = 2 ; l <= n ; l <<= 1 ) { int m = l / 2; for ( Complex *p = a ; p != a + n ; p += l ) { for ( int i = 0 ; i < m ; ++ i ) { Complex t = omega [n / l * i] * p [m + i] ; p [m + i] = p [i] - t ; p [i] += t ; } } } } void dft ( Complex *a, const int& n ) { transform ( a, n, omega ) ; } void idft ( Complex *a, const int& n ) { transform ( a, n, omegaInverse ) ; for ( int i = 0 ; i < n ; ++ i ) a [i] /= n ; } } fft ;