Fourier Transform

1. 傅里叶变换

傅里叶变换起源于18世纪研究振动与热传导的问题。欧拉、达朗贝尔和伯努利提出用正弦波描述振动，奠定了思想基础。

19世纪初，法国数学家傅里叶在研究热传导时提出任何函数都可展开为正弦与余弦级数，这就是傅里叶级数，并逐渐发展为傅里叶变换。

1.1 函数的正交

如果是向量的正交，就是向量的内积为$0$。

而函数的正交，我们以区间 $[-\pi,\pi]$ 为例；在 $[0,2\pi]$ 或任意长度为 $2\pi$ 的区间上结论相同

定义定积分：

$$ \langle f,g\rangle:=\int_{-\pi}^{\pi} f(x)g(x) dx $$

表示函数的内积，它的值如果为$0$，我们就称这两个函数在这个周期内正交。

1.2 三角函数正交性

有如下函数

$$ 1, cosx, sinx, cos2x, sin2x, cos3x, sin3x, … $$

其中 $1$ 也可以看成 $cos0x$ 。

首先，我们知道，对于任何非零整数 $k$，都有

$$ \begin{aligned} \int_{-\pi}^{\pi}\cos kxdx & =0\\ \int_{-\pi}^{\pi}\sin kxdx & =0 \end{aligned} $$

因为它们在完整周期上的积分为零。所以 $1$ 与上述其他三角函数都是正交的。

现在讨论三角函数之间的三种情况：

$\cos mx$ 与 $\cos nx$ 的正交性

若 $m\neq n$ 并且 $m\neq 0, n\neq 0$

$$ \begin{aligned} \int_{-\pi}^{\pi}\cos mx\cos nxdx &=\tfrac12\int_{-\pi}^{\pi}\big[\cos(m-n)x+\cos(m+n)x\big]dx \\ &=0 \end{aligned} $$

$\sin mx$ 与 $\sin nx$ 的正交性

若 $m\neq n$ 并且 $m\neq 0, n\neq 0$

$$ \begin{aligned} \int_{-\pi}^{\pi}\sin mx\sin nxdx &=\tfrac12\int_{-\pi}^{\pi}\big[\cos(m-n)x-\cos(m+n)x\big]dx\\ &=0 \end{aligned} $$

$\sin mx$ 与 $\cos nx$ 互相正交

若 $m\neq 0, n\neq 0$

$$ \begin{aligned} \int_{-\pi}^{\pi}\sin mx\cos nxdx &=\tfrac12\int_{-\pi}^{\pi}\big[\sin(m+n)x+\sin(m-n)x\big]dx\\ &=0 \end{aligned} $$

所以上述任意三角函数之间，都是正交的。

现在来计算一下，与自身内积的结果：

$$ \int_{-\pi}^{\pi}\cos^2 nxdx =\tfrac12\int_{-\pi}^{\pi}\big[1+\cos(2n x)\big]dx =\tfrac12\big[2\pi+0\big]=\pi. $$

$$ \int_{-\pi}^{\pi}\sin^2 nxdx =\tfrac12\int_{-\pi}^{\pi}\big[1-\cos(2n x)\big]dx =\tfrac12\big[2\pi-0\big]=\pi. $$

$$ \int_{-\pi}^{\pi} 1\cdot 1,dx=2\pi $$

总结：

任意$m\neq n$的三角函数之间正交。
三角函数在 $[-\pi,\pi]$ 周期的内积是 $\pi$
$1$ 在 $[-\pi,\pi]$ 周期的内积是 $2\pi$

1.3 周期为 $2\pi$ 的傅里叶级数展开

考虑一个周期函数，它的周期为$T = 2\pi$，即满足$f(x)=f(x+2\pi)$。

那么它可以表示为三角级数的和：

$$ \begin{aligned} f(x) & = \sum_{n=0}^{\infin}{a_n \cos{nx}}+\sum_{n=0}^{\infin}{b_n \sin{nx}} \\ & = a_0 + \sum_{n=1}^{\infin}{a_n \cos{nx}}+\sum_{n=1}^{\infin}{b_n \sin{nx}} \end{aligned} $$

如何求出每个 $a_n$ 的值？可以利用正交性来求。

在周期内进行内积，左右同时乘以 $1$ 并积分：

$$ \begin{aligned} \int_{-\pi}^{\pi}{1 \cdot f(x) dx} &= \int_{-\pi}^{\pi}{a_0 dx} + \int_{-\pi}^{\pi}{\sum_{n=1}^{\infin}{a_n \cos{nx} dx}} + \int_{-\pi}^{\pi}{\sum_{n=1}^{\infin}{b_n \sin{nx} dx}} \end{aligned} $$

由于三角函数周期性，后面的三角函数都是积分结果都是 $0$，所以有：

$$ \begin{aligned} \int_{-\pi}^{\pi}{f(x) dx} &= 2 \pi a_0 \\ a_0 &= \frac{1}{2\pi} \int_{-\pi}^{\pi}{f(x) dx} \end{aligned} $$

注意，这里有的教科书上，是令

$$ \begin{aligned} a_0 &= \frac{1}{\pi} \int_{-\pi}^{\pi}{f(x) dx} \end{aligned} $$

然后在公式阶段使用的是：

$$ \begin{aligned} f(x) & = \frac{a_0}{2} + \sum_{n=1}^{\infin}{(a_n \cos{nx}+ b_n \sin{nx})} \end{aligned} $$

它们只是定义习惯，实际上没啥区别。我们后面采用这个定义的方式。

