傅里叶变换：不只是多项式乘法

# 傅里叶变换：不只是多项式乘法 ## 动机 我们的信息科学与技术导论课程期末需要做一个project，笔者选择的题目是跟信号处理有关的，而谈到信号处理就必然离不开傅里叶变换。笔者作为一个前OIer，曾经接触过用于求多项式乘法的FFT、NTT算法，并大致了解过生成函数、多项式相关的理论，同时又对数学分析有着浓厚的兴趣，于是就花了一些时间简单学习了一下傅里叶分析，初步了解了它在信号处理中的应用。之后，我决定将我的笔记和我看到的一些资料整理成这篇随笔。 ## 前言 音频信号是什么？ 如果你用过Cool Edit，Fruity Loops，Cubase之类的音乐制作软件，你应该看到过类似这样的图像：  这是我用MATLAB读入了一个音频文件之后直接调用`plot`函数画出的图像，它的横坐标是时间，纵坐标是信号的强度，也就是说它是关于时间的函数$f(t)$的图像，称为**时域**图像。 现在我们来看一个简单一点的信号：$f(t)=\sin\displaystyle\frac{t}{100}$，它的时域图像如下  利用数学知识容易求得，这个信号的频率是$\frac1{200\pi}$，注意这里的单位并不是赫兹。这里我没有规定任何一个量的单位，这不重要。 问题来了，一个正弦信号的频率容易求出，那如果是复杂的信号呢？例如像第一张图中那样的信号，它肯定不是像正弦那样只包含一个频率的简单信号，那么你能看出它包含了哪些频率，每种频率的强度是多少吗？ 这说明我们需要一种办法，将一个信号的时域函数$f(t)$变换为**频域**函数$F(\omega)$，即我们需要一张频域图像，横轴是不同的频率，纵轴代表信号强度。这里，我对第一张图的那个信号作了相应的变换，得到的频域图像如下：  可以看到，这个信号中掺杂了各种频率的信号，但在$2500Hz$左右的频率的能量尤其突出。 本文将会介绍这种变换的数学原理，并简单地介绍一些应用。由于笔者水平有限，应用部分仅仅给出两个简单的例子，有兴趣的读者可以自行查阅资料进一步研究。 本文需要一些数学分析的基础，包括积分、级数等基础知识。 本文侧重于定理的应用，因而略去了一些证明。 作者水平有限，如有错误欢迎指出，不胜感激。 ## 傅里叶级数 傅里叶分析来自于一个有趣的猜想。数学家、物理学家傅里叶在研究一维杆状物体的热流问题时，求解了一个偏微分方程组，发现其解可以用一个三角级数表示： $$ f(x)=\displaystyle\sum\limits_{n\geqslant 1}b_n\sin\displaystyle\frac{n\pi}lx. $$ 这意味着，对于一个事先给定的定义在$[0,l]$上的函数$f(x)$，不管这个函数的定义如何，甚至可能有间断点，它也许都能表示成无穷多个三角函数的和 $$ f(x)=\displaystyle\frac{a_0}{2}+\displaystyle\sum\limits_{n\geqslant 1}\left(a_n\cos\displaystyle\frac{n\pi}{l}x+b_n\sin\displaystyle\frac{n\pi}{l}x\right), $$ 也就是无穷多个不同频率、不同振幅的正弦、余弦曲线的叠加？如果真的是这样，那么这些系数又如何确定？ 为了正式介绍傅里叶级数，我们首先需要一些基础知识。 ### 周期函数 首先是周期函数$f(x)$，设定义域为$\mathbb R$，周期为$T$，则$f(x)=f(x+T)$对任意$x$都成立。这是周期函数的基本定义，但我们还可以证明一个有趣的性质： $$ \displaystyle\int_a^{a+T}f(x)\mathrm dx=\displaystyle\int_0^Tf(x)\mathrm dx, $$ 即周期函数在一个周期上的积分与起点无关。 **证** 设$k\in\mathbb Z$满足$a\leqslant kT<a+T$，则 $$ \begin{aligned} \displaystyle\int_a^{a+T}f(x)\mathrm dx&=\displaystyle\int_a^{kT}f(x)\mathrm dx+\displaystyle\int_{kT}^{a+T}f(x)\mathrm dx\\ &=\displaystyle\int_a^{kT}f(x)\mathrm dx+\displaystyle\int_{(k-1)T}^af(x)\mathrm dx\\ &=\displaystyle\int_{(k-1)T}^{kT}f(x)\mathrm dx\\ &=\displaystyle\int_0^Tf(x)\mathrm dx. \end{aligned} $$ 证毕。 这个性质在几何上十分显然，读者可以挑一个周期函数画画看。 注意，对于一个周期为$2l$的函数，将它沿$x$轴按一定比例缩放就可以得到一个周期为$2\pi$的函数，所以周期具体是多少并没有本质上的影响。 ### 周期延拓 对于一个定义在$[-l,l]$上的函数$f(x)$，我们可以将它的图像向左、向右复制，将它变成一个定义在$\mathbb R$上的周期为$2l$的函数，这就是**周期延拓**。但是注意，在区间的端点处（即$l$的奇数倍处）规定函数值为$\displaystyle\frac{f(l)+f(-l)}2$。形式化地，定义为 $$ F(x)=\begin{cases}f(x-2nl),&(2n-1)l<x<(2n+1)l,\\\displaystyle\frac{f(l)+f(-l)}2,&x=(2n+1)l\end{cases}. $$ 如果函数$f(x)$不是定义在$[-l,l]$上的，而是$[a,b]$上的，我们可以令$2l=b-a$，并将它平移到$[-l,l]$上，再作周期延拓即可。因此，以下所有讨论虽然针对周期函数，但只要将结果限制在一个给定的区间上，其结果对于定义在这个区间上的函数是成立的，只是在端点处的结果要视具体情况而定。 ### 偶延拓、奇延拓 对于定义在$[0,l]$上的函数，我们还可以进行所谓的偶延拓和奇延拓，就是先将它延拓到$[-l,l]$上的偶函数或奇函数，再作周期延拓，就产生了定义在全体实数$\mathbb R$上的周期函数。 ### 三角函数的正交性 函数列 $$ 1,\cos x,\sin x,\cos 2x,\sin 2x,\cdots,\cos nx,\sin nx,\cdots $$ 是所谓的**三角函数系**，它们具有共同的周期$2\pi$，并且对于任何的$m,n\in\mathbb N$，有如下三个式子成立 $$ \begin{aligned} \displaystyle\frac{1}{\pi}\displaystyle\int_{-\pi}^\pi\cos mx\cos nx\mathrm dx&=(1+[m=0])[m=n],\\ \displaystyle\frac{1}{\pi}\displaystyle\int_{-\pi}^\pi\sin mx\sin nx\mathrm dx&=[m=n],\\ \displaystyle\frac{1}{\pi}\displaystyle\int_{-\pi}^\pi\cos mx\sin nx\mathrm dx&=0. \end{aligned} $$ 证明留给读者。提示：利用积化和差公式。 这三个式子描述了这样一个事实：对于来自三角函数系中的任意两个不同的函数$f(x),g(x)$，在区间$[-\pi,\pi]$上由如下方法定义的**内积** $$ \displaystyle\int_{-\pi}^{\pi}f(x)g(x)\mathrm dx $$ 是$0$。类比向量的内积与正交性的关系，我们说三角函数系是**正交**的。 ### $\mathbb R$上的周期函数的傅里叶级数 下面让我们做一些大胆的尝试。注意，这里的过程**并不严谨**，只是一个猜想。伟大的成果往往来自于大胆的猜想。 假设周期为$2\pi$，定义域为$\mathbb R$的函数$f(x)$真的可以写成如下三角级数 $$ f(x)=\displaystyle\frac{a_0}2+\displaystyle\sum\limits_{n\geqslant 1}\left(a_n\cos nx+b_n\sin nx\right), $$ 这里为了方便，我们将第$0$项单独拿出来。 现在我们试着求出系数$a_n$和$b_n$。我们可以在两边同乘以$\cos mx(m\in\mathbb N^*)$： $$ f(x)\cos mx=\displaystyle\frac{a_0}2\cdot\cos mx+\displaystyle\sum\limits_{n\geqslant 1}a_n\cos nx\cos mx+\displaystyle\sum\limits_{n\geqslant 1}b_n\sin nx\cos mx $$ 接下来，我们试着取$[-\pi,\pi]$上的积分。注意到 $$ \displaystyle\int_{-\pi}^{\pi}\displaystyle\frac{a_0}2\cos mx\mathrm dx=0 $$ 是十分显然的，又由三角函数的正交性， $$ \displaystyle\int_{-\pi}^{\pi}a_n\cos nx\cos mx\mathrm dx=\pi a_n[m=n]=\pi a_m[m=n], $$ $$ \displaystyle\int_{-\pi}^{\pi}b_n\sin nx\cos mx\mathrm dx=0. $$ 所以 $$ \displaystyle\int_{-\pi}^{\pi}f(x)\cos mx\mathrm dx=\pi a_m. $$ 但这里是有问题的：我们不经意间交换了无限求和符号和积分符号的顺序，这不是一定成立的。事实上，只有当这个函数项级数**一致收敛**时，这样的操作才正确。但现在我们为了方便，不妨假设它是一致收敛的，于是我们得到了 $$ a_m=\displaystyle\frac{1}{\pi}\displaystyle\int_{-\pi}^{\pi}f(x)\cos mx\mathrm dx. $$ 同理，有 $$ b_m=\displaystyle\frac{1}{\pi}\displaystyle\int_{-\pi}^{\pi}f(x)\sin mx\mathrm dx. $$ 这里，我们将常数项记为$\displaystyle\frac{a_0}2$而不是$a_0$，是为了让$a_0$的表达式与其它的$a_m$一致。我们称由上式给出的$a_n,b_n$为**傅里叶系数**。为了保证傅里叶系数的存在性，我们再假设$f(x)$在$[-\pi,\pi]$上是**可积**的。在这种情况下，我们就可以计算出傅里叶系数的值，然后构造如下的三角级数 $$ f(x)\sim\displaystyle\frac{a_0}2+\displaystyle\sum\limits_{n\geqslant 1}\left(a_n\cos nx+b_n\sin nx\right), $$ 这称为$f(x)$的**傅里叶级数**。 我们已经作了太多的假设，就这样不明不白地导出了这些式子，难免心里发慌。幸好，**狄利克雷**(Dirichlet)大师给出了如下的定理：设函数$f(x)$以$2\pi$为周期， - 如果在任何有限区间上$f(x)$是**分段可微**的，那么它的傅里叶级数在整个数轴上都**收敛**，并且 $$ \displaystyle\frac{a_0}2+\displaystyle\sum\limits_{n\geqslant 1}\left(a_n\cos nx+b_n\sin nx\right)=\displaystyle\frac{f(x+0)+f(x-0)}{2}, $$ 其中$f(x+0)$和$f(x-0)$分别是函数$f$在$x$点处的右极限和左极限。 - 如果在任何有限区间上$f(x)$是**分段可微**的，并且还处处连续，那么它的傅里叶级数在整个数轴上**一致收敛**于$f(x)$。 本文略去了狄利克雷定理（以及其他关于收敛性的定理）的证明，因为我们更侧重于应用。有兴趣的读者可以自行查阅相关资料了解。 现在让我们看一个例子：设有函数$f(x)$以$2\pi$为周期，在一个周期$[-\pi,\pi]$内的定义如下： $$ f(x)=\begin{cases}x,&0\leqslant x<\pi\\0,&-\pi\leqslant x<0\end{cases}. $$ 注意这个函数在任何有限区间上都分段可微，但它并不是处处连续的，因此它满足狄利克雷定理的第一条，将会收敛于$\displaystyle\frac{f(x+0)+f(x-0)}{2}$。当$x$是$\pi$的奇数倍时，这个结果是$\displaystyle\frac{\pi}2$；否则就是$f(x)$。 按照定义计算傅里叶系数，这里只需要用到非常基本的定积分计算，过程略去。我们得到 $$ a_0=\displaystyle\frac{\pi}2,a_n=\displaystyle\frac{(-1)^n-1}{n^2\pi},b_n=\displaystyle\frac{(-1)^{n+1}}{n}. $$ 于是我们有 $$ \displaystyle\frac{\pi}4+\displaystyle\sum\limits_{n\geqslant 1}\left(\displaystyle\frac{(-1)^n-1}{n^2\pi}\cos nx+\displaystyle\frac{(-1)^{n+1}}{n}\sin nx\right)=\begin{cases}f(x),&x\neq(2k-1)\pi,\\\displaystyle\frac{\pi}2,&x=(2k-1)\pi.\end{cases} $$ 这个式子看起来复杂，却又暗藏玄机。如果我们令$x=0$，则$\sin nx$都是$0$，$\cos nx$都是$1$，又注意到$(-1)^n-1$在$n$为偶数时为零，将式子作适当整理可得 $$ \displaystyle\sum\limits_{n\geqslant 1}\displaystyle\frac{1}{(2n-1)^2}=\displaystyle\frac{\pi^2}8. $$ 进而又可以求出另外几个级数，例如 $$ \displaystyle\sum\limits_{n\geqslant 1}\displaystyle\frac{1}{n^2}=\displaystyle\sum\limits_{n\geqslant 1}\displaystyle\frac{1}{(2n-1)^2}+\displaystyle\sum\limits_{n\geqslant 1}\displaystyle\frac{1}{(2n)^2}=\displaystyle\frac{\pi^2}8+\displaystyle\frac14\displaystyle\sum\limits_{n\geqslant 1}\displaystyle\frac{1}{n^2}, $$ 移项可得 $$ \displaystyle\sum\limits_{n\geqslant 1}\displaystyle\frac{1}{n^2}=\displaystyle\frac{\pi^2}6. $$ 这是著名的**二阶调和数**$H_\infty^{(2)}$，或者是**黎曼$\zeta$函数**在$z=2$处的函数值。这个式子的证明方式有很多，这是其中之一。在我的博客 [为了求正切的麦克劳林展开式，我复习了伯努利数](https://blog.csdn.net/qq_39677783/article/details/109827370) 中提到了另一种证法。 以上是周期为$2\pi$的周期函数的傅里叶级数，对于周期为$2l$的周期函数，我们当然也有类似的结果（并且狄利克雷定理仍然成立）： $$ f(x)\sim\displaystyle\frac{a_0}2+\displaystyle\sum\limits_{n\geqslant 1}\left(a_n\cos\displaystyle\frac{n\pi}lx+b_n\sin\displaystyle\frac{n\pi}lx\right), $$ 其中傅里叶系数为 $$ a_n=\displaystyle\frac1l\displaystyle\int_{-l}^lf(x)\cos\displaystyle\frac{n\pi}lx\mathrm dx,b_n=\displaystyle\frac1l\displaystyle\int_{-l}^lf(x)\sin\displaystyle\frac{n\pi}lx\mathrm dx. $$ ### 有限区间上的函数的傅里叶级数 如果函数$f(x)$仅仅定义在有限区间上，我们只要对它作平移和周期延拓，就可以将它变成定义在$\R$上的周期函数。这里给出有限区间$[-l,l]$上的函数的狄利克雷定理： 设函数$f(x)$是定义在$[-l,l]$上的分段可微函数，那么$f(x)$的傅里叶级数**收敛**于 $$ \begin{cases}\displaystyle\frac{f(x+0)+f(x-0)}{2},&x\in(-l,l),\\\displaystyle\frac{f(-l+0)+f(l-0)}{2},&x=\pm l.\end{cases} $$ 如果再增加$f(x)$在$[-l,l]$上连续的条件，并且$f(-l)=f(l)$（这是为了保证周期延拓后仍然连续），那么其傅里叶级数就在$[-l,l]$上**一致收敛**于$f(x)$。 特别地，我们提一下所谓的**正弦级数**、**余弦级数**的概念。对于定义在$[0,l]$上的函数，如果将它作偶延拓，就会得到一个偶函数$f_{\text{even}}(x)$，这时$f_{\text{even}}(x)\sin\displaystyle\frac{n\pi}lx$是奇函数，在对称区间$[-l,l]$上的积分恒为零，因此傅里叶系数$b_n$全都是零，傅里叶级数就只剩下了余弦函数，这就是所谓的**余弦级数**。类似地，作奇延拓之后得到的函数的傅里叶级数是**正弦级数**。 ### 傅里叶级数的复数形式 我们有著名的欧拉公式 $$ e^{ix}=\cos x+i\sin x, $$ 那么 $$ \cos x=\displaystyle\frac{e^{ix}+e^{-ix}}{2},\sin x=\displaystyle\frac{e^{ix}-e^{-ix}}{2i}. $$ 设定义在区间$[-l,l]$上的函数$f(x)$可以展开成傅里叶级数 $$ f(x)=\displaystyle\frac{a_0}2+\displaystyle\sum\limits_{n\geqslant 1}\left(a_n\cos n\omega x+b_n\sin n\omega x\right), $$ 其中$\omega=\displaystyle\frac{\pi}l$，并且系数 $$ a_n=\displaystyle\frac1l\displaystyle\int_{-l}^lf(x)\cos n\omega x\mathrm dx, $$ $$ b_n=\displaystyle\frac1l\displaystyle\int_{-l}^lf(x)\sin n\omega x\mathrm dx. $$ 现在将由复指数定义的三角函数代入傅里叶级数，有 $$ f(x)=\displaystyle\frac{a_0}2+\displaystyle\sum\limits_{n\geqslant 1}\displaystyle\frac{a_n-ib_n}2e^{in\omega x}+\displaystyle\sum\limits_{n\geqslant 1}\displaystyle\frac{a_n+ib_n}{2}e^{-in\omega x}. $$ 现在我们试图统一形式，方法是将第二个和式的正整数$n$进行翻转，变成负整数，于是指数上的负号就消失了： $$ \displaystyle\sum\limits_{n\geqslant 1}\displaystyle\frac{a_n+ib_n}{2}e^{-in\omega x}=\displaystyle\sum\limits_{n\leqslant-1}\displaystyle\frac{a_{-n}+ib_{-n}}{2}e^{in\omega x}. $$ 合并得 $$ f(x)=\displaystyle\sum\limits_{n=-\infty}^{\infty}F_ne^{in\omega x}, $$ 其中 $$ F_0=\displaystyle\frac{a_0}2,F_{\pm n}=\displaystyle\frac12(a_n\mp ib_n). $$ 这里的$n$取一切整数，而不仅仅是正整数。 ## 卷积 现在介绍一个至关重要的概念：卷积。 设函数$f(x),g(x)$在$\mathbb R$上可积并且绝对可积，那么积分 $$ \displaystyle\int_{-\infty}^{+\infty}f(t)g(x-t)\mathrm dt $$ 定义了一个新的函数，称为$f(x)$和$g(x)$的**卷积**，记作$(f*g)(x)$。 卷积满足诸多性质，如交换律 $$ f*g=g*f, $$ 结合律 $$ (f*g)*h=f*(g*h), $$ 对加法的分配律 $$ (f+g)*h=f*h+g*h. $$ 读者自证不难。 以上是对于实数集上的可积函数定义的卷积，而事实上还有很多不同类型的卷积。例如，对于数列$\{a_n\},\{b_n\}$，定义其卷积数列为 $$ (a*b)_n=\displaystyle\sum\limits_{k=0}^na_kb_{n-k}=\displaystyle\sum\limits_{i+j=n}a_ib_j, $$ 这是一个非常有趣的卷积，因为如果我们考虑这些数列的生成函数$A(x)=\displaystyle\sum\limits_{n\geqslant0}a_nx^n,B(x)=\displaystyle\sum\limits_{n\geqslant 0}b_nx^n$，有 $$ A(x)B(x)=\displaystyle\sum\limits_{n\geqslant 0}\displaystyle\sum\limits_{i+j=n}a_ib_jx^n, $$ 也就是说，序列的卷积对应了生成函数的乘积。 卷积相关的内容还不止于此，但限于篇幅，我们这里只介绍这两种与本文内容有关的卷积。 ## 傅里叶变换 考虑定义在$[-l,l]$上的$f(x)$的傅里叶级数的复数形式 $$ f(x)\sim\displaystyle\sum\limits_{n=-\infty}^{\infty}F_ne^{in\omega x,} $$ 其中$\omega=\displaystyle\frac{\pi}l,$并且系数 $$ F_{\pm n}=\frac{1}{2l}\int_{-l}^lf(x)e^{\mp in\omega x}\mathrm dx. $$ 将系数代入傅里叶级数，并记$\lambda_n=n\omega=\dfrac{n\pi}l,\Delta\lambda_n=\lambda_n-\lambda_{n-1}=\omega$，则 $$ \begin{aligned} f(x)&\sim\displaystyle\sum\limits_{n=-\infty}^{\infty}\frac{1}{2l}\int_{-l}^lf(\xi)e^{-in\omega(\xi-x)}\mathrm d\xi\\ &=\displaystyle\sum\limits_{n=-\infty}^{\infty}\frac{1}{2\pi}\int_{-l}^lf(\xi)e^{-i\lambda_n(\xi-x)}\mathrm d\xi\Delta\lambda_n \end{aligned} $$ 这看起来像是一个黎曼和。虽然这里的$\lambda_n$并不是有限区间上的分割，但是在$l\to+\infty$时，上式的一个合理的极限应该是这个积分 $$ \begin{aligned} f(x)&\to\frac{1}{2\pi}\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}f(\xi)e^{-i\lambda(\xi-x)}\mathrm d\xi\mathrm d\lambda\\ &=\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{i\lambda x}\int_{-\infty}^{+\infty}f(\xi)e^{-i\lambda\xi}\mathrm d\xi\mathrm d\lambda. \end{aligned} $$ 注意，如果将$x$看做时间，$f(x)$看做一个跟时间有关的函数，那么$\lambda$就是频率。 但是这又是一个猜想，它真的严谨吗？于是一个新的狄利克雷定理应运而生： 如果定义在$\mathbb R$上的函数$f(x)$在任何有限区间上是分段可微的，并且在$\mathbb R$上绝对可积，那么对任何$x\in\mathbb R$，有 $$ \frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{i\lambda x}\int_{-\infty}^{+\infty}f(\xi)e^{-i\lambda\xi}\mathrm d\xi\mathrm d\lambda=\frac{f(x+0)+f(x-0)}{2}. $$ 如果再增加$f(x)$连续的条件，那么上式等于$f(x)$。 定理的证明略去。 现在，我们将这个二重积分的内层积分提出来，令作$F(\lambda)$，就得到了如下公式： $$ F(\lambda)=\int_{-\infty}^{+\infty}f(\xi)e^{-i\lambda\xi}\mathrm d\xi, $$ $$ f(x)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}F(\lambda)e^{i\lambda x}\mathrm d\lambda. $$ 我们惊讶地发现，这个$F(\lambda)$就是我们要的频域函数，它称为$f(x)$的**傅里叶变换**，这一组公式称为**傅里叶变换的反演公式**。经过如此艰苦卓绝的数学推导，我们终于看到了结果。 现在，我们可以考虑它在信号处理上的实际意义了。如果用关于时间$t$的函数$f(t)$来表示一个信号，那么它的频域函数就是 $$ F(\omega)=\int_{-\infty}^{+\infty}f(t)e^{-i\omega t}\mathrm dt, $$ 而如果我们知道了频域函数$F(\omega)$，也可以反推出时域函数 $$ f(t)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}F(\omega)e^{i\omega t}\mathrm d\omega. $$ 我们可以将$F$记作$\mathrm{FT}[f]$（或$\mathrm{FT}[f(t)]$，下同），将$f$记作$\mathrm{IFT}[F]$，以此展示它们之间紧密的联系。 现在研究傅里叶变换所具有的性质。显然有线性性 $$ \mathrm{FT}[\alpha f+\beta g]=\alpha\mathrm{FT}[f]+\beta\mathrm{FT}[g] $$ 以及所谓的**频移特性**： $$ \mathrm{FT}[f](\omega+\omega_0)=\mathrm{FT}[f(t)e^{-i\omega_0t}](\omega), $$ **时移特性**： $$ \mathrm{FT}[f(t+t_0)](\omega)=\left(e^{i\omega t_0}\mathrm{FT}[f]\right)(\omega). $$ 读者自证不难。 除此之外，我们还有著名的**时域卷积定理**和**频域卷积定理**： $$ \mathrm{FT}[f*g]=\mathrm{FT}[f]\mathrm{FT}[g], $$ $$ \frac{1}{2\pi}\mathrm{IFT}[F*G]=\mathrm{IFT}[F]\mathrm{IFT}[G]. $$ 这里，我们只给出后者的证明，前者是类似的，读者自证不难。 **证** 记$\mathrm{IFT}[F]=f,\mathrm{IFT}[G]=g$，那么 $$ \begin{aligned} \mathrm{IFT}[F*G] &=\frac{1}{2\pi}\int_{-\infty}^{+\infty}(F*G)(\omega)e^{i\omega t}\mathrm d\omega\\ &=\frac{1}{2\pi}\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}F(u)G(\omega-u)\mathrm du e^{i\omega t}\mathrm d\omega\\ &=\frac{1}{2\pi}\int_{-\infty}^{+\infty}F(u)\mathrm du\int_{-\infty}^{+\infty}G(\omega-u)e^{i\omega t}\mathrm d\omega, \end{aligned} $$ 令$v=\omega-u$，则上式 $$ \begin{aligned} &=\frac{1}{2\pi}\int_{-\infty}^{+\infty}F(u)\mathrm du\int_{-\infty}^{+\infty}G(v)e^{i(u+v)t}\mathrm dv\\ &=\frac{1}{2\pi}\int_{-\infty}^{+\infty}F(u)e^{iut}\mathrm du\int_{-\infty}^{+\infty}G(v)e^{ivt}\mathrm dv\\ &=\frac{1}{2\pi}\cdot2\pi\mathrm{IFT}[F]\cdot 2\pi\mathrm{IFT}[G]\\ &=2\pi\mathrm{IFT}[F]\mathrm{IFT}[G]. \end{aligned} $$ 于是频域卷积定理得证。 除了这种傅里叶变换之外，对于定义在$[0,+\infty)$上的函数，作偶延拓或者奇延拓后再作类似的变换，就可分别得到函数的**余弦变换**和**正弦变换**。 向认真读到这里的读者表示钦佩，现在可以稍微休息一下啦~ ## 离散的情形 ### 采样 前面的变换都很美好，但我们得面对现实：你没法在计算机中存储连续定义的函数，更不能对它进行计算。那该如何表示一个信号？ 一个听起来很自然的想法：取定义域上的一列等间距的点，将这些点处的函数值保存下来，于是函数变成了数列。这就是采样。 但是这产生了两个问题： - 间距$T$（即**采样周期**）取何值，才能使得丢失的信息尽可能少？ - 采样对傅里叶变换的影响是什么？ 为了在数学上很好地定义采样，我们首先给出所谓的**单位冲激函数**$\delta(t)$的定义。单位冲激函数又称**狄拉克函数**，它由如下两个式子定义： $$ \begin{cases} \delta(t)=0,t\neq 0,\\ \displaystyle\int_{-\infty}^{+\infty}\delta(t)\mathrm dt=1. \end{cases} $$ 这是一个非常特殊的函数，它在非零处的值都是零，而零处的值并没有给出。但是它在整个实数轴上的积分是$1$。读者也许会猜想它在$0$处的值为无穷大，但这似乎与一般函数的定义不符。事实上它不是经典意义下的函数，而是一个**广义函数**，它是**单位阶跃函数** $$ u(t)=\begin{cases}1,&t>0\\0,&t<0\end{cases} $$ 的**广义导**。正因为它的特殊性，我们在使用它的时候必须结合其实际意义来考虑。例如，定义其傅里叶变换为 $$ \mathrm{FT}[\delta](\omega)=1 $$ 是合理的，因为考虑它与某个函数$f(t)$的卷积，可以使用分部积分法得到 $$ \begin{aligned} \int_{-\infty}^{+\infty}\delta(\tau)f(t-\tau)\mathrm d\tau &=\int_{-\infty}^{+\infty}f(t-\tau)\mathrm du(\tau)\\ &=f(t-\tau)u(\tau)\bigg|_{-\infty}^{+\infty}-\int_{-\infty}^{+\infty}u(\tau)\mathrm df(t-\tau)\\ &=f(-\infty)-\int_0^{+\infty}\mathrm df(t-\tau)\\ &=f(-\infty)-f(-\infty)+f(t)\\ &=f(t). \end{aligned} $$ 也就是说，$\delta(t)$与任何函数$f(t)$的卷积仍然得到$f(t)$本身。注意这里我们用了一个不太严谨的记号$f(-\infty)$，但无论如何它被抵消了。那么根据时域卷积定理，$\mathrm{FT}[\delta](\omega)$与其它傅里叶变换$F(\omega)$的乘积也应当是$F(\omega)$本身，所以$FT[\delta](\omega)=1$是合理的。再根据时移特性，我们还能得到$\delta(t-t_0)$的傅里叶变换，它是$e^{-i\omega t_0}$。那么对一个复指数信号$e^{-i\omega t_0}$进行逆变换就是 $$ \frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{-i\omega(t_0-t)}\mathrm d\omega=\delta(t-t_0). $$ 现在给出采样的数学定义。对一个信号$f(t)$采样后的结果是 $$ f^*(t)=f(t)\sum\limits_{n=-\infty}^{\infty}\delta(t-nT) $$ 它表示为$f(t)$与一个周期冲激串$s(t)=\displaystyle\sum\limits_{n=-\infty}^{\infty}\delta(t-nT)$的乘积。其中，$T$是**采样周期**（或采样间隔），即每隔$T$时间采一个点。通常记$f_s=\dfrac1T$为**采样率**（或采样频率），即每秒钟采了$f_s$个点；记$\omega_s=\dfrac{2\pi}T$为**采样角频率**。 为了计算$f^*(t)$的傅里叶变换，可以使用频域卷积定理。首先计算$s(t)$的傅里叶系数，即如果 $$ s(t)=\sum\limits_{n=-\infty}^{\infty}F_ne^{in\omega_st},\omega_s=\frac{2\pi}{T} $$ 则 $$ \begin{aligned} F_n&=\frac1T\int_{-T/2}^{T/2}\sum\limits_{k=-\infty}^{\infty}\delta(t-kT)e^{-in\omega_st}\mathrm dt\\ &=\frac1T\int_{-T/2}^{T/2}\delta(t)e^{-in\omega_st}\mathrm dt\\ &=\frac1T\mathrm{FT}[\delta](n\omega_s)\\ &=\frac1T. \end{aligned} $$ 于是 $$ s(t)=\frac1T\sum\limits_{n=-\infty}^{\infty}e^{in\omega_st} $$ 那么$s(t)$的傅里叶变换$S(\omega)$是 $$ \begin{aligned} S(\omega)&=\int_{-\infty}^{+\infty}s(t)e^{-i\omega t}\mathrm dt\\ &=\int_{-\infty}^{+\infty}\frac1T\sum\limits_{n=-\infty}^{\infty}e^{-it(\omega-n\omega_s)}\mathrm dt\\ &=\frac1T\sum\limits_{n=-\infty}^{\infty}\int_{-\infty}^{+\infty}e^{-it(\omega-n\omega_s)}\mathrm dt\\ &=\frac{2\pi}{T}\sum\limits_{n=-\infty}^{\infty}\delta(n\omega_s-\omega)\\ &=\frac{2\pi}{T}\sum\limits_{n=-\infty}^{\infty}\delta(\omega-n\omega_s), \end{aligned} $$ 它仍然是周期冲激串。根据频域卷积定理，$f^*(t)$的傅里叶变换是 $$ \begin{aligned} F^*(\omega)&=\frac{1}{2\pi}\int_{-\infty}^{+\infty}F(\omega-u)S(u)\mathrm du\\ &=\frac{1}{2\pi}\int_{-\infty}^{+\infty}F(\omega-u)\frac{2\pi}{T}\sum\limits_{n=-\infty}^{\infty}\delta(u-n\omega_s)\mathrm du\\ &=\frac1T\sum\limits_{n=-\infty}^{\infty}\int_{-\infty}^{+\infty}F(\omega-u)\delta(u-n\omega_s)\mathrm du. \end{aligned} $$ 如果令$v=u-n\omega_s$，那么后面的积分是 $$ \int_{-\infty}^{+\infty}F(\omega-n\omega_s-v)\delta(v)\mathrm dv, $$ 这是$F(\omega-n\omega_s)$与狄拉克函数的卷积，前面已经知道这个卷积的结果就是$F(\omega-n\omega_s)$自己，于是 $$ F^*(\omega)=\frac1T\sum\limits_{n=-\infty}^{\infty}F(\omega-n\omega_s). $$ 这是一个至关重要的式子，它表明：在时域对信号以$T$为周期进行采样，相当于在频域按$\omega_s$为周期进行延拓。  这张图片来自Oppenheim的《离散时间信号处理》。图(a)是一个简单的信号傅里叶变换后的结果，其中$\Omega_N$称为**截止频率**。图(b)就是我们刚刚给出的$S(\omega)$，其中$\Omega_s$就是我们的采样角频率$\omega_s$。由于采样的傅里叶变换就是将原信号的傅里叶变换与一个周期冲激串求卷积，并除以$2\pi$，当$\Omega_s\geqslant 2\Omega_N$时采样后的傅里叶变换如图(c)所示，当$\Omega_s<2\Omega_N$时采样后的傅里叶变换如图(d)所示。可以看到，图(c)的频谱在每一个整数倍的$\Omega_s$上仍然保持一个与原信号的频谱完全一样的副本（并在幅度上除以$T$），而图(d)的频谱发生了**混叠**，不再含有原频谱的完整信息。于是我们有如下定理： 对于一个信号$x(t)$，假设其截止频率为$\omega_N$，现对其按角频率$\omega_s=\dfrac{2\pi}{T}$进行采样，则当且仅当 $$ \omega_s\geqslant 2\omega_N $$ 时，这个信号能唯一地由它的样本 $$ x^*(n)=x(nT),n\in\mathbb Z $$ 所确定。这个定理称为**奈奎斯特采样定理**（Nyquist-Shannon Sampling Theorem）。频率$\omega_N$一般称为**奈奎斯特频率**，而频率$2\omega_N$称为**奈奎斯特率**。我们称$\omega_s<2\omega_N$的采样为**欠采样**，这会使信号的频谱发生混叠。一个直观的例子是高速旋转的车轮看上去向相反方向转。 到此，我们对信号进行采样的想法已经有了基本的理论支持。那么如何从一个采样后的信号恢复出原信号呢？一个简单的办法是先对采样后的信号作傅里叶变换，在频域上截取$[-\omega_s,\omega_s]$这个区间，再作逆傅里叶变换即可。 ### 离散时间傅里叶变换 现在我们有一个采样后的信号 $$ \begin{aligned} f^*(t)&=f(t)\sum\limits_{n=-\infty}^{\infty}\delta(t-nT)\\ &=\displaystyle\sum\limits_{n=-\infty}^{\infty}f(nT)\delta(t-nT). \end{aligned} $$ 试着对它进行傅里叶变换，有 $$ \begin{aligned} F^*(\omega)&=\int_{-\infty}^{+\infty}\displaystyle\sum\limits_{n=-\infty}^{\infty}f(nT)\delta(t-nT)e^{-i\omega t}\mathrm dt\\ &=\displaystyle\sum\limits_{n=-\infty}^{\infty}f(nT)\int_{-\infty}^{+\infty}\delta(t-nT)e^{-i\omega t}\mathrm dt. \end{aligned} $$ 注意，我们交换了求和符号和积分符号，这要求级数的一致收敛性，但这里我们忽略了这一点。 注意到冲激信号$\delta(t)$的傅里叶变换是常函数$1$，即 $$ \int_{-\infty}^{+\infty}\delta(t)e^{-j\omega t}\mathrm dt=1. $$ 所以 $$ \begin{aligned} F^*(\omega)&=\displaystyle\sum\limits_{n=-\infty}^{\infty}f(nT)\int_{-\infty}^{+\infty}\delta(t)e^{-i\omega(t+nT)}\mathrm dt\\ &=\displaystyle\sum\limits_{n=-\infty}^{\infty}f(nT)e^{-i\omega nT}\int_{-\infty}^{+\infty}\delta(t)e^{-i\omega t}\mathrm dt\\ &=\displaystyle\sum\limits_{n=-\infty}^{\infty}f(nT)e^{-i\omega nT}. \end{aligned} $$ 这个变换称为**离散时间傅里叶变换**（Discrete Time Fourier Transform）。注意，这个函数是以$\omega_s=\dfrac{2\pi}T$为周期的函数，因为 $$ F^*(\omega+\omega_s)=\displaystyle\sum\limits_{n=-\infty}^{\infty}f(nT)e^{-i(\omega+\omega_s)nT}=\displaystyle\sum\limits_{n=-\infty}^{\infty}f(nT)e^{-i\omega nT}\cdot e^{-i2\pi n}, $$ 而$e^{-i2\pi n}=1$，所以$F^*(\omega+\omega_s)=F^*(\omega)$。 但是我们想一想，采样实际上是把连续的函数（这里指定义域连续）变成了离散的序列，那么我们完全可以将周期$T$看做“$1$个单位时间”，而$nT$就变成了“第$n$个采样点”，不妨设$x(n)=f(nT)$，那么我们有 $$ X(\omega)=\displaystyle\sum\limits_{n=-\infty}^{\infty}x(n)e^{-i\omega n}. $$ 这里将$T$看作$1$，也就是将采样频率$f_s$看作$1$，这是所谓的**频率归一化**。 既然有了DTFT，自然也应该有**逆离散时间傅里叶变换**IDTFT。根据定义有 $$ f(nT)=\frac{1}{\omega_s}\int_{-\omega_s/2}^{\omega_s/2}F(\omega)e^{i\omega nT}\mathrm d\omega. $$ 频率归一化后就是 $$ x(n)=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(\omega)e^{i\omega n}\mathrm d\omega. $$ ### 离散傅里叶变换 有了DTFT，我们可以将离散的信号转换为频域信号了。但是这又带来了问题。一方面，输入的时域信号$x(n)$是一个无限长的序列，这仍然无法在计算机中存储。另一方面，输出的$X(\omega)$仍然是一个定义域连续的函数。 第一个问题很好解决，因为实际情况下信号都是有限长的，$n$不会取全体整数$\mathbb Z$而只能取某个范围（不妨记为$[0,L-1]$）中的整数。 第二个问题则需要动动脑子。也许我们可以考虑，在频谱的一个区间$[0,\omega_s]$中只求出等距的$N$个点的值，这感觉有点像在频谱上做“采样”。既然我们已经做了频率归一化，$\omega_s$就是$2\pi$。那么第$k$个频率是 $$ k\cdot\frac{2\pi}{N}, $$ 于是，$X(\omega)$就变成了$X(k)$，即所谓的“第$k$个频率”： $$ X(k)=\displaystyle\sum\limits_{n=0}^{L-1}x(n)e^{-ink\frac{2\pi}{N}},k=0,1,2,\cdots,N-1. $$ 如果令$W_N^k=e^{-ik\frac{2\pi}{N}}$，那么 $$ X(k)=\displaystyle\sum\limits_{n=0}^{L-1}x(n)W_N^{kn}. $$ 但是要注意，这个$W_N^k$并不是真正的第$k$个$N$次单位根，而是其共轭，不过因为它们具有非常类似的性质，所以我们采用了这个记号。不用$\omega_N^{kn}$也是为了防止跟其它的表示频率的$\omega$混淆起来。 到这里问题还没有完全解决，因为我们还不知道$N$是什么。一个合理的猜想是$N=L$，因为如果$N>L$，我们可以在输入序列$x$的末尾补上$N-L$个$0$，这不会影响结果。而如果$N<L$，会发现 $$ \begin{aligned} X(k)&=\displaystyle\sum\limits_{n=0}^{N-1}x(n)W_N^{kn}+\displaystyle\sum\limits_{n=N}^{L-1}x(n)W_N^{kn}\\ &=\displaystyle\sum\limits_{n=0}^{N-1}x(n)W_N^{kn}+\displaystyle\sum\limits_{n=0}^{L-N-1}x(N+n)W_N^{k(n+N)}, \end{aligned} $$ 其中$W_N^{k(n+N)}=W_N^{kn}$。因此 $$ X(k)=\displaystyle\sum\limits_{n=0}^{N-1}x(n)W_N^{kn}+\displaystyle\sum\limits_{n=0}^{L-N-1}x(N+n)W_N^{kn}. $$ 观察发现，这相当于将输入序列中的每个元素$x(n)$都变为$x(n\bmod N)$，并叠加，再进行变换，这是所谓的**时域混叠**。因此$N$必须大于等于$L$，那么为了方便，取$N=L$。于是最终我们得到了 $$ X(k)=\displaystyle\sum\limits_{n=0}^{N-1}x(n)W_N^{kn},k=0,1,2,\cdots,N-1. $$ 这就是所谓的**离散傅里叶变换**（Discrete Fourier Transform），也就是我们在计算机中真正采用的傅里叶变换。 但到这里还没结束，因为我们还需要其逆变换。根据连续函数（指定义域连续）的傅里叶变换的反演公式，我们有理由猜想 $$ x(k)=\frac{1}{N}\displaystyle\sum\limits_{n=0}^{N-1}X(n)W_N^{-kn}. $$ 现在尝试证明。注意到一件有趣的事：从$x$到$X$的变换，可以由这个$N\times N$的矩阵 $$ A=\begin{bmatrix} 1&1&1&\cdots&1\\ 1&W_N^1&W_N^2&\cdots&W_N^{N-1}\\ 1&W_N^2&W_N^4&\cdots&W_N^{2(N-1)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&W_N^{N-1}&W_N^{2(N-1)}&\cdots&W_N^{(N-1)(N-1)} \end{bmatrix}, $$ 即$A_{ij}=W_N^{ij}$（下标从$0$开始）表示： $$ \mathbf{X}=A\mathbf{x}. $$ 现在要求逆变换，就是要求矩阵$A$的逆矩阵$B$。注意到当$B$满足$B_{ij}=\displaystyle\frac1NW_N^{-ij}$时，有 $$ \displaystyle\sum\limits_{n=0}^{N-1}B_{in}A_{nj}=\displaystyle\sum\limits_{n=0}^{N-1}\frac{1}{N}W_N^{-in}W_N^{nj}=\frac{1}{N}\displaystyle\sum\limits_{n=0}^{N-1}W_N^{n(j-i)}. $$ 注意，这里的$i$和$j$是两个$[0,N-1]$中的整数，而不是虚数单位。上式是一个等比数列求和的形式，当$j=i$时每一项都是$1$，结果为$1$。当$j\neq i$时，它等于 $$ \frac1N\cdot\frac{1-W_N^{N(j-i)}}{1-W_N^{j-i}}=0. $$ 因此，乘积$BA$的第$i$行第$j$列的元素的确等于$[i=j]$，所以$B$是$A$的逆矩阵，则由$B$确定的变换 $$ x(k)=\frac1N\displaystyle\sum\limits_{n=0}^{N-1}X(n)W_N^{-kn} $$ 的确是DFT的逆变换IDFT。 现在我们有了离散傅里叶变换，可以类比傅里叶变换的性质，有以下的**频移特性** $$ \mathrm{DFT}[x](k+k_0)=\mathrm{DFT}[x(n)W_N^{k_0n}](k) $$ 以及**时移特性** $$ \mathrm{DFT}[x(n+n_0)]=W_N^{-kn_0}\mathrm{DFT}[x], $$ 注意这里的移动都是循环移动，例如$x(n+n_0)=x((n+n_0)\bmod N)$。 在离散卷积的意义下，我们自然也有对应的**时域卷积定理**和**频域卷积定理**，只是有一个地方要注意：我们的信号长度不是无穷大了。比方说考虑两个序列的DFT之积： $$ \begin{aligned} \left(\mathrm{DFT}[x]\mathrm{DFT}[y]\right)(k) &= \displaystyle\sum\limits_{i=0}^{N-1}x(i)W_N^{ik}\displaystyle\sum\limits_{j=0}^{N-1}y(j)W_N^{jk}\\ &= \displaystyle\sum\limits_{m=0}^{2(N-2)}W_N^{mk}\displaystyle\sum\limits_{i=0}^{\min\{m,N-1\}}x(i)y(m-i). \end{aligned} $$ 会发现，它们的乘积有$2N-1$项，但是每一个序列都只有$N$项，于是当$m>N-1$时，后面将不再是完整的卷积形式。但是如果只取乘积的前$N$项，我们有 $$ \left(\mathrm{DFT}[x]\mathrm{DFT}[y]\right)_N(k)=\displaystyle\sum\limits_{m=0}^{N-1}(x*y)(m)W_N^{mk}=\mathrm{DFT}_N[x*y](k). $$ 这是时域卷积定理，类似地也有频域卷积定理。 ### 快速傅里叶变换 朴素的计算DFT的算法显然是$O(N^2)$的，但是我们有更快的办法。 考虑单位根的性质：对于偶数$N$有 $$ W_N^k=W_{N/2}^{k/2}=-W_N^{k+N/2}, $$ 如果$N$是奇数，就在后面补一个$0$即可。 我们有 $$ \begin{aligned} X(k)&=\displaystyle\sum\limits_{n=0}^{N-1}x(n)W_N^{kn}\\ &=\displaystyle\sum\limits_{0\leqslant 2n+1\leqslant N-1}x(2n+1)W_N^{2kn+k}+\displaystyle\sum\limits_{0\leqslant 2n\leqslant N-1}x(2n)W_N^{2kn}\\ &=W_N^k\displaystyle\sum\limits_{n=0}^{N/2-1}x(2n+1)W_{N/2}^{kn}+\displaystyle\sum\limits_{n=0}^{N/2-1}x(2n)W_{N/2}^{kn}. \end{aligned} $$ 也就是说，将原序列$x(n)$的奇子列和偶子列分别作为两个新的序列进行DFT，再合并结果即可。合并的时间复杂度为$O(N)$，因此该算法运行时间为 $$ T(N)=2T(N/2)+O(N). $$ 根据主定理（或者其它办法），易知时间复杂度为$O(N\log N)$。 至于逆傅里叶变换，只要将$W_N^k$换成$W_N^{-k}$，并将结果除以$N$即可。 ### 短时傅里叶变换 有时候，我们可能有一段很长的信号，对它直接作傅里叶变换会得到整个信号的频谱，但这样我们就丢失了时间上的信息：我们并不知道这个频谱中的某些频率是何时出现的。为此我们可以每次截取信号的一小段作傅里叶变换，得到不同时段的傅里叶变换的结果。 分段的过程本质上是给一段信号乘上一个取值不为零的函数，这个函数称为**窗函数**。最简单的窗函数是矩形窗函数$w(n)=1$，它相当于直接截取了这段信号。 实际应用中需要考虑信号的特殊性，为了提高频率分辨率，可能需要更复杂的窗函数，例如汉宁窗（Hanning Window） $$ w(n)=\frac12\left(1-\cos\frac{2\pi(n-1)}{N}\right), $$ 还有高斯窗、布莱克曼窗等等。 这样的变换称为**短时傅里叶变换**（Short Time Fourier Transform, STFT）。 ## 简单的应用 应用主要是两块，一块是在信号处理上的应用，一块是在数学上的应用。 ### 电话按键音识别 很久以前有这样一个新闻：一段记者采访360董事长周鸿祎的视频在网上流传开来，视频中有一串电话按键音引起了南京大学一名学生的注意。他通过这段电话按键音识别出了周鸿祎的手机号码，还因此获得了李开复创新工场的邀请。  而事实上这个电话按键音识别一点都不复杂。电话按键音是由两个频率的正弦波叠加而成的： | Frequency (Hz) | 1209 | 1336 | 1477 | | :------------: | :--: | :--: | :--: | | 697 | 1 | 2 | 3 | | 770 | 4 | 5 | 6 | | 852 | 7 | 8 | 9 | | 941 | * | 0 | # | 所以只要截取一段按键音频信号，用FT变换到频域，答案就十分显然了。 ### 频分复用 可以在一条音频线路上同时传输多个音频信号。具体的方法是：将这些信号变换到频域，平移到互不干扰的位置进行叠加，再变换回时域进行传输。例如在采样率为44100Hz的线路上，可以同时传输三个采样率为14000Hz的信号。 ### 多项式乘法 考虑两个多项式$A(x)=\displaystyle\sum\limits_{k=0}^{n-1}a_kx^k,B(x)=\displaystyle\sum\limits_{k=0}^{n-1}b_kx^k$，我们希望求出它们的乘积 $$ A(x)B(x)=\displaystyle\sum\limits_{k=0}^{2n-2}c_kx^k, $$ 其中 $$ c_k=\displaystyle\sum\limits_{i+j=k}a_ib_j. $$ 朴素的算法时间复杂度当然是$O(n^2)$。为了寻求更快的速度，我们首先要改变多项式的表示方式。原来我们用序列$a_k,k=0,1,\cdots,n-1$表示一个多项式的各项系数，从而确定一个多项式。但是依拉格朗日插值法，任意$n$个横坐标不同的点都能唯一确定一个次数不超过$n-1$的多项式，于是我们也可以用这样的**点-值对** $$ \left(x_0,A(x_0)\right),\left(x_1,A(x_1)\right),\cdots,\left(x_{n-1},A(x_{n-1})\right) $$ 来表示一个多项式，这样的表示方法称为**点值表示法**。 如果我们还能求出$B(x)$在$x_0,x_1,\cdots,x_{n-1}$处的值，也就是把$B(x)$也表示为 $$ \left(x_0,B(x_0)\right),\left(x_1,B(x_1)\right),\cdots,\left(x_{n-1},B(x_{n-1})\right), $$ 那么，计算$A(x)$和$B(x)$的乘积的点值表示是$O(n)$的，只要将每一个$A(x_k)B(x_k)$乘起来即可。注意，由于乘积多项式的次数是$2n-2$，需要$2n-1$个点值对，因此我们也得求出$A(x)$和$B(x)$的$2n-1$个点值对。 所以我们的步骤是：先求出$A(x)$和$B(x)$的点值表示，相乘，再变回系数表示。但是如何快速求出点值表示呢？ 回顾DFT： $$ X(k)=\displaystyle\sum\limits_{n=0}^{N-1}x(n)W_N^{nk}, $$ 不难发现，它实际上就是在求多项式$P(z)=\displaystyle\sum\limits_{n=0}^{N-1}x(n)z^n$在$z=W_N^0,W_N^1,\cdots,W_N^{N-1}$处的$N$个值，也就是说对序列进行DFT就可以得到其点值表示。于是利用快速傅里叶变换算法，我们可以在$O(n\log n)$的时间内求出点值表示，用$O(n)$的时间计算乘积，再用$O(n\log n)$的时间将乘积的点值表示变换为系数表示，于是总时间复杂度$O(n\log n)$。 事实上多项式乘法本质上是在求卷积，根据时域卷积定理，我们可以对其做DFT变换到频域，在频域上做乘积，再变换回时域，这和我们刚才的推导是一致的。 ### 范德蒙德矩阵求逆 刚才在计算IDFT时我们用到了这样一个矩阵 $$ A=\begin{bmatrix} 1&1&1&\cdots&1\\ 1&W_N^1&W_N^2&\cdots&W_N^{N-1}\\ 1&W_N^2&W_N^4&\cdots&W_N^{2(N-1)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&W_N^{N-1}&W_N^{2(N-1)}&\cdots&W_N^{(N-1)(N-1)} \end{bmatrix}. $$ 这是一个**范德蒙德矩阵**。一般地，设有一列数$\alpha_1,\alpha_2,\cdots,\alpha_m$，矩阵 $$ \begin{bmatrix} 1&\alpha_1&\alpha_1^2&\cdots&\alpha_1^{n-1}\\ 1&\alpha_2&\alpha_2^2&\cdots&\alpha_2^{n-1}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\alpha_m&\alpha_m^2&\cdots&\alpha_m^{n-1} \end{bmatrix} $$ （或其转置）是一个**范德蒙德矩阵**。当$m=n$时这是一个方阵，具有行列式 $$ \prod\limits_{1\leqslant j<i\leqslant n}(\alpha_i-\alpha_j). $$ 对这样一个矩阵求逆，本质上是一个多项式插值/多点求值问题，具体可以参考[gzy神仙的博客](https://www.cnblogs.com/gzy-cjoier/p/9741950.html)。 **** 参考资料： 程艺, 陈卿, 李平《数学分析讲义》第一册、第二册 NOI WC2018 《傅里叶变换及其在OI中的应用》by 王逸松 傅里叶变换 - Mathworks 中国（MATLAB官方文档） 胡广书《数字信号处理——理论、算法与实现》第三版 Alan V. Oppenheim, Ronald W. Schafer 《离散时间信号处理》第三版 上海科技大学《信息科学与技术导论》课程课件 by 娄鑫