前言
本篇会以多项式乘法作为引子,进一步了解快速傅立叶变换的原理和实现,从而深入理解傅立叶变换以及广义傅立叶变换的原理,为后面进一步推导和理解球谐函数打一下基础。
多项式乘法
已知
\left\{ \begin{array}{lr} f(x) = \sum_{i=0}^{i \lt n} a_i \cdot x^i = a_0 + a_1x^1 + a_2x^2 + \cdots + a_{n-1}x^{n-1} \\ g(x) = \sum_{i=0}^{i \lt m} b_i \cdot x^i = b_0 + b_1x^1 + b_2x^2 + \cdots + b_{m-1}x^{m-1} \end{array} \right. 求
F(x) = f(x) \cdot g(x) = \sum_{i=0, j=0}^{i \lt n, j \lt m}a_ib_jx^{i+j} \tag{1.1} 朴素的乘法
从
((1.1))不难看出,如果我们把每一项展开两两相乘的话,拢共需要计算
(n \cdot m)次。
举个具体的例子。
设
\left\{ \begin{array}{lr} f(x) = 1 + 2x + 3x^2 , & (a_i \in \{1, 2, 3\}, n = 3) \\ g(x) = 10 + 20x + 30x^2 + 40x^3, & (b_i \in \{10, 20, 30, 40\}, m = 4) \end{array} \right. \tag{1.2} 则有
\begin{array}{l} F(x) & = (1 + 2x + 3x^2) \cdot (10 + 20x + 30x^2 + 40x^3) \\ & = 10 + 20x + 30x^2 + 40x^3 \\ & + 20x + 40x^2 + 60x^3 + 80x^4 \\ & + 30x^2 + 60x^3 + 90x^4 + 120x^5 \\ & = 10 + 40x + 100x^2 + 160x^3 + 170x^4 + 120x^5, & (a_i \in \{10, 40, 100, 160, 170, 120\}, n = 6) \end{array} 显然,总共进行了
((n \cdot m = 3 \cdot 4 = 12))次计算,最终得到的结果有
(((n-1) + (m-1) + 1 = 2 + 3 + 1= 6))项,即最高次数为5次方。
这便是最朴素的多项式乘法了。接下来要做的事情就是让这件朴素的事情花哨起来,这就不得不从多项式的表示方式入手了。
多项式的表示
多项式可以用系数表示法
和点值表示法
来表示。
系数表示法 对于上面的多项式f(x)
,可以使用一个系数向量来唯一确定,即
\vec{a} = (a_0, a_1, a_2, \ldots, a_{n-1}) \tag{1.3} 点值表示法 我们如果可以取f(x)
上的n个自变量地方的值,那么这n个点值对
也可以唯一确定这个多项式。
设
(f_i = f(x_i), i \in \{0, 1, 2, \cdots, n-1\}),则有
\vec{(x, f)} = \{(x_0, f_0), (x_1, f_1), \ldots, (x_{n-1}, f_{n-1}) \} \tag{1.4} 之所以这n个点值对可以唯一确定这个多项式,是因为我们将其分别代入f(x)
,便可得到一个未知数为
((a_0, a_1, \cdots, a_{n-1})) 的n元一次线性方程组
,即
\left\{ \begin{array}{lr} f_0 &=& a_0 &+& a_1x_0 &+& a_2x_0^2 &+& \cdots &+& a_{n-1}x_0^{n-1} \\ f_1 &=& a_0 &+& a_1x_1 &+& a_2x_1^2 &+& \cdots &+& a_{n-1}x_1^{n-1} \\ & \vdots \\ f_{n-1} &=& a_0 &+& a_1x_{n-1} &+& a_2x_{n-1}^2 &+& \cdots &+& a_{n-1}x_{n-1}^{n-1} \\ \end{array} \right. 使用求解多元一次线性方程组的方法,如Gauss消元、LU分解、Crammer法则等[1],便可求出该多项式所有的系数
((a_0, a_1, a_2, \cdots, a_{n-1}))。
通过这两种表示方式,我们可以寻找到求多项式乘法的另一种方式。
另一种乘法方式
现在我们已知的是
(f(x))和
(g(x))的系数表达式,即
\left\{ \begin{array}{lr} f(x) \Leftrightarrow \vec{a} = (a_0, a_1, a_2, \ldots, a_{n-1}) \\ g(x) \Leftrightarrow \vec{b} = (b_0, b_1, b_2, \ldots, b_{m-1}) \end{array} \right. 我们需要求的是
(F(x))的系数表达式,即
\[ F(x) \Leftrightarrow \vec{c} = (c_0, c_1, c_2, \ldots, c_{k-1}) \text{, (k = (n-1) + (m-1) + 1)} \]
其中
(k = (n-1) + (m-1) + 1),即最高次项次数 + 1
,即
(F(x))的项数(包括系数为0的低次项)。
现在来考虑另一种思路,如果我们知道
(F(x))上的
(k)个点的值,组成
(k)个点值对,那么我们便可以通过上面提到的解多元一次线性方程组,来得到
(F(x))的系数表达式
(\vec{c})。设这
(k)个点值对为
\vec{(x, F)} = \{(x_0, F_0), (x_1, F_1), \ldots, (x_{k-1}, F_{k-1}) \} \tag{1.5} 很显然,对于任意定义域内的
(x_i),都有
F(x_i) = f(x_i) \cdot g(x_i) 那么
((1.5))便可转化为
\begin{array}{l} \vec{(x, F)} & = & \{(x_0, f_0) \cdot (x_0, g_0), (x_1, f_1) \cdot (x_1, g_1), \ldots, (x_{k-1}, f_{k-1}) \cdot (x_{k-1}, g_{k-1})\} \\ & = & \vec{(x, f)} * \vec{(x, g)} \tag{1.6} \end{array} (注意这里的
(*)并非向量积,而是向量的乘法,即各个分量相乘,得到相同项数的结果向量。)
那么我们只要分别在
(f(x))和
(g(x))上取
(k)个点值对,组成它们的点值表示法,然后将它们相乘,便可得到被积函数
(F(x))的点值表示法,然后再解多元一次线性方程组求出
(F(x))的系数表达式。
还拿上面的
((1.2))为例。
设
\left\{ \begin{array}{lr} f(x) \Leftrightarrow \vec{a} = (1, 2, 3) &, n = 3\\ g(x) \Leftrightarrow \vec{b} = (10, 20, 30 ,40) &, m = 4 \end{array} \right. 则要求的
(F(x))为
\[ F(x) \Leftrightarrow \vec{c} = (c_0, c_1, c_2, \ldots, c_{k-1}) \text{, (k = (n-1) + (m-1) + 1) = 6} \]
那么我们接下来按上面的步骤来求一下
(F(x))。
第一步 取
(k)个点,即
(x \in \{0, 1, 2, 3, 4, 5\}),求得
(f(x))和
(g(x))的点值表达式,即
\begin{array}{l} \vec{(x, f)} = \{(0, f(0)), \ldots, (5, f(5)) \} = \{(0, 1), &(1, 6), &(2, 17), &(3, 34), &(4, 57), &(5, 86) \} \\ \vec{(x, g)} = \{(0, g(0)), \ldots, (5, g(5)) \} = \{(0, 10), &(1, 100), &(2, 490), &(3, 1420), &(4, 3130), &(5, 5860) \} \end{array} 该步需要执行(
(kn+km))次乘法,是因为通过秦九韶算法
[2],每一个
(f(x_i))都只需要n次乘法得出,每一个
(g(x))都只需要m次乘法得出。
根据系数表达式求其点值表达式的过程叫做求值
。
第二步 根据上一步求出来的
(f(x))和
(g(x))的点值表达式,求出
(F(x))的点值表达式,即
\begin{array}{l} \vec{(x, F)} & = & \vec{(x, f)} * \vec{(x, g)} \\ & = & \{(0, f(0)*g(0)), \ldots, (5, f(5)*g(5)) \} \\ & = & \{(0, 10), (1, 600), (2, 8330), (3, 48280), (4, 178410), (5, 503960) \} \end{array} 这一步只需要(
(k))次乘法。
第三步 将上同求出来的点值表示法,通过解多元一次线性方程组的方法转为系数表示法即可。
将
((x, F))的元素分别代入
(F(x))中,得
\left\{ \begin{array}{l} 10 &=& c_0 \\ 600 &=& c_0 &+& c_1 &+& c_2 &+& c_3 &+& c_4 &+& c_5 \\ 8330 &=& c_0 &+& c_12 &+& c_22^2 &+& c_32^3 &+& c_42^4 &+& c_52^5 \\ 48280 &=& c_0 &+& c_13 &+& c_23^2 &+& c_33^3 &+& c_43^4 &+& c_53^5 \\ 178410 &=& c_0 &+& c_14 &+& c_24^2 &+& c_34^3 &+& c_44^4 &+& c_54^5 \\ 503960 &=& c_0 &+& c_15 &+& c_25^2 &+& c_35^3 &+& c_45^4 &+& c_55^5 \\ \end{array} \right. 使用高斯消元法
[3]转换成行阶梯式解方程组,可得
F(x) \Leftrightarrow \vec{c} = (10, 40, 100, 160, 170, 120) 与上面朴素的乘法得到了相同的结果。
高斯消元的时间复杂度为
(O(k^3)),由于在计算行阶梯式时每行的计算互不干扰,所以可使用多线程算法将复杂度降到
(O(k^2))根据点值表达式求其系数表达式的过程叫做插值
。
至此,我们通过求值(将系数转为点值)
、点值相乘
、插值(将点值转为系数)
,三个步骤,消耗了比朴素的乘法
((O(nm)))更高的复杂度
(O(kn+km) + O(k) + O(k^3))。
仿佛有些得不偿失-_-。
但是,这给优化时间复杂度提供了可操作的空间。因为如果我们选择一些特殊的点
来进行插值
和求值
的话,便可将其复杂度降低到比朴素乘法更低的程度。
而单位复数根
就是一种可以达成这种效果的特殊的点
。
复数
定义
复数[4]是一个数域(C
),它是实数域(R
)的延伸,额外包含带有虚数单位i
的数。其中虚数单位i
是-1的平方根
,即
((i^2 = -1))。
对于任意复数t
都可记为
t = a + bi \tag{2.1} 其中a、b都是实数,a称为t的实部,b称为t的虚部。所以,对于任意复数t
,都可以由唯一有序对表示,即
t \Leftrightarrow (a, b) \tag{2.2} 在二维笛卡尔坐标系中,如果我们将x轴作用于实部a,将y轴作用于虚部b,便可在该坐标系下唯一确定一点与唯一一个复数一一对应,这个坐标系称为复平面。
我们也可以通过极角坐标系
的方式表示复平面上的点。设
(r)为与坐标系原点的距离,
(\theta)为与x轴正方向的夹角且
(\theta \in [0, 2\pi)),则有序对
((r, \theta))也可唯一确定复平面上的一点,即可唯一确定一个与之对应的复数。如图所示。
复平面
即有
t \Leftrightarrow (r, \theta) \tag{2.3} 此时
\left\{ \begin{array}{l} a = r \cdot \cos{\theta} \\ b = r \cdot \sin{\theta} \end{array} \right. 即
t = r(\cos{\theta} + i \cdot \sin{\theta}) \tag{2.4} 同理也有
\left\{ \begin{array}{l} r = \sqrt{a^2 + b^2} \\ \theta = \arctan{\frac{b}{a}} (a \neq 0) or \frac{\pi}{2}(a = 0) \end{array} \right. 欧拉公式[5]将复数与三角函数关联起来,它的定义如下。
e^{i\theta} = \cos{\theta} + i \cdot \sin{\theta} \tag{2.5} 特别地,当
(\theta = \pi)时,有欧拉恒等式
e^{i\pi} + 1 = 0 \tag{2.6} 那么复数t可根据欧拉公式表示为指数形式,即
t = re^{i \theta} \tag{2.7} 运算
设
(t_1 \Leftrightarrow (a_1, b_1) \Leftrightarrow (r_1, \theta_1), t_2 \Leftrightarrow (a_2, b_2) \Leftrightarrow (r_2, \theta_2))。
加法
t_1 + t_2 = (a_1+a_2) + (b_1 + b_2)i \Leftrightarrow (a_1+a_2, b_1+b_2) \tag{2.8} 乘法
\begin{array}{l} t_1 \cdot t_2 &=& (a_1+b_1i) \cdot (a_2+b_2i) \\ &=& (a_1a_2-b_1b_2) + (a_1b_2+a_2b_1)i \\ &\Leftrightarrow& (a_1a_2-b_1b_2, a_1b_2+a_2b_1) \tag{2.9} \end{array} 除此之外,由欧拉公式易得极坐标系下的乘法运算结果。
\begin{array}{l} t_1 \cdot t_2 &=& r_1e^{i \theta_1} \cdot r_2e^{i \theta_2} \\ &=& r_1r_2e^{i(\theta_1+\theta_2)} \\ &\Leftrightarrow& (r_1r_2, \theta_1 + \theta_2) \tag{2.10} \end{array} 同理易得指数运算的结果如下。
\begin{array}{l} t^n &=& (re^{i \theta})^n \\ &=& r^ne^{i \cdot n \theta} \\ &\Leftrightarrow& (r^n, n \theta) \tag{2.11} \end{array} 单位复数根
单位复数是指
(r = 1)的复数,即落在复平面以坐标原点为圆心的单位圆上的复数。
单位复数根[6]的定义如下。
定义
n次单位复数根
是指n次幂
为1的复数,即
\omega^n = 1 , n \in \{1, 2, 3, \ldots \} \tag{3.1} 设
(\omega \Leftrightarrow (r, \theta),其中\theta \in [0, 2\pi)),则根据公式
((2.11))有
\omega^n \Leftrightarrow (r^n, n \theta) 故有
\left\{ \begin{array}{l} r^n = 1 \\ n \theta \mod 2\pi = 0 \end{array} \right. \Rightarrow \left\{ \begin{array}{l} r = 1 \\ n \theta = 2\pi \cdot k, k \in \{0, 1, 2, \ldots, n-1\} \end{array} \right. 即对于给定的次数n
,在
(\theta \in [0, 2\pi))上共有n
个单位复数根,它们是
(1, 2\pi \frac{0}{n}),(1, 2\pi \frac{1}{n}),(1, 2\pi \frac{2}{n}), \ldots ,(1, 2\pi \frac{n-1}{n}) 记为
\omega_n^{k} \Leftrightarrow (1, 2\pi \frac{k}{n}), k \in \{0, 1, 2, \ldots, n-1\} \tag{3.2} 注意,如果k可以取到n,那么此时有
(\omega_n^{n} \Leftrightarrow (1, 2\pi) \Leftrightarrow (1, 0)),它与
(\omega_n^0)重复了,都是复平面上x轴正方向上的单位根,即
\omega_n^n = \omega_n^0 \Leftrightarrow (1, 0) \tag{3.3} 所以说n次单位复根总共只有n个,剩下的就是重复的了,这n个根在复平面上如图所示。
单位复数根
单位复数根动态
性质
单位复数根有一些后面计算需要用到的性质。
相消引理
\omega_{dn}^{dk} = \omega_n^k,d \in Z_+ \tag{3.4} 由公式
((3.2))可证。
\begin{array}{l} \omega_{dn}^{dk} & \Leftrightarrow & (1, 2\pi \frac{dk}{dn}) \\ & = & (1, 2\pi \frac{k}{n}) & \Leftrightarrow & \omega_n^k \end{array} 折半引理
\omega_n^{k + \frac{n}{2}} = - \omega_n^k \text{, n为偶数} \tag{3.5} 同样由公式
((3.2))可得
\begin{array}{l} \omega_n^{k + \frac{n}{2}} &\Leftrightarrow& (1, 2\pi \frac{k+\frac{n}{2}}{n}) \\ &=& (1, 2\pi \frac{k}{n} + \pi) \\ &\Leftrightarrow& \cos(2\pi \frac{k}{n} + \pi) + i\sin(2\pi \frac{k}{n} + \pi) \\ &=& -\cos(2\pi \frac{k}{n}) - i\sin(2\pi \frac{k}{n}) \\ &\Leftrightarrow& -(1, 2\pi \frac{k}{n}) \\ &\Leftrightarrow& -\omega_n^k \end{array} 在复平面上表现为旋转了180度,即方向相反的两个向量。
求和引理
\sum_{k=0}^{n-1} \omega_n^k = 0 可由等比数列求和公式证明。
\begin{array}{l} \sum_{k=0}^{n-1} \omega_n^k &=& \sum_{k=0}^{n-1} e^{2\pi \frac{k}{n}i} \\ &=& \frac{e^{\frac{2\pi n i}{n}} - 1}{e^{\frac{2\pi i}{n}}-1} \\ &=& \frac{1 - 1}{e^{\frac{2\pi i}{n}}-1} &=& 0 \end{array} 在复平面上表现为所有的向量相加等于0。或者说所有向量终点组成的多边形的重心为坐标原点。
求和引理增强版
\[ \sum_{k=0}^{n-1} (\omega_n^k)^t = \left\{ \begin{array}{l} n & \text{, t \neq 0} \end{array} \right. \tag{3.6} \]
其中
(t \in Z)。此公式也可由等比数列求和证明。
当
(t=0)时,有
(\sum_{k=0}^{n-1} (\omega_n^k)^t = \sum_{k=0}^{n-1}1 = n)。
当
(t \neq 0)时,有
\sum_{k=0}^{n-1} (\omega_n^k)^t = (\omega_n^t)^0 + (\omega_n^t)^1 + (\omega_n^t)^2 + \cdots + (\omega_n^t)^{n-1} 显然这是一个首项为1,公比为
(\omega_n^t)的等比数列的前n项和。根据等比数列[9]的求和公式,可得
\sum_{k=0}^{n-1} (\omega_n^k)^t = \frac{(\omega_n^t)^n - 1}{\omega_n^t - 1} = \frac{(\omega_n^n)^t - 1}{\omega_n^t - 1} = \frac{1 - 1}{\omega_n^t - 1} = 0 得证。
离散傅立叶变换DFT
离散傅立叶变换(DFT, Discrete Fourier Transform)
[7]是指对于已知其系数表达式的n项多项式f(x)
,分别代入这n个n次单位复数根
(\omega_n^k),得到n个点值对,形成该多项式的点值表达式,即
f(x) \Leftrightarrow \vec{(x, f)} = \{(\omega_n^0, f(\omega_n^0)), (\omega_n^1, f(\omega_n^1)), \ldots, (\omega_n^k, f(\omega_n^k)), \ldots, (\omega_n^{n-1}, f(\omega_n^{n-1})) \} 由
((3.1)),若直接代入的话,时间复杂度为
(O(n^2))。但是我们可以利用单位复数根的性质,将问题转化为规模更小的子问题,使用分治策略处理问题,可以使复杂度降至
(O(nlogn)),这也是使用单位复数根的原因。
我们首先来考虑一个项数为
(n)的多项式f(x),其中
(n)为2的幂,即
(n = 2^k, k \in N)。其实这已经考虑了所有的多项式,因为项数不是2的幂的多项式,总可以添加系数为0的项,一直加到项数为2的幂为止即可。
f(x) = \sum_{i=0}^{i \lt n} a_i \cdot x^i = a_0 + a_1x^1 + a_2x^2 + \cdots + a_{n-1}x^{n-1} 将其各项按照次数的奇偶性分组,即
f(x) = \underbrace{(a_0 + a_2x^2 + \cdots + a_{n-2}x^{n-2})}_{\text{偶次数项}} + \underbrace{(a_1x + a_3x^3 + \cdots + a_{n-1}x^{n-1})}_{\text{奇次数项}} 提取奇次数项中的公因数
(x),可得两个关于
(x^2)的多项式,分别记为
(A_1(x^2))和
(A_2(x^2)),即
\[ f(x) = \underbrace{(a_0 + a_2x^2 + \cdots + a_{n-2}x^{n-2})}_{\text{偶次数项A_1(x^2)}} + x \cdot \underbrace{(a_1 + a_3x^2 + \cdots + a_{n-1}x^{n-2})}_{\text{奇次数项A_2(x^2)}} \]
则有 \[ \left\{ \begin{array}{l} A_1(x) = (a_0 + a_2x + a_4x^2 + \cdots + a_{n-2}x^{\frac{n-2}{2}}),\text{共\frac{n}{2}项} \\ A_2(x) = (a_1 + a_3x + a_5x^2 + \cdots + a_{n-1}x^{\frac{n-2}{2}}),\text{共\frac{n}{2}项} \end{array} \right. \Rightarrow \begin{array}{l} f(x) = A_1(x^2) + x \cdot A_2(x^2) \tag{3.7} \end{array} \]
下面我们开始代入n个单位复数根
(\omega_n^k, k \in \{0, 1, 2, \ldots, n-1\})。
我们将k的取值区间分为左右两个部分,每个部分
(\frac{n}{2})个元素来分别讨论。
左侧部分 对于左边的
(\frac{n}{2})个数,即
(k \in \{0, 1, 2, \ldots, \frac{n}{2} -1 \}),由公式
((3.7))可得
f(\omega_n^k) = A_1(\omega_n^{2k}) + \omega_n^k \cdot A_2(\omega_n^{2k}) 由相消引理
((3.4))可得
f(\omega_n^k) = A_1(\omega_{\frac{n}{2}}^{k}) + \omega_n^k \cdot A_2(\omega_{\frac{n}{2}}^{k}) \tag{3.8} 右侧部分 对于右边的
(\frac{n}{2})个数,我们将左边的
(k)值向右偏移
(\frac{n}{2}),即仍然取
(k \in \{0, 1, 2, \ldots, \frac{n}{2}-1\}),通过代入
(k + \frac{n}{2}),取遍右边的值,即
(k + \frac{n}{2} \in \{\frac{n}{2}, \frac{n}{2}+1, \ldots, n\})。
此时由公式
((3.7))有
f(\omega_n^{k+\frac{n}{2}}) = A_1(\omega_n^{2k+n}) + \omega_n^{k+\frac{n}{2}} \cdot A_2(\omega_n^{2k+n}) 由折半引理
((3.5))可得
(\omega_n^{k+\frac{n}{2}} = -\omega_n^k)。
同时有
(\omega_n^{2k+n} = \omega_n^{2k} \cdot \omega_n^n = \omega_n^{2k}),再由相消引理
((3.4))可得
(\omega_n^{2k+n} = \omega_{\frac{n}{2}}^{k})。
所以上式可转化为
f(\omega_n^{k+\frac{n}{2}}) = A_1(\omega_{\frac{n}{2}}^k) - \omega_n^k \cdot A_2(\omega_{\frac{n}{2}}^k) \tag{3.9} 综上所述,我们的
(k)值只需要取左侧的
(\frac{n}{2})个,即可通过
((3.8)(3.9))得到所有的
(n)个点值对。
\left\{ \begin{array}{l} f(\omega_n^k) &=& A_1(\omega_{\frac{n}{2}}^{k}) + \omega_n^k \cdot A_2(\omega_{\frac{n}{2}}^{k}) \\ f(\omega_n^{k+\frac{n}{2}}) &=& A_1(\omega_{\frac{n}{2}}^k) - \omega_n^k \cdot A_2(\omega_{\frac{n}{2}}^k) \end{array} \right. \tag{3.10} 与此同时,我们只需要求出项数为
(\frac{n}{2})的多项式
(A_1(x))的
(\frac{n}{2})个单位复数根
(\omega_{\frac{n}{2}}^k, k \in \{0, 1, 2, \ldots, \frac{n}{2}-1\})所表示的
(\frac{n}{2})个点值对即可,
(A_2(x))同理。
从上面这句描述可以发现,这是一个与原问题完全相同,只不过规模从
(n)变为
(\frac{n}{2})的一个子问题,所以我们完全可以使用分治的策略逐层划分为规则减半的子问题。
上述过程可以很容易转化成代码实现,一个有效的参考如下。
// complex
struct cp {
double r, i;
cp(double rr = 0, double ii = 0):r(rr), i(ii) {}
cp operator + (cp x) { return cp(r+x.r, i+x.i); }
cp operator - (cp x) { return cp(r-x.r, i-x.i); }
cp operator * (cp x) { return cp(r*x.r-i*x.i, r*x.i+i*x.r); }
};
/*
* @input a -- 系数组成的数组,长度为项数,a[i]表示第i项的系数。
* @input n -- 项数。
* @output y -- 单位复数根处的值,y[i]表示第i个单位复数根处的值f(w_n^i)。
*/
void dft(in const cp a[], in int n, out cp y[]) {
// 只有一个元素时为常数项,系数即为值
if(n == 1) {
y[0] = a[0];
return;
}
// 第一步,系数分左右
cp *ae = new cp[n>>1], *ao = new cp[n>>1];
for(int i = 0; i < n/2; ++i) {
ae[i] = a[i<<1]; // A_1
ao[i] = a[i<<1|1]; // A_2
}
// 第二步,分治策略递归求规则更小的子问题
cp *ye = new cp[n>>1], *yo = new cp[n>>1];
dft(ae, n>>1, ye); // 求n/2个单位复根在A_1(x)的值
dft(ao, n>>1, yo); // 求n/2个单位复根在A_2(x)的值
// 第三步,合并结果,公式(3.10)
cp wn = cp(cos(2*PI/n), sin(2*PI/n)); // 单位复数根的迭代增量
cp w(1, 0); // k = 0时的单位复根
for(int i = 0; i < n/2; ++i) {
y[i] = ye[i] + w*yo[i];
y[i + n/2] = ye[i] - w*yo[i];
w = w*wn; // 相当于单位复根w_n^k中的k++
}
delete ae; delete ao; delete ye; delete yo;
}
显然它的时间复杂度为
(T(n) = 2T(\frac{n}{2}) + O(n) = O(nlogn)),可通过主定理[8]求出。
离散傅立叶逆变换IDFT
离散傅立叶逆变换(IDFT, Inverse Discrete Fourier Transform)
是指对于已知点值表达式的多项式f(x)
,求其系数表达式
(\vec{a})的过程。
IDFT
同样可以使用单位复数根加速运算。
假设对于多项式
(f(x) \Leftrightarrow \vec{a}),
(\vec{a})是我们需要求的系数表达式,已知其n个单位复数根所对应的点值对
(\vec{(x, f)}),即
(x_i, f_i) = (\omega_n^i, f(\omega_n^i)),i \in \{0, 1, 2, \ldots, n-1\} 我们以这n个值
(\vec{f} = (f_0, f_1, \ldots, f_{n-1}))作为系数构造新的n项多项式g(x)
,即
f_i = f(\omega_n^i) = \sum_{j=0}^{n-1}a_j(\omega_n^i)^j 那么g(x)
的表达式为
g(x) = \sum_{i=0}^{n-1}f_ix^i 我们在g(x)
上取n个点计算其值,这n个点满足如下条件,即上面单位复数根的共轭复数
。
\omega_n^{-k}, k \in \{0, 1, 2, \ldots, n-1\} 对于任意
(k \in \{0, 1, 2, \ldots, n-1\}),我们计算其值具体是多少。
g(\omega_n^{-k}) = \sum_{i=0}^{n-1}f_i(\omega_n^{-k})^i 代入
(f_i),得
\begin{array}{l} g(\omega_n^{-k}) &=& \sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i \\ &=& \sum\limits_{j=0}^{n-1}(a_j\sum\limits_{i=0}^{n-1}(\omega_n^i)^j)(\omega_n^{-k})^i) \\ &=& \sum\limits_{j=0}^{n-1}(a_j\underbrace{\sum\limits_{i=0}^{n-1}(\omega_n^i)^{j-k}}_{G(j, k)}) \\ &=& a_0G(0, k) + a_1G(1, k) + \cdots + a_{n-1}G(n-1, k) \end{array} 设
(j-k = t),则可根据上面求和引理增强版
((3.6))得
\[ G(j, k) = \left\{ \begin{array}{l} n & \text{, j \neq k} \end{array} \right. \]
那么上面的g(x)
只有
(j = k)即
(G(k, k))项有值,即
g(\omega_n^{-k}) = a_0G(0, k) + a_1G(1, k) + \cdots + a_{n-1}G(n-1, k) = a_kG(k, k) = a_k \cdot n 那么对任意
(k \in \{0, 1, 2, \ldots, n-1\}),都有
a_k = \frac{g(\omega_n^{-k})}{n} 所以我们只需要求出来
(\omega_n^{-k})处的这n个g(x)
的值,便是f(x)
的系数表达式了,而求这n个g(x)
的值的方法,就跟DFT
一样使用分治策略即可,只不过把
(\omega_n^{-k})就可以了,证明过程也与DFT
一样,此处不作赘述。
总结一下流程。
1.已知f(x)
在单位复数根
(\omega_n^k)处的这n个值,将这n个值作为系数构造多项式g(x)
。
2.求g(x)
在
(\omega_n^{-k})这n个地方的值,便可得到f(x)
的系数表达式
(\vec{a})。
这个过程与DFT
只有两个差别,就是将
(\omega_n^k)变成
(\omega_n^{-k}),同时最终的结果除以
(n)即可。所以稍微改一下DFT
的代码即可实现。
/*
* @input a -- 单位复数根处的值组成的数组,长度为项数,a[i]表示第w_n^i处的的值。
* @input n -- 项数。
* @output y -- 系数*n,y[i]表示第i个系数*n,所以该结果要除以n才是真正的系数。
*/
void idft(in const cp a[], in int n, out cp y[]) {
// 只有一个元素时为常数项,系数即为值
if(n == 1) {
y[0] = a[0];
return;
}
// 第一步,系数分左右
cp *ae = new cp[n>>1], *ao = new cp[n>>1];
for(int i = 0; i < n/2; ++i) {
ae[i] = a[i<<1]; // A_1
ao[i] = a[i<<1|1]; // A_2
}
// 第二步,分治策略递归求规则更小的子问题
cp *ye = new cp[n>>1], *yo = new cp[n>>1];
idft(ae, n>>1, ye); // 求n/2个单位复根在A_1(x)的值
idft(ao, n>>1, yo); // 求n/2个单位复根在A_2(x)的值
// 第三步,合并结果,公式(3.10)
cp wn = cp(cos(-2*PI/n), sin(-2*PI/n)); // 与dft唯一的修改 k -> -k
cp w(1, 0); // k = 0时的单位复根
for(int i = 0; i < n/2; ++i) {
y[i] = ye[i] + w*yo[i];
y[i + n/2] = ye[i] - w*yo[i];
w = w*wn; // 相当于单位复根w_n^k中的k++
}
delete ae; delete ao; delete ye; delete yo;
}
可以使用额外的参数将DFT
和IDFT
整合起来。
/*
* @input a -- 系数组成的数组,长度为项数,a[i]表示第i项的系数。
* @input n -- 项数。
* @input d -- 1表示DFT,-1表示IDFT。
* @output y -- DFT时表示单位复数根处的值,y[i]表示第i个单位复数根处的值f(w_n^i)。
* @output y -- IDFT时表示系数*n,y[i]表示第i个系数*n,所以该结果要除以n才是真正的系数。
*/
void fft(in const cp a[], in int n, out cp y[], int d = 1) {
// 只有一个元素时为常数项,系数即为值
if(n == 1) {
y[0] = a[0];
return;
}
// 第一步,系数分左右
cp *ae = new cp[n>>1], *ao = new cp[n>>1];
for(int i = 0; i < n/2; ++i) {
ae[i] = a[i<<1]; // A_1
ao[i] = a[i<<1|1]; // A_2
}
// 第二步,分治策略递归求规则更小的子问题
cp *ye = new cp[n>>1], *yo = new cp[n>>1];
fft(ae, n>>1, ye); // 求n/2个单位复根在A_1(x)的值
fft(ao, n>>1, yo); // 求n/2个单位复根在A_2(x)的值
// 第三步,合并结果,公式(3.10)
cp wn = cp(cos(d*2*PI/n), sin(d*2*PI/n)); // 单位复数根的迭代增量
cp w(1, 0); // k = 0时的单位复根
for(int i = 0; i < n/2; ++i) {
y[i] = ye[i] + w*yo[i];
y[i + n/2] = ye[i] - w*yo[i];
w = w*wn; // 相当于单位复根w_n^k中的k++
}
delete ae; delete ao; delete ye; delete yo;
}
迭代版实现
我们可以将上面的分治策略改为迭代实现。
注意到每次实施分治策略的时候,对于每一层递归,我们都将奇数项分在一起,偶数项分在一起。这相当于考察这一项在这一层递归中所对应的二进制位是否为1,如果是1,则说明它是奇数,分到奇数堆里面去,反正分到偶数堆里面去。
举个例子,当
(n=8)时,对于一开始处于第三项的元素,它的二进制位是011
,那么第一层递归的时候,考察它最右边的一位的值是1,所以它会被分到奇数堆里去,即下标会往右放,假设它最终被放到的下标的二进制位是xyz
,那么x位必然为1,因为这一轮只有那4个偶数的x位才是0。
那么继续下一轮递归,此时考察它的中间的二进制位,这个二进制位也是1,那么它又被往右放了,且它最终被放到的下标的y位也是1。
继续下一轮,此时考察它的最左边的二进制位,这个二进制位是0,那么它也被往左放了,所以它最终下标的z位是0。
那么011
到了最后一轮,便会被放到下标为110
的位置,即它原来二进制位的反转结果。
如果我们确定了各个元素最终所在的下标,那么在每一层合并的时候都不会打乱上一层的下标顺序,就可以在原数组位置进行合并,如下代码所示。
// 二进制反转
int rev(int id, int n) {
int ret = 0;
for(int i = 0; (1<<i) < n; ++i) {
ret <<= 1;
if(id & (1<<i)) ret |= 1;
}
return ret;
}
// 另一种二进制反转
void rev(cp y[], int n) {
for(int i = 1, j = (n>>1); i < n-1; ++i) {
if(i < j) swap(y[i], y[j]); int k = (n>>1);
while(j >= k) { j -= k; k >>= 1; } if(j < k) j += k;
}
}
void FFT(cp a[], int n, cp y[], int d) {
// 二进制反转 -- 先将元素放到正确的下标
//for(int i = 0; i < n; ++i) y[i] = a[i]; rev(y, n); //另外一种二进制反转
for(int i = 0; i < n; ++i) y[rev(i, n)] = a[i];
// 从最底层开始迭代
for(int i = 1; (1<<i) <= n; ++i) { // 层数
int m = (1<<i); // 当前层元素个数
cp wn = cp(cos(d*2*PI/m), sin(d*2*PI/m));
for(int j = 0; j < n; j += m) { // 每块的起始下标
cp w = cp(1, 0);
for(int k = 0; k < (m>>1); ++k) { // 遍历每块的1/2个元素
// 在原数组位置进行合并 -- 蝴蝶操作
cp u = y[j+k];
cp t = w*y[j+k + (m>>1)];
y[j+k] = u + t;
y[j+k + (m>>1)] = u - t;
w = w*wn;
}
}
}
// IDFT直接获得系数
if(d == -1) for(int i = 0; i < n; ++i) y[i].r /= n, y[i].i /= n;
}
使用FFT实现多项式乘法
通过上面的fft,我们可以将插值
和求值
操作的复杂度都降到
(O(nlogn)),那么我们就可以使用fft来高效的实现多项式的乘法了。
// 返回大于x且离x最近的2的幂
int near(int x) {
int i = 0;
while(x > (1<<i)) ++i;
return 1<<i;
}
inline int dcmp(long double a) { if(a < -eps) return -1; return (a > eps); }
// 多项式乘法
void poly_mul(double a[], int na, double b[], int nb, out double c[], out int& nc) {
int i, j, n = na > nb ? na : nb;
n = 1 << ((int)ceil(log(2*n)/log(2)-eps)); // n = near(na+nb-1);
// f(x)和g(x)分别求值 -- O(nlogn)
cp *x = new cp[n], *ya = new cp[n], *yb = new cp[n], *yc = new cp[n];
for(i = 0; i < n; ++i) x[i] = cp(i < na? a[i] : 0, 0); FFT(x, n, ya, 1);
for(i = 0; i < n; ++i) x[i] = cp(i < nb? b[i] : 0, 0); FFT(x, n, yb, 1);
// 值相乘 -- O(n)
for(i = 0; i < n; ++i) yc[i] = ya[i]*yb[i];
// 插值 --得到最终结果 -- O(nlogn)
FFT(yc, n, x, -1);
for(i = 0; i < n; ++i) c[i] = x[i].r;
for(nc = n; nc > 0 && dcmp(c[nc-1]) == 0; --nc);
delete x; delete ya; delete yb; delete yc;
}
从代码可以看出,我们同样像文章开始那样通过求值(将系数转为点值)
、点值相乘
、插值(将点值转为系数)
,三个步骤,但复杂度更低的
(O(nlogn))优化了多项式的乘法问题。这也是FFT
最典型的应用之一了。
傅立叶级数
理解了离散傅立叶变换之后,我们可以继续更进一步理解更一般形式的傅立叶变换了。
接下来我们会先了解如何将周期函数通过傅立叶级数
展开,然后再了解傅立叶级数
是如何“进化”出傅立叶变换
的,最后再了解又是怎么“究级进化”成广义傅立叶变换
的。
在了解傅立叶级数之前,我们需要先了解一些信号统计和数学相关的前置知识。
周期信号和正弦波
周期信号广泛存在于我们的日常生活中,如简谐运动[10]、电磁波、引力波等。我们可以使用波动方程
来描述周期信号对应的波动。其中最经典的要数正弦波
[11],其波动方程如下。
f(t) = A \cdot \sin(\omega t + \varphi) \tag{4.1}
正弦波
A是振幅,是振动的幅度。
(\omega)是角频率(角速度),若频率为
(k),则有
(\omega = 2\pi k)。
(\varphi)是相位,表示波在起始时间处的状态。
周期
(T = \frac{1}{k} = \frac{2\pi}{\omega})我们可以线性组合多个不同参数的正弦波,以组成新的波动函数,从而形成新的波形。比如
f(t) = sin(t) + 2sin(3t) 这是一个自变量为时间
的函数,显然,它会随着时间的偏移而波动,同时具有周期性,但是可能有不同于正弦波的波形,所以我们说它是一个新的波动函数。
值得注意的地方是,上面的波动函数,其实是一个时域
上的函数,即它的自变量是时间
,它的图像的x轴是以时间为单位的。
那么我们是不是可以以频率
为自变量,或者说以角频率
(\omega)为自变量,统计这个波在各个频率上的振幅。比如上面的函数,它在
(\omega = 1)的地方的振幅为1,在
(\omega = 3)的地方的振幅为2。那么当我们建立以
(\omega)为x轴,以振幅为y轴的坐标系时,它的图像就变成了下图所示。
频域
那么这个图像所表示的函数解析式便是一个离散的函数,即
\[ F(\omega) = \left\{ \begin{array}{l} 1 & \text{, \omega = 3} \\ \end{array} \right. \]
如果
(\omega)的取值也是连续的,那么
(F(\omega))就变成了一个连续的函数,这个函数的作用域是频率,所以这个函数我们又称之为频域
上的函数。此时原来时域
上的函数
(f(t))也就变成了由无数个正弦波线性组合而成的了。
这就是喜闻乐见的时域
与频域
的概念。
频域与时域(图网侵删)
除此之外,还有一个值得思考的问题,是不是所有的周期信号(周期函数),都可以通过正弦波的线性组合来表示?或者更进一步,是不是所有的函数都可以通过正弦波的线性组合来表示?这个问题也许就是傅立叶当时思考的问题,我们后面再展开讨论。
级数
级数[12]是一个有穷或无穷序列的和。比如我们常见的等差数列和、等比数列和等。
s = \sum_{n=0}^{\infty}u_n = u_0 + u_1 + u_2 + \cdots 其中笔者认为比较常用且需要掌握的便是泰勒级数。
泰勒级数
泰勒级数[13]是指在数学上,对于一个在实数或复数a
邻域上,以实数作为变量或以复数作为变量的函数,并且是无穷可微的函数f(x),它的泰勒级数是以下这种形式的幂级数
\sum_{n=0}^{\infty} \frac{f^{(n)}(a)}{n!} (x-a)^n 其中
(f^{(n)}(a))表示
(f(x))在点a
处的n阶导数
。如果
(a = 0),则该级数又称麦克劳林级数
。
如
(e^x)的麦克劳林级数为
e^x = \sum_{n=0}^{\infty}\frac{x^n}{n!} = 1 + x + \frac{x^2}{2} + \frac{x^3}{3} + \cdots 如果
(f(x))在a
处
(n+1)次可导,那么对于这个区间的任意
(x),都有
\begin{array}{l} f(x) &=& \sum_{i=0}^{n} \frac{f^{(i)}(a)}{i!} (x-a)^i + R_n(x) \\ &=& f(a) + \frac{f'(a)}{1!}(x-a) + \frac{f^{(2)}a}{2!}(x-a)^2 + \cdots + \frac{f^{(n)}a}{n!}(x-a)^n + R_n(x) \end{array} 该多项式称为
(f(x))在a处的泰勒展开式
,其中
(R_n(x))是泰勒公式的余项,它是
((x-a)^n)的高阶无穷小。
可以看出,通过泰勒公式
,我们可以将
(f(x))以多项式的方式近似出来,这个多项式称为泰勒多项式
。
下面我们开始研究傅立叶级数
。
定义
对于满足一定条件的任意周期函数
(f(t)),设其周期为
(T),角频率为
(\omega = \frac{2\pi}{T}),则
(f(t))都可以使用一系列的正弦函数
展开,这里有一个比较炫酷的可视化效果网站[14]。
f(t) = A_0 + \sum_{n=1}^{\infty}A_n \sin(n\omega t + \varphi_n) \tag{9.1} 可以看出公式
((9.1))中
(t)为变量,
(\omega)已知,振幅
(A_n)和初相
(\varphi_n)是我们需要求的,但目前这种形式不太好求。
由三角恒等式两角和公式
(\sin(\alpha + \beta) = \sin\alpha\cos\beta + \cos\alpha\sin\beta),可得
\begin{array}{l} A_n \sin(n\omega t + \varphi_n) &=& \overbrace{A_n\sin(\varphi)}^{a_n}\cos(n\omega t) + \overbrace{A_n\cos(\varphi)}^{b_n}\sin(n\omega t) \\ &=& a_n\cos(n\omega t) + b_n\sin(n\omega t) \end{array} 同时设
(A_0 = \frac{a_0}{2})(为了结果与
(a_n)统一格式),故有
f(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty} [a_n\cos(n\omega t) + b_n\sin(n\omega t)] \tag{9.2} 这便是傅立叶级数
的一种形式,最终求出的结果如下。
\[ \left\{ \begin{array}{l} a_n = \frac{2}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}} f(t) \cdot \cos(n\omega t) dt & \text{, n \in \{1, 2, \ldots \}} \\ \end{array} \right. \tag{9.3} \]
公式
((9.2)(9.3))即为将
(f(t))使用傅立叶级数
展开的结果。
其中f(t)需要满足的条件如下。
f(t)
在一个周期内连续,或在一个周期内只有有限个第一类间断点
(左右极限都存在的间断点)。
f(t)
至多只有有限个极值点
。
上述条件又称狄利克雷充分条件
。当满足上述条件时,我们说f(t)
的傅立叶级数是收敛的,且有
若点t处f(t)
是连续的,则级数收敛于f(t)
。
若点t处f(t)
是间断的,则级数收敛于
(\frac{f(t-o) + f(t+o)}{2}),其中
(o)为无穷小。
下面来看上述结果是怎么推导出来的。
推导
观察公式
((9.2)),我们能取到的所有的
(\cos)函数值和
(\sin)函数值,由以下集合组成。
\{1, \cos nx, \sin nx | n = 1, 2, 3 \ldots \} = \{1, \cos x, \cos 2x, \ldots, \sin x, \sin 2x, \ldots \} 这个集合里的元素在区间
([-\pi, \pi])正交
,即两两乘积在这个区间的积分(专业术语叫做内积
)为0,如下所示。
\left\{ \begin{array}{l} &\int_{-\pi}^\pi \cos nx \mathrm{d}x = 0 &, n=1,2,3,\cdots\\ &\int_{-\pi}^\pi \sin nx \mathrm{d}x = 0 &, n=1,2,3,\cdots\\ &\int_{-\pi}^\pi \cos mx \cos nx \mathrm{d}x = 0 &, m\neq n, m, n=1,2,3,\cdots\\ &\int_{-\pi}^\pi \sin mx \sin nx \mathrm{d}x = 0 &, m\neq n, m, n=1,2,3,\cdots\\ &\int_{-\pi}^\pi \sin mx \cos nx \mathrm{d}x = 0 &, m, n=1,2,3,\cdots \end{array} \right. 可通过三角恒等式
证明,不作赘述。
我们称这个集合定义了一个坐标系三角函数系
,集合里的元素都是该坐标系的原始轴(类比于迪卡尔坐标系的xyz),我们称之为基
或基底函数
。那么我们求的
(a_n)和
(b_n),就是f(t)
在这个坐标系下的坐标,即f(t)
在各个基底
上的投影
。
下面开始正式求解。
第一步 积分
设
(u = \omega t,F(u) = f(t)),则由换元积分法
,可得在
(u \in [-\pi, \pi]),即
(t \in [-\frac{T}{2}, \frac{T}{2}])上的积分为
\int\limits_{-\pi}^{\pi} F(u) \mathrm{d}u = \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) (\omega t)' \mathrm{d}t = \omega \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) \mathrm{d}t 所以有
\begin{array}{l} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) \mathrm{d}t &=& \frac{1}{\omega} \int\limits_{-\pi}^{\pi}F(u) \mathrm{d}u \\ &=& \frac{1}{\omega} \int\limits_{-\pi}^{\pi} \frac{a_0}{2} \mathrm{d}u + \frac{1}{\omega} \int\limits_{-\pi}^{\pi}( \sum\limits_{n=1}^{\infty} [a_n\cos(nu) + b_n\sin(nu)] ) \mathrm{d}u \\ &=& \frac{1}{\omega} \int\limits_{-\pi}^{\pi} \frac{a_0}{2} \mathrm{d}u + \frac{1}{\omega} \sum\limits_{n=1}^{\infty}( a_n\underbrace{\int_{-\pi}^{\pi}\cos(nu) \mathrm{d}u}_{=0} + b_n\underbrace{\int_{-\pi}^{\pi}\sin(nu) \mathrm{d}u}_{= 0} ) \\ &=& \frac{1}{\omega} \int\limits_{-\pi}^{\pi} \frac{a_0}{2} \mathrm{d}u \\ &=& \frac{1}{\omega} \frac{a_0}{2} \cdot [u]|_{-\pi}^{\pi} = \frac{T}{2\pi} \frac{a_0}{2} \cdot 2\pi = \frac{T}{2} a_0 \end{array} 所以有
a_0 = \frac{2}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) \mathrm{d}t \tag{9.4} 第二步 乘
(\cos) 再积分
对公式
((9.2))两边同乘以
(\cos(k\omega t)),其中
(k = 1, 2, 3, \ldots),得
f(t)\cos(k\omega t) = \frac{a_0}{2}\cos(k\omega t) + \sum\limits_{n=1}^{\infty} [a_n\cos(n\omega t)\cos(k\omega t) + b_n\sin(n\omega t)\cos(k\omega t)] 设
(u = \omega t,F(u) = f(t)),则由换元积分法
,可得在
(u \in [-\pi, \pi]),即
(t \in [-\frac{T}{2}, \frac{T}{2}])上的积分为
\int\limits_{-\pi}^{\pi} F(u)\cos(ku) \mathrm{d}u = \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) \cos(k\omega t) (\omega t)' \mathrm{d}t = \omega \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) \cos(k\omega t) \mathrm{d}t 所以有
\begin{array}{l} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t)\cos(k\omega t) \mathrm{d}t &=& \frac{1}{\omega} \int\limits_{-\pi}^{\pi}F(t)\cos(ku) \mathrm{d}u \\ &=& \frac{1}{\omega} \frac{a_0}{2}\underbrace{\int_{-\pi}^{\pi} \cos(ku) \mathrm{d}u}_{=0} + \frac{1}{\omega}\int\limits_{-\pi}^{\pi}( \sum\limits_{n=1}^{\infty} [a_n\cos(nu)\cos(ku) + b_n\sin(nu)\cos(ku)] ) \mathrm{d}u \\ &=& \frac{1}{\omega} \sum\limits_{n=1}^{\infty}( a_n\underbrace{\int_{-\pi}^{\pi}\cos(nu)\cos(ku) \mathrm{d}u}_{n \neq k 时=0} + b_n\underbrace{\int_{-\pi}^{\pi}\sin(nu)\cos(ku) \mathrm{d}u}_{= 0} ) \\ &=& \frac{1}{\omega} a_k\int\limits_{-\pi}^{\pi}\cos(ku)\cos(ku) \mathrm{d}u \\ &=& \frac{1}{\omega} \frac{a_k}{2}\int\limits_{-\pi}^{\pi}(1+\cos(2ku)) \mathrm{d}u \\ &=& \frac{1}{\omega} \frac{a_k}{2}( \int\limits_{-\pi}^{\pi}1 \mathrm{d}u + \underbrace{\int_{-\pi}^{\pi}\cos(2ku) \mathrm{d}u}_{=0} ) \\ &=& \frac{1}{\omega} \frac{a_k}{2} \cdot [u]|_{-\pi}^{\pi} = \frac{T}{2\pi}\frac{a_k}{2} \cdot 2\pi = \frac{T}{2} a_k \end{array} 所以有
a_k = \frac{2}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t)\cos(k\omega t) \mathrm{d}t, k \in \{1, 2, 3 \ldots \} 至此,我们求出了
(a_1, a_2, a_3, \ldots)的结果如上式所示,同时由第一步结果
(a_0)也满足上式,故有
a_n = \frac{2}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t)\cos(n\omega t) \mathrm{d}t, n \in \{0, 1, 2, 3 \ldots \} \tag{9.5} 第三步 乘
(\sin) 再积分
第三步与第二步类似,只不过将
(\cos)变成
(\sin)即可。
对公式
((9.2))两边同乘以
(\sin(k\omega t)),其中
(k = 1, 2, 3, \ldots),得
f(t)\sin(k\omega t) = \frac{a_0}{2}\sin(k\omega t) + \sum\limits_{n=1}^{\infty} [a_n\cos(n\omega t)\sin(k\omega t) + b_n\sin(n\omega t)\sin(k\omega t)] 设
(u = \omega t,F(u) = f(t)),则由换元积分法
,可得在
(u \in [-\pi, \pi]),即
(t \in [-\frac{T}{2}, \frac{T}{2}])上的积分为
\int\limits_{-\pi}^{\pi} F(u)\sin(ku) \mathrm{d}u = \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) \sin(k\omega t) (\omega t)' \mathrm{d}t = \omega \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) \sin(k\omega t) \mathrm{d}t 所以有
\begin{array}{l} \int\limits_{\frac{T}{2}}^{\frac{T}{2}}f(t)\sin(k\omega t) \mathrm{d}t &=& \frac{1}{\omega} \int\limits_{-\pi}^{\pi}F(t)\sin(ku) \mathrm{d}u \\ &=& \frac{1}{\omega} \frac{a_0}{2}\underbrace{\int_{-\pi}^{\pi} \sin(ku) \mathrm{d}u}_{=0} + \frac{1}{\omega} \int\limits_{-\pi}^{\pi}( \sum\limits_{n=1}^{\infty} [a_n\cos(nu)\sin(ku) + b_n\sin(nu)\sin(ku)] ) \mathrm{d}u \\ &=& \frac{1}{\omega} \sum\limits_{n=1}^{\infty}( a_n\underbrace{\int_{-\pi}^{\pi}\cos(nu)\sin(ku) \mathrm{d}u}_{=0} + b_n\underbrace{\int_{-\pi}^{\pi}\sin(nu)\sin(ku) \mathrm{d}u}_{n \neq k 时=0} ) \\ &=& \frac{1}{\omega} b_k\int\limits_{-\pi}^{\pi}\sin(ku)\sin(ku) \mathrm{d}u \\ &=& \frac{1}{\omega} \frac{b_k}{2}\int\limits_{-\pi}^{\pi}(1-\cos(2ku)) \mathrm{d}u \\ &=& \frac{1}{\omega} \frac{b_k}{2}( \int\limits_{-\pi}^{\pi}1 \mathrm{d}u - \underbrace{\int_{-\pi}^{\pi}\cos(2ku) \mathrm{d}u}_{=0} ) \\ &=& \frac{1}{\omega} \frac{b_k}{2} \cdot [u]|_{-\pi}^{\pi} = \frac{T}{2\pi} \frac{b_k}{2} \cdot 2\pi = \frac{T}{2} b_k \end{array} 所以有
b_k = \frac{2}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t)\sin(k\omega t) \mathrm{d}t, k \in \{1, 2, 3 \ldots \} 至此,我们求出了
(b_1, b_2, b_3, \ldots)的结果如上式所示,即
b_n = \frac{2}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t)\sin(n\omega t) \mathrm{d}t, n \in \{1, 2, 3 \ldots \} \tag{9.6} 至此,
((9.5)(9.6))得到了与一开始
((9.3))一致的结果,这便是求解傅立叶级数的整个过程。
复数形式
根据欧拉公式
((2.5)) 即
(e^{i\theta} = \cos{\theta} + i \cdot \sin{\theta})可得
\left\{ \begin{array}{l} \cos{\theta} = \frac{e^{i\theta} + e^{-i\theta}}{2} \\ \sin{\theta} = \frac{e^{i\theta} - e^{-i\theta}}{2i} = \frac{-i(e^{i\theta}-e^{-i\theta})}{2}\\ \end{array} \right. 代入
((9.2))得
\begin{array}{l} f(t) &=& \frac{a_0}{2} + \sum\limits_{n=1}^{\infty} [a_n\cos(n\omega t) + b_n\sin(n\omega t)] \\ &=& \frac{a_0}{2} + \sum\limits_{n=1}^{\infty} [ a_n\frac{e^{in\omega t} + e^{-in\omega t}}{2} + b_n\frac{-i(e^{in\omega t} - e^{-in\omega t})}{2} ] \\ &=& \underbrace{\frac{a_0}{2}}_{c_0} + \sum\limits_{n=1}^{\infty} [ \underbrace{\frac{a_n-ib_n}{2}}_{c_n} e^{in\omega t} + \underbrace{\frac{a_n+ibn}{2}}_{c_{-n}} e^{i(-n)\omega t} ] \\ &=& c_0 + \sum\limits_{n=1}^{\infty} [ c_n e^{in\omega t} + c_{-n} e^{i(-n)\omega t} ] \\ \end{array} \tag{9.7} 所以有
\[ \begin{array}{l} c_0 &=& \frac{a_0}{2} &=& \frac{1}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) \mathrm{d}t \\ c_n &=& \frac{a_n-ib_n}{2} &=& \frac{1}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) [\cos((-n)\omega t) +i\sin((-n)\omega t)] \mathrm{d}t = \frac{1}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) e^{i(-n)\omega t} \mathrm{d}t & \text{, n = 1, 2, 3, \ldots} \\ \end{array} \]
显然这三个公式可以统一为以下形式。
\[ c_n = \frac{1}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) e^{-in\omega t} \mathrm{d}t \text{, n = \ldots, -3, -2, -1, 0, 1, 2, 3, \ldots} \tag{9.8} \]
代入
((9.7)),得
f(t) = \sum_{n = -\infty}^{\infty}c_ne^{in\omega t} \tag{9.9} 这便是傅立叶级数
的复数形式了,十分优雅。
显然我们如果已知
(f(t)),便可求出
(c_n)。如果已知
(c_n),也可求出
(f(t)),这也是时域和频域的相互转换,十分优雅。
不管是
((9.3))中的
(a_n)或
(b_n),还是
((9.8))中的
(c_n),如果我们把它们都看作是自变量为角频率
(\omega),因变量为振幅的频域
函数的话,那么每一个
(n)都可以确认一个角频率
(\omega_n),显然这个频域函数的定义域显然是离散的,即有
c_n = F(\omega_n) \tag{9.10} 其中
(\omega_n)的取值是离散的。
所说如果f(x)
的T逐渐增大的话,它的角频率
(\omega)就会逐渐减小,那么相邻的
(\omega_n)之间的距离就会减小,也就是说它的定义域的取值间隔越来越小,它的定义域在x轴上就会变得越来越密集。
考虑一种极端情况,当
(T \rightarrow \infty)时,
(\omega_n)之间的间距
(\rightarrow)无穷小,此时相邻的
(\omega_n)之间的距离为无穷小,即此时定义域将覆盖整个x轴,即
(F(\omega_n))变成了自变量为频率的连续函数,且此时f(x)
也不再是周期函数。而傅立叶变换
考虑的就是这种情况。
傅立叶变换
定义
由上面的分析可以看出,傅立叶级数只对周期函数有效。而傅立叶变换[16]就是将傅立叶级数推广到非周期函数的方法。
观察复数形式的傅立叶级数
((9.8)(9.9)),我们上面也讨论了,当
(T \rightarrow \infty)时,
(\omega_n)之间的间距
(\rightarrow)无穷小,所以我们从
((9.10))得到了一个__连续__的频域函数,且
(f(t))也并非周期函数了。
此时将有
\left\{ \begin{array}{l} f(t) &=& \color{}{} & \int_{-\infty}^{\infty}F(\omega)e^{i\omega t} \mathrm{d}\omega & \text{, 傅立叶逆变换(频域->时域),对频域(角频率)微元的积分}\\ F(\omega) &=& \frac{1}{2\pi} & \int_{-\infty}^{\infty}f(t)e^{-i\omega t} \mathrm{d}t & \text{, 傅立叶变换(时域->频域),对时域(时间)微元的积分}\\ \end{array} \right. \tag{10.1} 证明
观察复数形式的傅立叶级数
((9.8)(9.9)),我们将
((9.8))代入
((9.9)),并取
(T \rightarrow \infty),得
f(t) = \lim_{T \rightarrow \infty} \sum_{n = -\infty}^{\infty} [ \frac{1}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(m) e^{-in\omega m} \mathrm{d}m ]e^{in\omega t} 这里我们将
(c_n)中的积分微元
(t)换一个符号
(m),是为了与外层的
(f(t))的自变量
(t)作下区分。
显然,对于
(c_n)积分项,外层的
(e^in\omega t)是一个常数,可乘进积分项的每一项中,即
f(t) = \lim_{T \rightarrow \infty} \sum_{n = -\infty}^{\infty} [ \frac{1}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(m) e^{in\omega (t-m)} \mathrm{d}m ] \tag{10.2} 设
(\omega_n = n\omega),即随着n的不同而在定义域坐标轴上选取的不同的位置。那么定义域上相邻两个自变量的距离为
\begin{array}{l} \Delta \omega &=& \omega_n - \omega_{n-1} = \omega = \frac{2\pi}{T} \\ &\Rightarrow& T = \frac{2\pi}{\Delta\omega} \end{array} 代入
((10.2)),得
\begin{array}{l} f(t) &=& \lim\limits_{\Delta\omega \rightarrow 0} \sum\limits_{n = -\infty}^{\infty} [ \frac{1}{2\pi} \Delta\omega \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(m) e^{i\omega_n (t-m)} \mathrm{d}m] \\ &=& \lim\limits_{\Delta\omega \rightarrow 0} \sum\limits_{n = -\infty}^{\infty} [ \frac{1}{2\pi} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(m) e^{i\omega_n (t-m)} \mathrm{d}m ] \Delta\omega \end{array} 这是一个
(\Delta\omega \rightarrow 0)的极限,显然这等价于一个微元为
(\mathrm{d}\omega)的积分,即
\begin{array}{l} f(t) &=& \int\limits_{n = -\infty}^{\infty} [ \frac{1}{2\pi} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(m) e^{i\omega (t-m)} \mathrm{d}m ] \mathrm{d}\omega \\ &=& \int\limits_{n = -\infty}^{\infty} [ \frac{1}{2\pi} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(m) e^{-i\omega m)} \mathrm{d}m ] e^{i\omega t} \mathrm{d}\omega \end{array} 所以有
\left\{ \begin{array}{l} F(\omega) &=& \frac{1}{2\pi}\int\limits_{-\infty}^{\infty}f(t)e^{-i\omega t} \mathrm{d}t \\ f(t) &=& \int\limits_{-\infty}^{\infty}F(\omega)e^{i\omega t} \mathrm{d}\omega \end{array} \right. 得证。
可以看到,离散的
(c_n)变成连关于
(\omega)的连续函数
(F(\omega))。当我们传入一个频率作为自变量时,就会得到关于这个频率的信息(振幅、初相等),所以又称
(F(\omega))为时域函数
(f(t))的频谱
。可以直观的看出
(f(t))的频率分布情况。
傅立叶变换
将傅立叶级数
扩展到非周期函数,这使我们具备了将任意非周期函数
(f(t))从连续的时域函数
得到连续的频域函数
,或者从连续的频域函数
得到连续的时域函数
的能力。
但是对于周期函数
,我们却无法使用傅立叶变换
来展开,而是只能用傅立叶级数
展开。
所以我们需找寻找一种可以统一傅立叶级数
和傅立叶变换
二者的公式,使傅立叶变换
既可以应用于非周期函数
,又可以应用于周期函数
。这便是广义傅立叶变换
。
广义傅立叶变换
在了解广义傅立叶变换
之前,我们需要先来看一个特殊的函数,即狄拉克δ函数
[17],小名冲激函数
。
冲激函数
定义 狄拉克δ函数
的直观定义如下。
\[ \delta(x) = \left\{ \begin{array}{l} + \infty & \text{, x \neq 0} \\ \tag{11.1} \end{array} \right. \]
同时满足
\int_{-\infty}^{\infty} \delta(x) \mathrm{d}x = 1 \tag{11.2} 筛选性质
(\delta)函数有一个很重要的性质,即筛选性质即
,即对于
(\forall)定义域为
((-\infty, \infty))的函数
(f(x)),有
\int_{-\infty}^{\infty}\delta(x - x_0) f(x) \mathrm{d}x = f(x_0) \tag{11.3} 根据定义
((11.2))可知,当且仅当
(x = x_0)时
(f(x))项不为0,故
(f(x))项可提出,得
\begin{array}{l} \int_{-\infty}^{\infty}\delta(x - x_0) f(x) \mathrm{d}x &=& f(x_0)\overbrace{\int_{-\infty}^{\infty}\delta(x - x_0) \mathrm{d}x}^{= 1} \\ &=& f(x_0) \end{array} 之所以叫筛选性质
,是因为它可以将
(f(x))在
(x_0)处的值筛选出来。
定义
对于
(\forall)弦函数
(f(t) = e^{i\omega_0 t}),可以通过傅立叶变换
得到其频域函数
(F(\omega))如下。
F(\omega) = \frac{1}{2\pi}\int_{-\infty}^{\infty}e^{i(\omega_0-\omega)t} \mathrm{d}t = \delta(\omega_0 - \omega) \tag{11.4} 同时频域函数
(F(\omega))还可通过傅立叶逆变换
变换回时域函数
(f(t)),即
f(t) = \int_{-\infty}^{\infty}\delta(\omega_0 - \omega)e^{i\omega t} \mathrm{d}\omega = e^{i \omega_0 t} \tag{11.5} 注意其中
(\omega_0)为原函数的角频率,
(\omega)为最终求得的频域函数的自变量。
与此同时,我们知道,对于任意周期函数,都可通过傅立叶级数
展开成弦函数的线程组合,即公式
((9.9)) f(t) = \sum_{n = -\infty}^{\infty}c_ne^{in\omega_0 t} \tag{9.9} 那么对于任意周期函数
(f(t)),代入公式
((11.4)),都可得到其频域函数
(F(\omega))如下。
F(\omega) = \sum_{n = -\infty}^{\infty}c_n \delta(n\omega_0 - \omega) \tag{11.6} 同理,对于任意周期函数的频域函数
(F(\omega)),都可通过公式
((11.5))对每一项进行傅立叶逆变换
,也都可以重新得到其时域函数
(f(t)),即公式
((9.9))。
公式
((11.4)(11.5)(11.6))即为广义傅立叶变换
,它使得周期函数也可以通过傅立叶变换
展开为
(\delta)函数的线程组合。且只要证明
((11.4)(11.5))成立,即可推导出
((11.6))。下面来证明。
证明
(\delta)函数为一个非周期函数,我们对它进行傅立叶变换,即对时域函数
(f(t) = \delta(t))进行傅立叶变换,可得其频域函数
(F(\omega))为
F(\omega) = \frac{1}{2\pi}\int_{-\infty}^{\infty}\delta(t)e^{-i\omega t} \mathrm{d}t 根据其筛选性质
((11.3))可得
\begin{array}{l} F(\omega) &=& \frac{1}{2\pi}\int_{-\infty}^{\infty}\delta(t)e^{-i\omega t} \mathrm{d}t \\ &=& \left. \frac{1}{2\pi} e^{-i\omega t} \right \vert_{t = 0} \\ &=& \frac{1}{2\pi} \end{array} 即
(\delta(t))的频谱为一个常函数
(F(\omega) = \frac{1}{2\pi})。我们对该常函数重新进行傅立叶逆变换
,便可得到原来的时域函数
(\delta(t)),即
\delta(t) = \int_{-\infty}^{\infty}\frac{1}{2\pi}e^{i\omega t} \mathrm{d}\omega 由于
(\omega)和
(t)的定义域都是
((-\infty, \infty)),所以对于上式,我们换个符号表示,即将
(\omega)用
(t)表示,将
(t)用
(\omega)表示,该式仍然成立,即
\[ \delta(\omega) = \frac{1}{2\pi}\int_{-\infty}^{\infty}e^{i\omega t} \mathrm{d}t \text{, \omega \in (-\infty, \infty)} \tag{11.7} \]
由
((11.7))取
(\omega = \omega_0 - \omega)显然公式
((11.4))成立。
由
(\delta)函数的筛选性质
((11.3)),易证公式
((11.5))也成立,即
\begin{array}{l} f(t) &=& \int_{-\infty}^{\infty}\delta(\omega_0 - \omega)e^{i\omega t} \mathrm{d}\omega \\ &=& \left. e^{i\omega t} \right \vert_{\omega = \omega_0} \\ &=& e^{i\omega_0 t} \end{array} 得证。
后记
至此,我们了解了离散傅立叶变换
、傅立叶级数
、傅立叶变换
以及广义傅立叶
变换的一些并不十分严谨的定义和证明。正如前文所述,本文的侧重点是对概念和思想的熟悉,以及推导过程中的重要步骤和结果,顺带着科普一些基础但很有趣的数学知识,所以有选择地忽略了很多细枝末节(如边界条件、级数的敛散性及相关证明、一些性质和应用等),读者可自行脑(goo)补(gle)。
个人认为最重要的是了解时域
和频域
概念及其互相转换的思想。我们工作中也许会碰到一些在时域中处理不了或者不好处理的问题,转换到频域上处理可能能找到一条更好的道路。
我们在傅立叶级数
的推导过程中,有提到过内积
的概念,以及三角函数系
的概念。其实这里的三角函数系
在泛函分析
中是一个希尔伯特空间
,即完备的内积空间
,它相当于将三维的欧几里德空间推广到n维乃至无限维但又不失完备性。而傅立叶变换
相于在n维的希尔伯特空间
的上向各个坐标轴(基底)的投影。
而广义傅立叶变换
与泛函数分析
就更密切了,因为
(\delta)函数本身就不是一个严格意义上的函数。它可以定义为一个分布
,即广义函数
,即给定一个
(f(x)),
(\delta(f(x)))可以得到某个值,即
(\delta函数)是测试函数
(f(x))的一个线性泛函。
(\delta)函数可以使用测度
的概念严谨定义,它其实是n维希尔伯特空间
的一个测度,一个函数相对于
(\delta)函数的积分便可认为是相对于这个测度的勒贝格积分
。
以上两段没有什么实际参考意义,纯粹是为了显得专业
而胡说八道,限于笔者自身水平,本文的绝大部分公式的推导和证明都是初等的微积分而已,应该很容易理解。不过如果读者掌握过《测度论》《泛函分析》等更专业的知识,也许就可以从另外一个更加高级的抽象层次理解和描述问题,无奈本人非数学专业出身,也没读过研-_-||,目前来说很难再从更高一层的视角作抽象和总结,希望以后可以找个时间填一下这些坑,此刻只能深表遗憾了。
参考
[1]多元线性方程组的解
[2]秦九韶算法
[3]高斯消元法
[4]复数
[5]欧拉公式
[6]单位复数根
[7]离散傅立叶变换
[8]主定理
[9]等比数列
[10]简谐运动
[11]正弦波
[12]级数
[13]泰勒级数
[14]Fourier-Visualizations
[15]An Interactive Introduction to Fourier Transforms
[16]傅立叶变换
[17]狄拉克δ函数