这篇文章参考了 知乎文章, 作者清华爷 @cronus-裁

傅里叶级数回忆

Bf={12π,1πcos⁡(nx),1πsin⁡(nx)}n=1∞B_f = \{\frac 1 {\sqrt{2\pi}},\frac 1{\sqrt{\pi}}\cos(nx),\frac 1{\sqrt\pi}\sin(nx)\}_{n=1}^\infty 是平方可积空间 L2([−π,π])L^2([-\pi,\pi]) 的一组函数的基，对于空间中的任意函数，我们都有 f=lim⁡SNf = \lim S_N
拟合公式 SN(x)=⟨f,1⟩2π+Σn=1N⟨cos⁡(nx),f⟩πcos⁡(nx)+Σn=1N⟨sin⁡(nx),f⟩πsin⁡(nx)S_N (x) = \frac{\langle f,1\rangle}{2\pi}+\Sigma_{n=1}^N \frac{\langle \cos(nx),f\rangle}{\pi}\cos(nx)+\Sigma_{n=1}^N\frac{\langle \sin(nx),f\rangle}{\pi}\sin(nx)
从上述这个公式我们可以看出，傅里叶级数是将函数投影到不同频率的 sin⁡(nx)\sin(nx) 和 cos⁡(nx)\cos(nx) 的空间上的，或者在信号领域，我们可以将这个 nn 写作 nωn\omega 表示基频的整数倍频率空间
分解的过程中，对于每一个 nω0n\omega_0 ，我们都得到了对应的幅值，这就组成了一个函数关系（自变量为 nω0n\omega_0 ，因变量为幅值，即相应频率信号的强度）我们称之为频谱函数。

傅里叶变换回忆

F(f(t))=F(ω)=∫−∞∞​f(t)e−iωtdt\mathcal{F}(f(t))=F(\omega)=\int_{-\infty}^\infty ​f(t)e^{-i\omega t}dt

这个公式表示将一个时域空间的函数通过傅里叶变换变成了 频域函数，表示在某一个频率处的强度

傅里叶变换和傅里叶级数的联系

傅里叶级数是傅里叶变换的离散形式，而傅里叶变换是傅里叶级数在周期趋于无穷大时的极限情况

离散傅里叶变换 DFT

傅里叶变换的目的就是得到信号的频谱密度函数（自变量是 𝜔 ，因变量是信号幅值在频域中的分布密度，即单位频率信号的强度），它实际上揭示了信号的强度和频率的关系

采样

首先我们说明狄拉克单位冲激函数对于采样的描述：
 $$\int_{-\infty}^{\infty}\delta(t-t_0)f(t)dt=f(t_0)$$
 那么我们在时域上将每一个通过狄拉克函数进行筛选
 的点叠加起来就可以得到时域上的采样函数了，我们称之为 梳状函数，表达式

\delta_s(t)=\sum_{n=-\infty}^\infty \delta(t-nT_s)$$ 那么我们用时域上的连续信号与之相乘（其实是离散形式的卷积）就得到了  $$f_s=\sum_{n=-\infty}^\infty f(t)\delta(t−nT_s)

时域离散化计算

对于采样得到的 xs(t)=∑n=−∞∞x(t)δ(t−nTs)x_s(t)=\sum_{n=-\infty}^\infty x(t)\delta(t−nT_s) 进行傅里叶变换。
根据傅里叶变换的定义:

X(ω)=∫−∞+∞x(t)e−jωtdtX(\omega)=\int_{-\infty}^{+\infty}x(t)e^{-j\omega t}dt

则 𝑋_𝑠(𝜔)=∫_{−∞}^∞(∑_{𝑛=−∞}^∞ 𝑥(𝑡)𝛿(𝑡−𝑛𝑇_𝑠))𝑒^{−𝑗𝜔𝑡}𝑑𝑡, 交换积分与求和顺序(这里用到了 Fubini 定理或者绝对收敛级数可交换性)，得到 𝑋_𝑠(𝜔)=∑_{𝑛=−∞}^∞∫_{−∞}^∞𝑥(𝑡)𝛿(𝑡−𝑛𝑇_𝑠)𝑒^{−𝑗𝑤𝑡}𝑑𝑡, 根据Dirac函数的筛选特性，

𝑋_𝑠(𝜔)=\sum_{𝑛=−∞}^∞𝑥(𝑛𝑇_𝑠)𝑒^{−𝑗𝑤𝑛𝑇_𝑠}

但是不幸的是，这个函数对于 ω\omega 还是一个连续函数

频域离散化计算