接下来求 $a_n$，依然是同时乘以对应的三角函数然后在周期内进行积分，假设我们要求的是 $n=m$ 时的 $a_m$，依据三角函数正交性：

$$ \begin{aligned} \int_{-\pi}^{\pi}{f(x) \cos{mx} dx} =& \int_{-\pi}^{\pi}{a_0 \cos{mx} dx} \\ &+ \int_{-\pi}^{\pi}{\sum_{n=1}^{\infin}{a_n \cos{nx} \cos{mx} dx}} \\ &+ \int_{-\pi}^{\pi}{\sum_{n=1}^{\infin}{b_n \sin{nx} \cos{mx} dx}} \\ =& 0 + \int_{-\pi}^{\pi}{a_m \cos^2{mx}dx} + 0 \\ =& a_m \pi \end{aligned} $$

所以有：

$$ \begin{aligned} a_n = \frac{1}{\pi}\int_{-\pi}^{\pi}{f(x) \cos{nx} dx} \end{aligned} $$

同理 $b_n$ 也可以求出：

$$ \begin{aligned} b_n = \frac{1}{\pi}\int_{-\pi}^{\pi}{f(x) \sin{nx} dx} \end{aligned} $$

最后，结论如下：

一个周期为 $2\pi$的函数 $f(x)$，可以通过三角函数展开：

$$ \begin{aligned} f(x) & = \frac{a_0}{2} + \sum_{n=1}^{\infin}{a_n \cos{nx}} + \sum_{n=1}^{\infin}{b_n \sin{nx}} \end{aligned} $$

其中：

$$ \begin{aligned} a_0 &= \frac{1}{\pi} \int_{-\pi}^{\pi}{f(x) dx} \\ a_n &= \frac{1}{\pi}\int_{-\pi}^{\pi}{f(x) \cos{nx} dx} \\ b_n &= \frac{1}{\pi}\int_{-\pi}^{\pi}{f(x) \sin{nx} dx} \end{aligned} $$

1.4 周期为 $2L$ 的傅里叶级数展开

上一小节探讨的是在 $2\pi$ 为周期的函数上的三角级数展开，而我们希望可以在任意周期上进行展开。

设原始函数 $f(t)$ 的周期是 $[-L, L]$，如果我们希望一个中间变量 $x$ 的周期是 $[-\pi, \pi]$，那么就可以令：

$$ \begin{aligned} x &= \frac{\pi}{L}t \\ t &= \frac{L}{\pi}x \end{aligned} $$

此时有，

$$ f(t) = f(\frac{L}{\pi}x) \overset{\triangle}{=} g(x) $$

根据在前面的推导中有：

$$ \begin{aligned} g(x) & = \frac{a_0}{2} + \sum_{n=1}^{\infin}{a_n \cos{nx}} + \sum_{n=1}^{\infin}{b_n \sin{nx}}\\ a_0 &= \frac{1}{\pi} \int_{-\pi}^{\pi}{g(x) dx} \\ a_n &= \frac{1}{\pi}\int_{-\pi}^{\pi}{g(x) \cos{nx} dx} \\ b_n &= \frac{1}{\pi}\int_{-\pi}^{\pi}{g(x) \sin{nx} dx} \end{aligned} $$

统统做如下的替换，

$$ \begin{aligned} x &= \frac{\pi}{L}t\\ \Rightarrow dx &= \frac{\pi}{L}dt \end{aligned} $$

以及

$$ nx = n\frac{\pi}{L}t = (\frac{n\pi}{L})t $$

于是有，

$$ \begin{aligned} f(t) & = \frac{a_0}{2} + \sum_{n=1}^{\infin}{a_n \cos{\frac{n\pi}{L}t}} + \sum_{n=1}^{\infin}{b_n \sin{\frac{n\pi}{L}t}}\\ a_0 &= \frac{1}{\pi} \int_{-\pi}^{\pi}{g(x) dx}\\ &= \frac{1}{\pi} \int_{-\pi}^{\pi}{f(\frac{L}{\pi}x) dx}\\ &= \frac{1}{\pi} \int_{-L}^{L}{f(t) \frac{\pi}{L} dt} \\ &= \frac{1}{L} \int_{-L}^{L}{f(t) dt} \end{aligned} $$

同理， $$ \begin{aligned} a_n &= \frac{1}{L}\int_{-L}^{L}{f(t) \cos{\frac{n\pi}{L}t} dt} \\ b_n &= \frac{1}{L}\int_{-L}^{L}{f(t) \sin{\frac{n\pi}{L}t} dt} \end{aligned} $$

在工程领域中，习惯从 $0$ 开始积分，记

$$ T=2L\\ \\ \omega = \frac{\pi}{L} = \frac{2\pi}{T} $$

于是三角级数展开又可以写为：

$$ \begin{aligned} f(t) & = \frac{a_0}{2} + \sum_{n=1}^{\infin}{a_n \cos{n\omega t}} + \sum_{n=1}^{\infin}{b_n \sin{n\omega t}}\\ a_0 &= \frac{2}{T} \int_{0}^{T}{f(t) dt}\\ a_n &= \frac{2}{T}\int_{0}^{T}{f(t) \cos{n\omega t} dt} \\ b_n &= \frac{2}{T}\int_{0}^{T}{f(t) \sin{n\omega t} dt} \end{aligned} $$

1.5 复数形式的傅里叶级数展开

根据欧拉公式

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

得到

$$ \begin{aligned} \cos\theta&=\frac{1}{2}(e^{i\theta}+e^{-i\theta})\\ \sin\theta&=-\frac{1}{2}i(e^{i\theta}-e^{-i\theta}) \end{aligned} $$

