Skip to content

Introduction to FFT

中文互联网中讲解 FFT 博客的作者大致分为两类:一类人大概是信号人,另一类人是 OI (信息竞赛)人。这两种博客往往内容其实没有太大的交集,虽然讲解的是同一个算法。

在我高中时期学习这个算法的时候就感到非常非常的困惑,因为我无法理解为什么上来就要代入一些特殊的点,也无法理解为什么逆变换和正变换这么的类似。更令人迷惑的是当你试图查找一些资料,大部分傅里叶变换的讲解都在讲时域、频域等等没听说过的事情,搞不懂到底这两者之间有什么联系。

也有过很多人问我或是交流这两者之间的关系,为了更好的阐述,我们写一篇博客(其实主要是防止自己再一次忘掉,我已经数不清忘掉多少次了)。

历史

热方程

为了更好的理解,我们先从历史说起。

在傅里叶之前欧拉就已经用三角函数的级数来解决一些弦振动问题了,但是傅里叶第一个提交论文声明:

任何周期函数都可以表示成三角函数的级数

当时大部人并不接受这种说法,因为它并没有很严格的证明。事实上这个结论也确实是错的,但在实际应用中大多数情况都满足 Dirichlet 条件(一个傅里叶变换存在的充分不必要条件),所以它依然是一个很强大实用的工具。

最早这个工具被傅里叶用来解决热方程,一种极其重要的偏微分方程:

\[ \frac{\partial u}{\partial t} = \nabla^2u \]

这个式子中 \(u\) 代表温度,\(t\) 代表时间,\(x\) 代表坐标。式子左边是温度随时间的变化,右侧可以理解为一个点和周围点温度的差的大小。大概意思是说如果你比周围的人热的多,那你温度就会快速的下降。

我们发现如果我们让 \(t=0\) 时初始条件的 \(u\) 为正余弦函数,那么这个方程将非常容易求解;同时这里的算符都满足线性性,所以一些三角函数的叠加也是简单的。

那如果我们能让现实里的 \(u_{t=0}\) 分解成一系列三角函数,事情就简单了。

信号与系统

满足线性性的系统还有很多,除了微分方程的解之外还有信号里的 LTI 系统。

此时傅里叶变换可以被理解为从时域(一个波,横坐标是时间,纵坐标是振幅)到频域(把一个波分解成若干个正余弦波,此时一个波可以表示为横坐标是频率,纵坐标是对应频率分量的强度)的变换。

快速乘法

讲了这么多,到底和快速乘法有什么关系?

后面我们会介绍到原域的卷积(乘积)对应傅里叶变换后的域上的乘积(卷积)。

而我们注意到乘法本身可以看成是各位数字的卷积,所以如果我们有一种可以快速求傅里叶变换的操作,那真是太棒了。我们可以先快速傅里叶变换,再对变换后的函数逐点乘积,最后再变换回去。这样做的好处是直接求卷积复杂度很高,而逐点乘积就简单多了。

快速傅里叶变换

历史上有很多人重新发明过很多次傅里叶变换,大家一般认为最早是高斯发现了这个算法,但通用的算法是 James Cooley 和 John Tukey 发明的,他们分析了算法的复杂度是 \(O(n\log n)\).

但是这篇文章和傅里叶的文章据说因为写的非常晦涩,所以大家一开始也看不懂。所以我们学起来很费劲也是有道理的 :)

缘起 - 无穷维向量空间

在高中数学中,我们应该学习了欧氏空间的向量和内积。实际上内积这个概念可以推广到很多空间,比如这里我们考虑函数空间。

对于复函数空间,我们可以定义内积为:

\[ \langle f, g \rangle := \int_{-\infty}^{+\infty} f(x)\bar{g}(x) \text{d} x \]

这里 \(\bar{\cdot}\) 表示共轭。

实际上函数的内积和有限维的向量内积形式上非常相似,可以把函数想象成一个无穷维的向量,
每个点的函数值就是这个向量在这一维的投影值。