对于连续信号 𝑥(𝑡)𝑥(𝑡) 进行 NN 次（N为有限值）采样，采样周期为 𝑇𝑠𝑇_𝑠 。然后对采样得到的信号进行时域上的周期延拓（延拓至正负无穷），这样我们就得到了一个周期为 𝑇0=𝑁𝑇𝑠𝑇_0=𝑁𝑇_𝑠 的函数。对于周期函数而言，它的频谱密度函数是离散化的，这样我们就成功把频域也进行了离散化。
在一个周期内，离散信号的表达式为 𝑥_𝑠(𝑡)=∑_{𝑛=0}^{𝑁−1}𝑥(𝑡)𝛿(𝑡−𝑛𝑇_𝑠)
离散信号的傅里叶变换 𝑋(𝑘𝜔_0)=\frac{1}{𝑇_0}∫_0^𝑇(∑_{𝑛=0}^{𝑁−1}𝑥(𝑡)𝛿(𝑡−𝑛𝑇_𝑠))𝑒^{−𝑗𝑘𝑤_0𝑡}𝑑𝑡
其中， 𝜔_0=2𝜋 / 𝑇_0
交换积分和求和次序， 𝑋(𝑘𝜔_0)=\frac{1}{𝑇_0}∑_{𝑛=0}^{𝑁−1}∫_0^𝑇𝑥(𝑡)𝛿(𝑡−𝑛𝑇_𝑠)𝑒^{−𝑗𝑘𝜔_0𝑡}𝑑𝑡
根据狄拉克函数的筛选特性

𝑋(𝑘𝜔_0)=\frac{1}{𝑁𝑇_𝑠}∑_{𝑛=0}^{𝑁−1}𝑥(𝑛𝑇_𝑠)𝑒^{−𝑗\frac{2𝜋}{𝑁}𝑘𝑛}

为更加简明，令 𝑋[𝑘]=𝑋(𝑘𝜔_0)𝑇_𝑠 ， 𝑥[𝑛]=𝑥(𝑛𝑇_𝑠) ，则

𝑋[𝑘]=\frac{1}{𝑁}∑_{𝑛=0}^{𝑁−1}𝑥[𝑛]𝑒^{−𝑗\frac{2𝜋}{𝑁}𝑘𝑛}

复杂度分析

首先对于变量 X[k]X[k]，我们需要计算 NN 次的 x[n]x[n] 和 e−j2πNkne^{-j\frac{2\pi}{N}kn} , k 的取值也是 (0,N−1)(0,N - 1)
所以整体的复杂度是 O(N2)O(N^2)
至少次数是有限次

快速傅里叶变换 FFT

数理原理

首先我们要解释一下上面 X[k]X[k] 的数理意义，nn 相当于是时域的 tt, kk 相当于是 频率 nω0n\omega_0 的 nn, X[k]X[k] 相当于是 频谱函数在 kω0k\omega_0 处的幅值大小
那么接下来我们分析一下 e−j2πNkne^{-j\frac{2\pi}{N}kn} 的等价性，如果可以化简这部分的算法，那么就能降低运算量。
不难发现，当 nn 是偶数的时候，

𝑒^{−𝑗\frac{2𝜋}{𝑁}2𝑛𝑘}=𝑒^{−𝑗\frac{2𝜋𝑛}{𝑁/2}𝑘}=𝑒^{−𝑗\frac{2𝜋𝑛}{𝑁/2}(𝑘+𝑁/2)}

也就是说周期是 N2\frac{N}{2}.
因此，我们可以将 n 分为奇偶数两类分别求和，其中对奇数部分提取一个 e−j2πN/2nke^{-j\frac{2\pi}{N / 2}nk} 的系数，则剩余部分就是偶数部分，同时定义函数 𝑊[𝑛,𝑘]=𝑒^{−𝑗\frac{2𝜋}{𝑁}𝑛𝑘} , 那么我们有

X[k] = ∑_{𝑛=0}^{𝑁/2−1}𝑥[2𝑛]𝑊[2𝑛,𝑘]+𝑒^{−𝑗\frac{2𝜋}{𝑁}𝑘}∑_{𝑛=0}^{𝑁/2−1}𝑥[2𝑛+1]𝑊[2𝑛,𝑘]

我们令偶数部分为 𝐸[𝑘]𝐸[𝑘] ，奇数部分提取公因数后的结果为 𝑂[𝑘]𝑂[𝑘]，则

𝑋[𝑘]=𝐸[𝑘]+𝑊[1,𝑘]𝑂[𝑘]𝑋[𝑘]=𝐸[𝑘]+𝑊[1,𝑘]𝑂[𝑘]

, 其中 𝑥[2𝑛] ,𝑥[2𝑛+1]𝑥[2𝑛] ,𝑥[2𝑛+1] 都是已知且固定的，且 𝑊[2𝑛,𝑘]𝑊[2𝑛,𝑘] 满足 𝑊[2𝑛,𝑘]=𝑊[2𝑛,𝑘+𝑁/2]𝑊[2𝑛,𝑘]=𝑊[2𝑛,𝑘+𝑁/2] 
那么， 𝐸[𝑘]=𝐸[𝑘+𝑁/2]𝐸[𝑘]=𝐸[𝑘+𝑁/2] ， 𝑂[𝑘]=𝑂[𝑘+𝑁/2]𝑂[𝑘]=𝑂[𝑘+𝑁/2]
𝑊[1,𝑘]=−𝑊[1,𝑘+𝑁/2]𝑊[1,𝑘]=−𝑊[1,𝑘+𝑁/2]（推导如下）