代入到公式中得到：

$$ \begin{aligned} f(t) & = \frac{a_0}{2} + \sum_{n=1}^{\infin}{a_n \cos{n\omega t}} + \sum_{n=1}^{\infin}{b_n \sin{n\omega t}}\\ &=\frac{a_0}{2} + \sum_{n=1}^{\infin}{[a_n\frac{1}{2}(e^{in\omega t}+e^{-in\omega t}) - \frac{1}{2}ib_n(e^{in\omega t}-e^{-in\omega t})]}\\ &=\frac{a_0}{2} + \sum_{n=1}^{\infin}{\frac{a_n-ib_n}{2}e^{in\omega t}} + \sum_{n=1}^{\infin}{\frac{a_n+ib_n}{2}e^{-in\omega t}}\\ &=\frac{a_0}{2}e^{0}+\sum_{n=1}^{\infin}{\frac{a_n-ib_n}{2}e^{in\omega t}} + \sum_{n=-\infin}^{-1}{\frac{a_{-n}+ib_{-n}}{2}e^{in\omega t}}\\ &=\sum_{-\infin}^{\infin}{c_n e^{in\omega t}} \end{aligned} $$

其中 $c_n$ 定义为：

$$ c_n = \begin{cases} \frac{a_0}{2}, & n = 0\\ \frac{a_n-ib_n}{2}, & n = 1,2,3,4,…\\ \frac{a_{-n}+ib_{-n}}{2}, & n = -1,-2,-3,-4,… \end{cases} $$

现在具体计算一下 $c_n$，

$$ \begin{aligned} \underset{n=1,2,3,…}{c_n} &= \frac{1}{2}[\frac{2}{T}\int_{0}^{T}{f(t) \cos{n\omega t} dt} -i\frac{2}{T}\int_{0}^{T}{f(t) \sin{n\omega t} dt}]\\ &= \frac{1}{T}\int_{0}^{T}{f(t) (\cos{n\omega t}-i\sin{n\omega t}) dt}\\ &= \frac{1}{T}\int_{0}^{T}{f(t) (\cos{(-n\omega t)} + i\sin{(-n\omega t)}) dt}\\ &= \frac{1}{T}\int_{0}^{T}{f(t) e^{-in\omega t} dt} \end{aligned} $$

$$ \begin{aligned} \underset{n=-1,-2,-3,…}{c_n} &= \frac{1}{2}[\frac{2}{T}\int_{0}^{T}{f(t) \cos{(-n\omega t)} dt} +i\frac{2}{T}\int_{0}^{T}{f(t) \sin{(-n\omega t)} dt}]\\ &= \frac{1}{T}\int_{0}^{T}{f(t) (\cos{(-n\omega t)} + i\sin{(-n\omega t)}) dt}\\ &= \frac{1}{T}\int_{0}^{T}{f(t) e^{-in\omega t} dt} \end{aligned} $$

$$ \begin{aligned} \underset{n=0}{c_n} &= \frac{1}{2}[\frac{2}{T} \int_{0}^{T}{f(t) dt}]\\ &= \frac{1}{T} \int_{0}^{T}{f(t) dt}\\ &= \frac{1}{T} \int_{0}^{T}{f(t) e^{-i0\omega t} dt}\\ &= \frac{1}{T} \int_{0}^{T}{f(t) e^{-in\omega t} dt} \end{aligned} $$

所以 $c_n$ 都可以写为一样的形式：

$$ c_n = \frac{1}{T} \int_{0}^{T}{f(t) e^{-in\omega t} dt} $$

总结：

复数形式的傅里叶变换中，对于周期函数 $f(t) = f(t+T)$，有

$$ \begin{aligned} f(t) &= \sum_{-\infin}^{\infin}{c_n e^{in\omega t}} \\ c_n &= \frac{1}{T} \int_{0}^{T}{f(t) e^{-in\omega t} dt} \end{aligned} $$

1.6 周期为 $\infin$ 的傅里叶级数展开

周期为 $T$ 的函数

$$ f_T(t) = f_T(t+T) $$

可以被解析为下列表示

这里用 $\omega_0$ 表示基频率

$$ \begin{aligned} f_T(t) &= \sum_{-\infin}^{\infin}{c_n e^{in\omega_0 t}} \\ c_n &= \frac{1}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}}{f_T(t) e^{-in\omega_0 t} dt}\\ \omega_0 &= \frac{2\pi}{T} \end{aligned} $$

现在的问题是，如果一个非周期函数，那么该怎么处理？

不是周期函数，换个角度来说，就是周期无穷大，就是无穷大之后才会重复。

我们直接把 $c_n$ 代入回去，同时令 $f(t) = \underset{T \to \infin}{f_T}$，

$$ \begin{aligned} f(t) &= \sum_{-\infin}^{\infin}{\frac{1}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}}{f(t) e^{-in\omega_0 t} dt} e^{in\omega_0 t}}\\ &= \sum_{-\infin}^{\infin}{\frac{\omega_0}{2\pi} \int_{-\infin}^{\infin}{f(t) e^{-in\omega_0 t} dt} e^{in\omega_0 t}} \end{aligned} $$

令 $\omega_0 = d\omega$，那么 $n\omega_0 = \omega$，所以