这里一切看起来非常合理,但是我们需要注意我们的函数应当在 \(L^2\) 上,否则内积不保证一定存在。而在 \(L^2\) 里这是 Cauchy Schwarz 保证的,两个平方可积的函数的内积一定是存在的。这在别的空间并不一定成立。比如我们考虑定义在 \([0,1]\) 上的函数 \(f(x) = \frac{1}{\sqrt{x}} \in L^1\)

同时如果我们对例如 \(L^3\) 定义内积,那么内积自然诱导得到的范数和自带的 \(l^3\) norm 并不匹配。

在物理上,一个平方可积的函数往往代表着能量有限,所以在物理上这样的数学假设一般都是很合理的。而且 Parseval's Theorem 告诉我们时域的能量是等于频域的能量的,这种能量守恒在物理上告诉我们这种变换可逆,所以不丢失信息。

Stone Weierstrass Theorem

SMT 这个定理证明了对于一个紧致的Hausdorff空间 \(X\)(不同的点可以找到不相交的邻域)
\(A\)\(C(X,R)\) 的subalgebra,且 \(A\) 包含一个 non-zero constant function,那么 \(A\)\(C(X,R)\) 上稠密当且仅当 \(A\) separates points. separates points 定义为对于任意一个点我都能找到两个函数使得在这个点上的函数值不相等。这里稠密指的是 supremum norm \(||f|| = \sup_{a <= x <= b} |f(x)|\).

这个定理告诉我们对于一个连续的周期函数(周期函数可以看成是一个 \(S^1\) 上的函数)
我们总是能用三角函数的线性组合无限逼它(因为三角函数的线性组合是稠密的)

这个定理很强大,也可以用在多项式函数上等等,但是它无法取代 Dirchlet 判定依据因为
Dirchlet 更一般一些。

函数空间与投影

既然函数可以构成一个函数空间,我们自然地思考这个空间的基是什么?随着各位数学家的推进,
大家发现三角函数或者说复指数函数是非常好的基。对于一个周期函数,我们总是可以用一族(各个倍频)复指数基函数的线性组合逼近这个周期函数。
具体来说对于周期函数 \(f(t)\) with period \(T\), i.e. frequency \(f_0 = \frac{1}{T}, w_0 = \frac{2\pi}{T} = 2\pi f_0\).
We have:

\[ F(n) = \langle f(t), e^{jn2\pi f_0 t} \rangle = \frac{1}{T} \int_T f(t) e^{-jn2\pi f_0 t} dt \]

这里负号是因为我们需要共轭;而前面的系数 \(\frac{1}{T}\) 是因为我们需要归一化。否则这个函数也可以认为是周期为 \(2T\),此时频率的分量居然就因此而变大很明显是不合理的。

我们已经有了各个倍频的投影 \(F(n)\),我们如何复原原本的函数?

\[ f(t) = F(0) + F(1)e^{j2\pi f_0 t} + F(-1)e^{-j2\pi f_0 t} + F(2) \cdots \]

注意到这里的 \(F(0)\) 就是我们的直流分量。

这样的离散的数列就是傅里叶级数。

从有限到无穷

我们知道傅里叶变换可以对不周期的函数进行,比如冲激函数那肯定不是周期函数,这是怎么做到的呢?
一个非周期的函数可以理解为它的周期无穷大。所以整理之后我们可以得到:

\[ F(\omega) = \int_{-\infty}^{\infty} f(t) e^{-j\omega t} dt \]
\[ f(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{j\omega t} d\omega \]

拉普拉斯变换

傅里叶变换毫无疑问有局限性,它无法处理发散的函数。这是因为它只考虑虚轴 \(j\omega\),也就是说只看纯频率。 而拉普拉斯变换则更一般,变换后的域不仅有频率还有衰减或者是增长率。 傅里叶变换只有虚轴而拉普拉斯变换是一个复平面。

有趣的是拉普拉斯早于傅里叶发明了这个数学工具,傅里叶后来发明了物理思想。赫维赛德在工程中为了解决微分方程自创了完全不严谨的数学方法,后来的数学家发现,好吧赫维赛德是对的,他用的东西其实就是拉普拉斯变换。

拉普拉斯当时并没有想到这个东西可能有深刻的物理意义。一开始他只是想推广离散数列的生成函数。生成函数可以让数列的卷积变成简单的生成函数的乘积,非常棒。那么连续的函数呢?

既然离散数列可以对应为一个生成函数 \(\sum a_n x^n\) 那么一个连续的函数应该也能有一个生成函数:\(\int f(t) x^t dt\)。这样概率密度函数的卷积计算也可以变成简单的乘法计算。

上面这个函数的生成函数如果我们带入 \(x = e^{-s}\) 就会发现这就是拉普拉斯变换的式子。

离散傅里叶变换

傅里叶变换的形式非常优美,可惜实际中我们没办法真的处理连续信号,所以我们需要离散傅里叶变换。

首先我们需要把傅里叶变换积分中的连续的时间变成离散的采样,我们每隔时间 \(T\) 采样一下,得到 \(f(0T), f(1T), f(2T), \dots f((N-1)T)\)。 然后我们把积分变成有限的近似求和。

\[ F(\omega) = \sum_{n=0}^{N-1} f(nT) e^{-j\omega n T} T \]

现在变换后的频域依然是连续的,实际中还是没办法处理。我们进一步对变换后的频域进行采样。考虑到我们时域上的采样频率是 \(f_s = 1/T\),我们已经丢失了更高的频率的信息(比如说两倍三倍采样频率的信号实际上我们采不出来)那么变换后得到的频域最多也就是 \([0, 2\pi / T]\)

我们依然把这个频率区间切割成 \(N\) 份,第 \(k\) 个频率是 \(w_k = k \cdot \frac{2\pi}{N T}\)

最终我们就获得了

\[ X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j (\frac{2\pi k}{NT}) (nT)} T \]

回到最初:快速乘法

让我们再回到那个问题,为什么快速傅里叶变换可以用来加速大数乘法?

这个问题的答案分为两个部分。首先是为什么傅里叶变换可以用来做大数乘法,以及为什么傅里叶变换可以变的快速。

我们先来看另一个东西:生成函数。我们考虑幂级数生成函数。给定两个序列 \(a_0, a_1, \cdots\)\(b_0, b_1, \cdots\),我们定义生成函数

\[ A(x) = \sum_{n >= 0} a_n x^n, \;\;\; B(x) = \sum_{n >= 0} b_n x^n \]

他们的离散卷积就是

\[ c_n = \sum_{k = 0}^n a_k b_{n-k} \]

现在看两个生成函数的乘积

\[ A(x)B(x) = \sum_{n >= 0} (\sum_{k=0}^n a_k b_{n-k}) x^n = \sum_{n >= 0} c_n x^n \]

我们发现这两个生成函数的乘积对应的函数的技术就是原来的系数的卷积。也就是说如果我们认为从数列到生成函数是一个变换的话,那么这个变换也满足卷积后变换等于变换后乘积。

既然如此生成函数是不是也可以用来做大数乘法呢?答案是肯定的,我们把各位数字当成是系数,然后把生成函数乘在一起就可以得到乘法的结果,不需要计算各位数字的卷积。然而这并没有意义,因为如果我们没办法求出生成函数的封闭形式,生成函数相乘无非也只能是把系数展开卷积。

所以把卷积变成乘积并不是傅里叶变换的神奇特点,这更像是指数自带的特性。卷积的平移运算的特征函数就是指数,而傅里叶变换相当于把一些特殊的值带入了生成函数,自然也具有卷积乘积对应的特性。

既然生成函数并不能加速大数乘法,为什么傅里叶变换可以呢?我们注意到生成函数求出来后获得封闭形式是很困难,没有系统化方法的。这导致两个生成函数相乘并没有让任务变简单。而离散傅里叶变换最终得到的是频域函数 \(X(\omega)\)\(N\) 个点。换句话说,傅里叶变换得到的结果本身就是点值表示的函数,两个函数的这 \(N\) 个点对应相乘就等于函数相乘后的这 \(N\) 个点。这使得函数相乘变成简单的点对点相乘。可以理解为生成函数的相乘也可以看成是频域函数的相乘,没有明显的区别。

Comments