𝑒^{−𝑗\frac{2𝜋}{𝑁}(𝑘+𝑁/2)}=𝑐𝑜𝑠(−\frac{2𝜋}{𝑁}(𝑘+𝑁/2))+𝑗\ 𝑠𝑖𝑛(−\frac{2𝜋}{𝑁}(𝑘+𝑁/2)) =𝑐𝑜𝑠(−\frac{2𝜋}{𝑁}𝑘−𝜋)+𝑗\ 𝑠𝑖𝑛(−\frac{2𝜋}𝑁𝑘−𝜋) =−𝑒^{−𝑗\frac{2𝜋}𝑁𝑘}

我们可以得出结论 𝑋[𝑘+𝑁/2]=𝐸[𝑘]−𝑊[1,𝑘]𝑂[𝑘]𝑋[𝑘+𝑁/2]=𝐸[𝑘]−𝑊[1,𝑘]𝑂[𝑘]
注意这里

  • E[k]=Σn=0N/2−1x[2n]W[2n,k]E[k] = \Sigma_{n = 0}^{N / 2 - 1} x[2n]W[2n,k]
  • O[k]=Σn=0N/2−1x[2n+1]W[2n,k]O[k] = \Sigma_{n = 0}^{N / 2 - 1} x[2n + 1]W[2n,k]

蝴蝶结构

其实在上面我们就可以看出，我们会将一个周期为 N 的数据组进行平均分（或者说这就是分治算法的精髓）
然后我们两边齐头并进，如下图
butterfly_structure
我们有 a0=x0+x8a_0 = x_0 + x_8, a1=x0−x8a_1 = x_0 - x_8
这个形式我们称之为 蝴蝶结构

逆快速傅里叶变换 IFFT

Evaluation

我们首先从多项式入手，我们对于一个已知系数的多项式，我们可以得到 其位于某些点的值。但是，如果我们已知 xx 点的值，我们如何快速推算 −x-x 的值？
从奇偶性的角度思考，我们可以将多项式分成 奇数次项 和 偶数次项。再回想上面做的内容，E[k]=E[−k]E[k] = E[-k], O[k]=−O[−k]O[k] = -O[-k]
总的多项式的等式

P(x)=E(x2)+xO(x2)P(x) = E(x^2) + xO(x^2)

我们由此推算 P(−x)=E(x2)−xO(x2)P(-x) = E(x^2) - x O(x^2)
这个操作最牛的地方在于它将一个高次函数降次到 一半的次数，那对于一个人或者电脑来说，次数越低肯定是越好算的。
如果我们想要重复执行一次这个操作，或者说，我们想要递归的执行这个操作直到 base case, 那么我们就需要将 x2nx^{2^n} 进行配对（就是有正有负，奇偶性的使用很依赖于这个条件），意思就是，我们要构造 ±x2n\pm x^{2^n} 都在多项式的 x 坐标上。那么，我们就需要使用 复数的思路。
假设一个复平面中我们有一个单位圆，如果我们要将一个圆 2n2^n 等分, 那么每个节点上面的数就满足 配对要求
所以我们定义变量 ω=e2πin\omega = e^{\frac{2\pi i}{n}}
然后就有求值等式

P(ωj)=E(ω2j)+ωjO(ω2j)P(\omega^j) = E(\omega^{2j}) + \omega^jO(\omega^{2j})

P(ωj+n/2)=E(ω2j)−ωjO(ω2j)P(\omega^{j+n/2}) = E(\omega^{2j}) - \omega^jO(\omega^{2j})

使用 python 表示算法思路，就是

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def FFT(P):
	n = len(P)
	if n == 1:  # base case
		return P
	w = power(e,2 * pi * i / n)
	P_e = [::2]
	P_o = [1::2]
	# recursion for the odd or even part respectively
	E = FFT(P_e)
	O = FFT(P_o)
	y = zeros(n)
	for j in range(n /2):
		y[j] = E[j] + w^j * O[j]
		y[j + n/2] = E[j] - w^j * O[j]
	return y

Interpolation

首先，大量平行的 evaluation 过程是 可以通过矩阵表示的
linear equations
matrix_solver
眼尖的你一定发现了，这就是范德蒙的矩阵，所以求解这类问题的算法可以查阅范德蒙的矩阵的性质
但是，我们也知道，插值法就是求值法的逆运算，也就是
ifft
你还会发现，两个矩阵只有 1/n1/n 以及倒数的关系，也就是说，我们可以通过已有点的值进行反向估计差值的系数集合