$$ \begin{aligned} f(t) &= \sum_{-\infin}^{\infin}{\frac{\omega_0}{2\pi} \int_{-\infin}^{\infin}{f(t) e^{-in\omega_0 t} dt} e^{in\omega_0 t}}\\ &= \int_{-\infin}^{\infin}{\frac{1}{2\pi}\int_{-\infin}^{\infin}{f(t) e^{-iomega t} dt} e^{i\omega t} d\omega}\\ &= \int_{-\infin}^{\infin}{\frac{1}{2\pi} \left( \int_{-\infin}^{\infin}{f(t) e^{-i\omega t} dt} \right) e^{i\omega t} d\omega}\\ &= \frac{1}{2\pi} \int_{-\infin}^{\infin}{\left( \int_{-\infin}^{\infin}{f(t) e^{-i\omega t} dt} \right) e^{i\omega t} d\omega} \end{aligned} $$

这里定义 $F(\omega)$，就得到傅里叶变换的最终形式

$$ \begin{aligned} F(\omega) &= \int_{-\infin}^{\infin}{f(t) e^{-i\omega t} dt} \end{aligned} $$

而逆变换是

$$ \begin{aligned} f(t) &= \frac{1}{2\pi} \int_{-\infin}^{\infin}{ F(\omega) e^{i\omega t} d\omega} \end{aligned} $$

1.7 对称形式、工程学派、频率形式

在上面的傅里叶变换中，我们把 $2\pi$ 放在逆变换中，实际上还有下列的形式：

对称形式（常用在数学里），把 $2\pi$ 均分到正变换和逆变换：

$$ F(\omega) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} f(t)e^{-i\omega t} dt $$

$$ f(t) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} F(\omega)e^{i\omega t} d\omega $$

工程学派（常见于电气工程/信号处理），把 $2\pi$ 放在正变换中，逆变换没有 $2\pi$：

$$ F(\omega) = \frac{1}{2\pi}\int_{-\infty}^{\infty} f(t)e^{-i\omega t} dt $$

$$ f(t) = \int_{-\infty}^{\infty} F(\omega)e^{i\omega t} d\omega $$

周期形式（用频率 $f$ 而不是角频率 $\omega$），如果令 $\omega = 2\pi f$，就可以写成：

$$ F(f) = \int_{-\infty}^{\infty} f(t)e^{-i 2\pi f t} dt $$

$$ f(t) = \int_{-\infty}^{\infty} F(f)e^{i 2\pi f t} df $$

2. 离散傅里叶变换

在前面的讨论中，我们涉及到是连续函数和无穷频率。

虽然前面已经学习了傅里叶变换的积分形式，但是实际工程中的最大问题是工程中的采样点不是连续函数。通常是离散的点。

2.1 离散傅立叶变换 DFT

对于一个离散的采样点，比如 N 个点的序列 $x[n]$，直接认为它的周期就是 N。

已知周期为 T 的傅里叶变换形式是：

$$ \begin{aligned} c_n &= \frac{1}{T}\int_{0}^{T}{f(t) e^{-in\omega t} dt}\\ f(t) &= \sum_{n=-\infin}^{\infin}{c_n e^{in\omega t}} \end{aligned} $$

现在令 $\omega = \frac{2\pi}{N}$，$T \to N$，同时把 $\frac{1}{T} \to \frac{1}{N}$ 这个常量移动到逆变换里面。同时原本的 n 用 k 表示，原本的 t 用序号 n 表示。

$$ \begin{aligned} X[k] &=\sum_{n=0}^{N-1}{x[n] \cdot e^{-j\frac{2\pi}{N}kn}}\\ x[n] &= \frac{1}{N} \sum_{k=0}^{N-1}{X[k]\cdot e^{j\frac{2\pi}{N}kn}} \end{aligned} $$

其中 $k=0,1,2,…,N-1$

2.2 证明 DFT 成立

把 $X[k]$ 代入逆变换公式：

$$ \begin{aligned} x[n] &= \frac{1}{N}\sum_{k=0}^{N-1} X[k],e^{j\frac{2\pi}{N}kn} \\ &= \frac{1}{N}\sum_{k=0}^{N-1}\left(\sum_{m=0}^{N-1} x[m],e^{-j\frac{2\pi}{N}km}\right)e^{j\frac{2\pi}{N}kn} \\ &= \frac{1}{N}\sum_{m=0}^{N-1} x[m]\sum_{k=0}^{N-1} e^{j\frac{2\pi}{N}k(n-m)} \end{aligned} $$

考虑内层和

$$ S_{n-m}=\sum_{k=0}^{N-1} e^{j\frac{2\pi}{N}k(n-m)} $$

令 $r=n-m$（整数）。若 $r\equiv 0\pmod N$（即 $n=m$），则每项都是 $1$，于是

$$ S_{0}=N $$

若 $r\not\equiv 0\pmod N$，则这是一个首项为 $1$、公比 $q=e^{j\frac{2\pi}{N}r}\neq1$ 的几何级数，

$$ S_{r}=\frac{1-q^{N}}{1-q}=\frac{1-e^{j2\pi r}}{1-q}=0 $$

因为 $e^{j2\pi r}=1$。因此

$$ S_{n-m}= \begin{cases} N,& n=m,\\ 0,& n\ne m. \end{cases} $$

把它代回去：

$$ x[n]=\frac{1}{N}\sum_{m=0}^{N-1} x[m]S_{n-m}=\frac{1}{N}\big(x[n]\cdot N\big)=x[n]. $$

这就证明了逆变换确实把由前向变换得到的 $X[k]$ 恢复回原序列 $x[n]$。换句话说，复指数在这组离散点上是正交的，且

$$ \frac{1}{N}\sum_{k=0}^{N-1} e^{j\frac{2\pi}{N}k(n-m)}=\delta_{n,m} $$

2.3 共轭对称性

$$ \begin{aligned} X[k] &=\sum_{n=0}^{N-1}{x[n] \cdot e^{-j\frac{2\pi}{N}kn}} \end{aligned} $$

$X[k]$ 的共轭值为：

$$ \begin{aligned} X^{\ast}[k] &=\left( \sum_{n=0}^{N-1}{x[n] \cdot e^{-j\frac{2\pi}{N}kn}} \right)^{\ast} \end{aligned} $$

由于 $x[n]$ 是实数，所以

$$ \begin{aligned} X^{\ast}[k] &=\sum_{n=0}^{N-1}{x[n] \cdot (e^{-j\frac{2\pi}{N}kn}})^{\ast}\\ &=\sum_{n=0}^{N-1}{x[n] \cdot (e^{j\frac{2\pi}{N}kn}})\\ &=\sum_{n=0}^{N-1}{x[n] \cdot (e^{j\frac{2\pi}{N}(kn-Nn)}})\\ &=\sum_{n=0}^{N-1}{x[n] \cdot (e^{-j\frac{2\pi}{N}(N-k)n}})\\ &=X[N-k] \end{aligned} $$

结论：

如果 $x[n]$ 是实数序列，则有

$$ X[k] = X^*[N-k], \quad k=1,2,\dots,N-1 $$

也就是说，DFT 的频谱是共轭对称的。

还有一个有意思的结论，那就是当 N 是偶数的时候， $k = 0$ 和 $k = N/2$ 都只有实数部分的。

$$ \begin{aligned} X[N/2] &= X^{\ast}[N-N/2] = N^{\ast}[N/2]\\ &= \Re{(X[N/2])} \end{aligned} $$

2.4 奈奎斯特采样定理

对长度为 $N$ 的序列 $x[n]$，其 DFT 是

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N}kn}, \quad k = 0,1,\dots,N-1 $$

这里 $k$ 对应的是频率索引，实际频率为

$$ f_k = \frac{k}{N} f_s $$

其中 $f_s$ 是采样率。

根据奈奎斯特采样定理：

连续信号在采样率 $f_s$ 下，最大能无失真恢复的频率是 $f_s/2$（奈奎斯特频率）。
也就是说，$0 \leq f \leq f_s/2$ 才是“有效”区间。

所以，DFT 给出的频率范围 $0 \leq f < f_s$，其中一半其实是镜像频率。

这个结论实际上和前面的共轭对称性也是符合的，有一半的频率实际上是冗余的。

同时从信息的角度来说，$N$ 个实数可以通过 $\frac{N}{2}$ 个复数来恢复信息，所以也是合理的。

2.5 周期性与负频率

离散傅里叶变换 (DFT) 在频率域上具有 周期性。其定义为

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-j\frac{2\pi}{N}kn}, \quad k=0,1,\dots,N-1 $$

如果将频率索引扩展到任意整数 $m$，则有

$$ X[k+mN] = \sum_{n=0}^{N-1} x[n] e^{-j\frac{2\pi}{N}(k+mN)n} = \sum_{n=0}^{N-1} x[n] e^{-j\frac{2\pi}{N}kn} e^{-j2\pi mn} $$

由于 $e^{-j2\pi mn} = 1$，因此

$$ X[k+mN] = X[k] $$

这表明 DFT 在频率索引上以 $N$ 为周期。换句话说，频谱在频率轴上是环形的，超过范围会自动折返。

虽然 DFT 的索引 $k$ 只定义在 $[0, N-1]$，但通过周期性，可以自然引入“负频率”的概念。具体地：

$$ X[-k] = X[N-k]\\ k=-1,-2,-3,… $$

结合前面的共轭对称性，可以发现一个有趣的现象：

$$ X[-k] = X[N-k] = X^{\ast}[k] $$

即正负频率也呈现共轭对称性。

因此，原本的索引

$$ k = 0, 1, 2, \dots, \tfrac{N}{2}-1, \tfrac{N}{2}, \tfrac{N}{2}+1, \dots, N-1 $$

可以等价地写作

$$ k = 0, 1, 2, \dots, \tfrac{N}{2}-1, \tfrac{N}{2}, -\tfrac{N}{2}+1, \dots, -2, -1 $$

3. 快速傅里叶变换

快速傅里叶变换（FFT）的思想最早可追溯到19世纪初。1805年，高斯在研究天体轨道时，已经提出过一种相当于 FFT 的高效算法，但他的工作仅以笔记形式保存，并未被广泛传播。之后一个多世纪里，离散傅里叶变换（DFT）仍以 $O(N^2)$ 的复杂度被直接计算。

直到1965年，IBM 的 James Cooley 与 Princeton 的 John Tukey 在研究数值天气预报时重新发现并系统化了这一方法，他们意识到可以通过分治与旋转因子对称性，大幅减少计算量。论文发表后引起轰动，使大规模频域分析首次成为现实。

有趣的是，FFT 并非“新算法”，而是对古老数学思想的复兴与推广。正是计算机硬件的兴起与实际需求的迫切，才让这项算法走出历史笔记，成为信号处理、通信、图像分析与科学计算中不可或缺的核心工具。

3.1 离散傅里叶变换的复杂度

依据前面，离散傅里叶变换的定义为：

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} kn}, \quad k=0,1,\dots,N-1 $$

直接计算每个 $X[k]$ 需要 $N$ 次乘法以及 $N$ 项求和，总共有 $N$ 个 $k$，因此时间复杂度为：

$$ \mathcal{O}(N^2) $$

当 $N$ 较大时，这样的复杂度难以接受。

快速傅里叶变换（FFT）通过利用对称性和递归分治思想，将复杂度降到：

$$ \mathcal{O}(N \log N) $$

这也是其在信号处理、图像处理、数值计算等领域广泛应用的原因。

3.2 旋转因子的性质

在 DFT 中，我们定义指数部分：

$$ W_N = e^{-j \frac{2\pi}{N}} $$

它被称为 旋转因子 (Twiddle Factor)，于是

$$ W_N^{kn} = e^{-j \frac{2\pi}{N}kn} $$

此时 DFT 公式可写为：

$$ X[k] = \sum_{n=0}^{N-1} x[n] W_N^{kn} $$

旋转因子具有一些重要性质：

周期性

$$ W_N^{a+N} = e^{-j \frac{2\pi}{N}(a+N)} = e^{-j \frac{2\pi}{N}(a)} e^{-j2\pi} = e^{-j \frac{2\pi}{N}(a)} = W_N^a $$

对称性

$$ W_N^{a+\frac{N}{2}} = e^{-j \frac{2\pi}{N}(a+\frac{N}{2})} = e^{-j \frac{2\pi}{N}(a)} e^{-j\pi} = - e^{-j \frac{2\pi}{N}(a)} = - W_N^a $$

缩放性

$$ W_N^{a} = e^{-j \frac{2\pi}{N}(a)} = e^{-j \frac{2\pi}{N/m}(a/m)} = W_{N/n}^{a/m} $$

这些性质可以使得计算量大大减少。

3.3 快速傅里叶变换推导式

下面我们就来看看旋转因子的这些性质是如何被使用的。

设原始序列为 $x[n]$，长度为 $N$（假设 $N=2^m$ 为 2 的幂）。如果实际工程中遇到不满足 2 的幂的点，可以直接补零到 $N=2^m$ 长度，复原的时候截断后面的零即可。

我们还是先回到原来的公式：

$$ X[k] = \sum_{n=0}^{N-1} x[n] W_N^{kn} $$

拆分成奇数部分和偶数部分

$$ \begin{aligned} X[k] &= \sum_{m=0}^{N/2-1} x[2m] \cdot W_N^{k \cdot 2m} + \sum_{m=0}^{N/2-1} x[2m+1] \cdot W_N^{k \cdot (2m+1)} \\ & k=1,2,3,…,N-1 \end{aligned} $$

这里定义一个专门用来取奇偶的函数，

$$ \begin{aligned} x_{even}(m)&=x(2m) & ,m=0,1,2,…,N/2 - 1\\ x_{odd}(m)&=x(2m+1) & ,m=0,1,2,…,N/2 - 1 \end{aligned} $$

同时提取出奇数部分的一个公共旋转因子：

$$ \begin{aligned} X[k] &= \sum_{m=0}^{N/2-1} x_{even}(m) \cdot W_N^{k \cdot 2m} + W_N^k \sum_{m=0}^{N/2-1} x_{odd}(m) \cdot W_N^{k \cdot (2m)}\\ &=\sum_{m=0}^{N/2-1} x_{even}(m) \cdot W_{N/2}^{k \cdot m} + W_N^k \sum_{m=0}^{N/2-1} x_{odd}(m) \cdot W_{N/2}^{km}\\ & k=1,2,3,…,N-1 \end{aligned} $$

定义

$$ \begin{aligned} X_{even}[k] &= \sum_{m=0}^{N/2-1} x_{even}(m) \cdot W_{N/2}^{k \cdot m} \\ X_{odd}[k] &= \sum_{m=0}^{N/2-1} x_{odd}(m) \cdot W_{N/2}^{km} \end{aligned} $$

那么

$$ \begin{aligned} X[k] = X_{even}[k] + W_N^k X_{odd}[k] & , & k=0,1,2,…,N-1\\ \end{aligned} $$

现在 $k$ 比 $m$ 要多一倍的取值范围，如果把 $k$ 也取一半，

$$ \begin{aligned} X[k] &= X_{even}[k] + W_N^k X_{odd}[k] & & k=0,1,2,…,N/2-1\\ X[k+N/2] &= X_{even}[k+N/2] + W_N^{k+N/2} X_{odd}[k+N/2] & & k=0,1,2,…,N/2-1\\ \end{aligned} $$

依据旋转因子的周期性和对称性（这里可以展开代入一下），后半部分又可以写为

$$ \begin{aligned} X[k+N/2] &= X_{even}[k] - W_N^{k} X_{odd}[k] & & k=0,1,2,…,N/2-1 \end{aligned} $$

所以原始序列可以被分解为 2 个部分，前 $N/2$ 部分是偶数部分加奇数部分（乘以一个旋转因子），后 $N/2$ 部分是偶数部分减去奇数部分（乘以一个旋转因子）。

还记得前面小节提到的共轭对称性吗？也是只需要算一半频率就可以求出另一半频率，这里实际上也是暗含了这种性质。

但是这里我们还可以继续进行分解，比如偶数部分又可以拆分为两个部分的和

$$ \begin{aligned} X_{even}[k] &= \sum_{m=0}^{N/2-1} x_{even}(m) \cdot W_{N/2}^{k \cdot m} & & k=0,1,2,…,N/2-1\\ X_{even}[k] &= \sum_{m=0}^{N/4-1} x_{even-even}(m) \cdot W_{N/4}^{k \cdot m} \\ &+ W_{N/2}^{k} \sum_{m=0}^{N/4-1} x_{even-odd}(m) \cdot W_{N/4}^{k \cdot m} & & k=0,1,2,…,N/2-1 \end{aligned} $$

以此类推：

$$ \begin{aligned} X_{even}[k] &= X_{even-even}(k) + W_{N/2}^{k} X_{even-odd}(k) & & k=0,1,2,…,N/4-1\\ X_{even-even}[k] &= X_{even-even-odd}(k) + W_{N/4}^{k} X_{even-even-odd}(k) & & k=0,1,2,…,N/8-1 \end{aligned} $$

可以不断进行递归，最终在 2 个值的时候，只有 1 个偶数和 1 个奇数。

最后在只剩 1 个值的这一层的时候，直接返回值本身即可，因为此时没有奇数序列：

$$ X_{\text{final-even-index}} = x_{\text{final-even-index}} $$

最终的结果就符合上面的推导。

实际上我们算一下 2 点的 DFT 就知道，就是一个加法和一个减法。

比如 $[x0, x1]$ 两个点进行 DFT，

$$ \begin{aligned} X_0 &= x_0 \cdot W_2^{0\cdot0} + x_1 \cdot W_2^{0\cdot1} = x_0 \cdot 1 + x_1 \cdot 1 &= x_0 + x_1\\ X_1 &= x_0 \cdot W_2^{1\cdot0} + x_1 \cdot W_2^{1\cdot1} = x_0 \cdot 1 + x_1 \cdot (-1) &= x_0 - x_1 \end{aligned} $$

3.4 位反转重排

对于 8 个数的 FFT，按照下标序号分组的话：

第一层：

0 2 4 6 / 1 3 5 7

第二层：

0 4 / 2 6 / 1 5 / 3 7

第三层：

0 / 4 / 2 / 6 / 1 / 5 / 3 / 7

8个数可以用 3 个 bit 来表示，仔细观察会发现，第一层是通过最低位的奇偶性来分层的，相同的分在一组；第二层是第二位的奇偶性来分组；最后一层是第高位的奇偶性来分组。

如果我们把位反转：

0 -> 000 -> 000 -> 0
1 -> 001 -> 100 -> 4
2 -> 010 -> 010 -> 2
3 -> 011 -> 110 -> 6
4 -> 100 -> 001 -> 1
5 -> 101 -> 101 -> 5
6 -> 110 -> 011 -> 3
7 -> 111 -> 111 -> 7

通过重排 bit ，获得在内存中的连续块，正好满足迭代 FFT 的计算顺序。

不过为了便于理解，我们后续的代码编写采用递归方式编写，不会进行该操作。

3.5 算法复杂度分析

在每一层合并时，需要做 $N$ 次额外的“蝶形运算”（复加法 + 旋转因子乘法）。

每一层根据奇偶性细分一半，所以总共有 $\log_{2}{N}$ 层。

综上，复杂度为：

$$ \mathcal{O}(N \log N) $$

3.6 IFFT

$$ x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2\pi}{N} kn}, \quad n=0,1,\dots,N-1 $$

可以表示为：

$$ \begin{aligned} x[n] &= \frac{1}{N} (\sum_{k=0}^{N-1} (X[k])^{\ast} (e^{j \frac{2\pi}{N} kn})^{\ast})^{\ast}\\ &= \frac{1}{N} (\sum_{k=0}^{N-1} (X[k])^{\ast} e^{ -j \frac{2\pi}{N} kn})^{\ast}\\ &= \frac{1}{N} (\text{FFT}((X[k])^{\ast}))^{\ast} \end{aligned} $$

如果用 Python 可以表示为：

x = np.conj(np.fft.fft(np.conj(X))) / N

这里注意的一点是，如果输入序列是实数，那么 IFFT 返回只有实部，所以虚数可以不取。

3.7 FFT Python Code

假定输入序列是 $N = 2^{m}$。

有如下代码：

import numpy as np

def fft(x):
    N = len(x)
    if N <= 1:
        return x

    X_even = fft(x[0::2])
    X_odd = fft(x[1::2])

    X = [0] * N
    for k in range(N//2):
        t = np.exp(-2j * np.pi * k / N) * X_odd[k]
        X[k] = X_even[k] + t
        X[k + N//2] = X_even[k] - t
    return X


def ifft(X):
    N = len(X)
    x_time = np.conj(fft(np.conj(X))) / N
    return x_time


if __name__ == "__main__":

    x = [1,t(x)
    print("FFT:")
    for k, v in enumerate(X):
        print(f"X[{k:2d}] = {v.real:8.2f} + {v.imag:8.2f}j")

    # IFFT
    x_rec = ifft(X)
    print("\nIFFT:")
    for n, v in enumerate(x_rec):
        print(f"x[{n:2d}] = {v.real:8.2f} + {v.imag:8.2f}j")   2, 3, 4, 5, 6, 7, 8]
    print("x[n]: ", x, "\n")

    # FFT
    X = ff

结果为：

x[n]:  [1, 2, 3, 4, 5, 6, 7, 8] 

FFT:
X[ 0] =    36.00 +     0.00j
X[ 1] =    -4.00 +     9.66j
X[ 2] =    -4.00 +     4.00j
X[ 3] =    -4.00 +     1.66j
X[ 4] =    -4.00 +     0.00j
X[ 5] =    -4.00 +    -1.66j
X[ 6] =    -4.00 +    -4.00j
X[ 7] =    -4.00 +    -9.66j

IFFT:
x[ 0] =     1.00 +    -0.00j
x[ 1] =     2.00 +     0.00j
x[ 2] =     3.00 +    -0.00j
x[ 3] =     4.00 +     0.00j
x[ 4] =     5.00 +    -0.00j
x[ 5] =     6.00 +    -0.00j
x[ 6] =     7.00 +     0.00j
x[ 7] =     8.00 +    -0.00j

3.8 2D FFT IFFT

2D DFT 的定义

设输入矩阵 $f[m,n]$ 大小为 $M \times N$，2D DFT 定义为：

$$ F[k,l] = \sum_{m=0}^{M-1}\sum_{n=0}^{N-1} f[m,n] , e^{-j 2\pi \left(\frac{k m}{M} + \frac{l n}{N}\right)} $$

注意指数是行方向 + 列方向的和。

因为指数是可分解的：

$$ e^{-j 2\pi (\frac{k m}{M} + \frac{l n}{N})} = e^{-j 2\pi \frac{k m}{M}} \cdot e^{-j 2\pi \frac{l n}{N}} $$

代入原式：

$$ F[k,l] = \sum_{m=0}^{M-1}\sum_{n=0}^{N-1} f[m,n] , e^{-j 2\pi \frac{k m}{M}} , e^{-j 2\pi \frac{l n}{N}} $$

这就意味着可以先对 n 做求和，再对 m 做求和：

$$ F[k,l] = \sum_{m=0}^{M-1} \left[ e^{-j 2\pi \frac{k m}{M}} \left( \sum_{n=0}^{N-1} f[m,n] , e^{-j 2\pi \frac{l n}{N}} \right) \right] $$

这也就是说，可以拆分为先对行做 FFT，得到结果后再在列方向做一次 FFT。（也可以先列后行反过来，不影响结果）

算法时间复杂度上，DFT 的时间复杂度为：

$$ \mathcal{O}(M^2 N^2) $$

拆分为 2 个 1D FFT 后为：

$$ \mathcal{O}(M N \log N + N M \log M) = \mathcal{O}(M N (\log M + \log N)) $$

这里给出 2D 的 FFT 代码：

import numpy as np
from fft1d import fft, ifft

def fft2d(matrix):
    temp = np.array([fft(row) for row in matrix])
    temp = np.array([fft(col) for col in temp.T]).T
    return temp

def ifft2d(matrix):
    temp = np.array([ifft(row) for row in matrix])
    temp = np.array([ifft(col) for col in temp.T]).T
    return temp

if __name__ == "__main__":
    x = np.array([
        [1, 2, 3, 4],
        [5, 6, 7, 8],
        [9,10,11,12],
        [13,14,15,16]
    ], dtype=float)

    print("x[i,j]:")
    print(x)

    # 2D FFT
    X = fft2d(x)
    print("\n2D FFT X:")
    print(np.round(X, 2))

    # 2D IFFT
    x_rec = ifft2d(X)
    print("\n2D IFFT x_rec:")
    print(np.round(x_rec.real, 2))

结果如下：

x[i,j]:
[[ 1.  2.  3.  4.]
 [ 5.  6.  7.  8.]
 [ 9. 10. 11. 12.]
 [13. 14. 15. 16.]]

2D FFT X:
[[136. +0.j  -8. +8.j  -8. +0.j  -8. -8.j]
 [-32.+32.j   0. +0.j   0. +0.j   0. +0.j]
 [-32. +0.j   0. +0.j   0. +0.j   0. +0.j]
 [-32.-32.j   0. +0.j   0. +0.j   0. +0.j]]

2D IFFT x_rec:
[[ 1.  2.  3.  4.]
 [ 5.  6.  7.  8.]
 [ 9. 10. 11. 12.]
 [13. 14. 15. 16.]]

4. 参考致谢

花了几天时间终于写完了，感觉得到了进化。

个人感觉大量的科普视频和参考文章都局限于过多的打比方和举例子。实际上从第一性原理直接彻底了解反而是最简单的。

按照杨振宁的话说：“当直觉与书本知识有冲突，是最好的学习机会，必须抓住，把本来的直觉错误想清楚，形成新的直觉”。

特别感谢下列视频作者：

DR_CAN¹的傅里叶变换推导过程。
nksunmoon²的 FFT 过程讲解。

DR_CAN , 2018, 纯干货数学推导_傅里叶级数与傅里叶变换, https://www.bilibili.com/video/BV1Et411R78v ↩︎
nksunmoon, 2022, 快速傅里叶变换原理, https://www.bilibili.com/video/BV1Ng411v7Dk ↩︎

1. 傅里叶变换#

1.1 函数的正交#

1.2 三角函数正交性#

1.3 周期为 $2\pi$ 的傅里叶级数展开#

1.4 周期为 $2L$ 的傅里叶级数展开#

1.5 复数形式的傅里叶级数展开#

1.6 周期为 $\infin$ 的傅里叶级数展开#

1.7 对称形式、工程学派、频率形式#

2. 离散傅里叶变换#

2.1 离散傅立叶变换 DFT#

2.2 证明 DFT 成立#

2.3 共轭对称性#

2.4 奈奎斯特采样定理#

2.5 周期性与负频率#

3. 快速傅里叶变换#

3.1 离散傅里叶变换的复杂度#

3.2 旋转因子的性质#

3.3 快速傅里叶变换推导式#

3.4 位反转重排#

3.5 算法复杂度分析#

3.6 IFFT#

3.7 FFT Python Code#

3.8 2D FFT IFFT#

4. 参考致谢#