【科普】音频开发与音乐制作【第三篇-时频分析上篇】

在上一篇中了解了音频在计算机中的表示方式后，我们将在这一篇中讨论一个非常重要的概念：傅里叶变换。

注：因为本文是科普角度为主，并不是所有内容都附有严谨数学推导，文章有失误或不准确之处敬请谅解~

时域与频域

在本篇之前，我们一直是以时域来讨论音频的。时间域（time domain）描述了振幅-时间的关系，自变量是时间，因变量是振幅。一个C调属和弦（g1-b1-d1）的时域图看起来是这样的：

这里为了简化情况方便理解，三个音都是用正弦波来演奏的。此时，眼前的波形不再是一个简单的正弦波了，而是三个音相加的结果。

计算机中存储的PCM采样都是时域图。然而在思考音乐时，我们更多想的是包含了哪些频率（不同音高的音）——乐谱上标的都是这些。为了做更多有用的操作，我们需要从频率角度思考问题。

1807年，约瑟夫·傅里叶提出了一个相当重要的观点：任何周期函数都可用正弦函数级数表示。任何函数，只要是周期性的，都可以看做是一系列正弦函数之和。这就是傅里叶级数（Fourier Series）。当然，后来的研究证明周期函数只有满足狄利克雷条件（Dirichlet Conditions）时才能展开成傅里叶级数，即一个周期内第一类间断点和极大极小值数目有限且函数可积。

频率域（frequency domain）告诉我们当前波形的频率成分，换句话说，波形是由那些频率的正弦波组成的。这个属和弦的频域图长这样：

这下就相当明了，可以看到该波形包含了振幅为0.4的391.995Hz、振幅为0.4的493.883Hz正弦波和一个振幅为0.3的587.330Hz正弦波。由此可知，以上时域图的函数表达式如下：

$ f(x)=0.3\sin(391.995\cdot2\pi x)+0.3\sin(493.883\cdot2\pi x)+0.4\sin(587.330\cdot2\pi x) $

要确定一个正弦函数$y=A\sin(\omega x+\phi )+k$，我们需要知道它的振幅$A$、周期（频率）$\omega$、相位偏移$\phi$和垂直方向偏移$k$。音频波形始终以x轴为中心，所以k为零，不用考虑。频域图描述了每个正弦波的频率$\omega$和振幅$A$。

相位偏移$\phi$则由相位谱描述。要想准确的还原波形的话需要它，此外别的一些领域也会用到。但在这里就先不细说了。

读到这里，你可能觉得频域图是个很陌生的东西……其实并不。下面这个也是一种频域图，只不过以观赏目的为主：

再举个傅里叶级数的例子。方波看起来和圆润的正弦波一点也不像，但实际上，如果我们将一个正弦波和它的奇次谐波叠加，就像这样：$\sin{x}+\frac{\sin3x}{3}+\frac{\sin5x}{5}+\frac{\sin7x}{7}+\frac{\sin9x}{9}+…$

或者写作：

$f(x)=\sum_{n=1}^{\infty }\frac{\sin((2n-1)x)}{(2n-1)}$

随着叠加谐波数量的增加，图像看起来会越来越像方波。

借用一下维基百科上的演示图：

很神奇有木有。

一个有趣的现象：当有不连续点的周期函数（比如上面的方波，每周期有两次跳变）进行傅里叶展开后，使用有限项进行合成时（比如上面图像用了一百项），不连续点附近会出现突起。随着级数的增多，突起的大小趋近于总跳变值的大约9%。这被称作吉布斯现象（Gibbs phenomenon）。

傅里叶变换

现在的问题是，如何根据函数的时域采样数据得到函数的频域分布？换句话说，将采样序列看做一个离散函数f(x)，如何知道f(x)里面各频率正弦波的强度？

关于傅里叶变换，3b1b的一期视频中有比较可视化直观的解释；但我还是建议从函数正交角度推导一下傅里叶变换的原理。无内鬼，来点数学。放心，不需要特别多高等数学知识，高中数学水平外加一点微积分常识就足够理解了。相信我，高中生也能懂。

正交函数

假设有两个向量$\vec{v}=(v_1,v_2)$ ，$\vec{u}=(u_1,u_2)$ ，它们之间可以进行点乘（Dot Product）运算：

$\vec{v}\cdot \vec{u}=v_1u_1+v_2u_2$

或者:

$\vec{v}\cdot \vec{u}=\left | \vec{v} \right | \left | \vec{u} \right | \cos\left ( \theta \right )$

后者是点乘的几何含义，$\theta$是两向量间的夹角。我们先不考虑几何含义，只考虑第一种表达方式。

可以看到点乘其实就是把向量的每一项乘起来，然后相加得到一个标量。上面是二维向量的情况。但不管多少维的向量都可以这么操作。对于向量$\vec{v}=(v_1,v_2, … ,v_n)$ ，$\vec{u}=(u_1,u_2, … ,u_n)$ ，点乘运算为：

$\vec{v}\cdot \vec{u}=v_1u_1+v_2u_2+…+v_nu_n$

如果两个向量点乘得到的内积为0，我们称这两个向量为正交（Orthogonal）的。正交的向量在空间上垂直，比如$(0,3)$和$(5,0)$是相互垂直的，$0\times3+5\times0=0$。

接下来可能有点怪。我们把函数看做向量，函数图像上的每一个点看做向量的一项。换句话说，每个连续函数都是无限维的向量，离散函数为无限维或有限维的向量。对于定义域为$[a,b]$的连续函数$f(x)$，我们可以定义如下向量：

$\vec{f}=\lim_{\Delta x \to 0} (f(a),f(a+\Delta x),f(a+2\Delta x),…,f(a+n\Delta x),…,f(b))$

如果将声音采样序列a[]看做一个离散函数a[x]=$f(x)$，那么$\Delta x$为采样间隔。其实计算机中的数组也是一种向量么，你看C++ STL库中的动态数组容器就叫vector。

看做向量的函数一样可以进行点乘运算，只要它们的定义域相同。运算方式和普通向量一样：把对应项相乘，再加起来。对于定义域都为$[a,b]$函数$f(x)$和$g(x)$，它们的点乘结果为：

$\lim_{\Delta x \to 0}f(a)g(a)+f(a+\Delta x)g(a+\Delta x)+f(a+2\Delta x)g(a+2\Delta x)+…+f(a+n\Delta x)g(a+n\Delta x)+…+f(b)g(b)$

事实上也就是求如下积分：

$\int_{a}^{b} f(x)g(x)\mathrm{d}x$

因此，如果说定义域为$[a,b]$的实函数f(x)与g(x)正交，那么

$\int_{a}^{b} f(x)g(x)\mathrm{d}x=0$

三角函数系的正交性质

傅里叶变换利用了三角函数系的正交性质。

首先，正弦函数和余弦函数在整数周期内积分为0。在下面的图像上，蓝绿区域面积显然相等，因此互相抵消掉了：

即$\int_{a}^{a+n\frac{2\pi}{\omega }} \sin{\omega t} \mathrm{d}t=0$，n为正整数。这里$\frac{2\pi}{\omega }$是正弦函数的周期。

那正弦函数$\sin m_1\omega t$和$\sin m_2\omega t$在相同的整数周期定义域内是否正交呢？（$m_1$，$m_2$为正整数，$m_1\ne m_2$）我们已经知道，这两个函数的内积为：

$\int_{b}^{a} \sin m_1\omega t\sin m_2\omega t\mathrm{d}x$

根据三角函数的积化和差公式（prosthaphaeresis formulas）：

$\sin \alpha \sin \beta=-\frac{1}{2}(-2\sin \alpha \sin \beta)$

$=-\frac{1}{2}\left[\cos\alpha\cos\beta-\sin \alpha \sin \beta-\right(\cos\alpha\cos\beta+\sin \alpha \sin \beta)]$

$=-\frac{1}{2}\left[\cos(\alpha+\beta)-\cos(\alpha-\beta)\right]$

得到：

$\int_{b}^{a} \sin m_1\omega t\sin m_2\omega t\mathrm{d}x=\frac{1}{2}\left [ \int_{a}^{a+n\frac{2\pi}{\omega }} \cos{(m_1+m_2)\omega t} \mathrm{d}t +\int_{a}^{a+n\frac{2\pi}{\omega }} \cos{(m_1-m_2)\omega t} \mathrm{d}t \right ]$

我们已经知道$m_1$，$m_2$为正整数，那么$(m_1+m_2)$，$(m_1-m_2)$也是整数。所以$\int_{a}^{a+n\frac{2\pi}{\omega }} \cos{(m_1+m_2)\omega t} \mathrm{d}t$和$\int_{a}^{a+n\frac{2\pi}{\omega }} \cos{(m_1-m_2)\omega t} \mathrm{d}t$也都为0。

因此，$\sin m_1\omega t$和$\sin m_2\omega t$在整数周期的相同定义域内正交。

但是稍等一下！这里我们讨论的是$m_1\ne m_2$，也就是两个不同的三角函数是否正交。如果$m_1 = m_2$，两个三角函数相同呢？

将$m_1 = m_2$作为条件带入到$\frac{1}{2}\left [ \int_{a}^{a+n\frac{2\pi}{\omega }} \cos{(m_1+m_2)\omega t} \mathrm{d}t +\int_{a}^{a+n\frac{2\pi}{\omega }} \cos{(m_1-m_2)\omega t} \mathrm{d}t \right ]$中，我们可以得到：

$\frac{1}{2}\left [ \int_{a}^{a+n\frac{2\pi}{\omega }} \cos{(m_1+m_2)\omega t} \mathrm{d}t +\int_{a}^{a+n\frac{2\pi}{\omega }} \cos{(0)\omega t} \mathrm{d}t \right ]$

当然啦，我们已经知道$\int_{a}^{a+n\frac{2\pi}{\omega }} \cos{(m_1+m_2)\omega t} \mathrm{d}t=0$。但是$\cos 0 = 1$，很明显对1做积分结果可就不是0了。$\frac{1}{2}\int_{a}^{a+n\frac{2\pi}{\omega }} \cos{(0)\omega t} \mathrm{d}t=n\frac{\pi}{\omega }$诶！

如果我们给正弦波加上振幅$A$，变成$A\cos(m_1-m_2)\omega t$，那么对其积分的结果也就变成了$An\frac{\pi}\omega$。换句话说，正比于弦波振幅和持续时间。

欧拉公式

现在，我们已经推导出了三角函数（准确来说是正弦函数）的一个重要性质：在整数周期内，不同的三角函数内积为0（正交），相同的三角函数内积正比于正弦波振幅和持续时间。

如果设一段波形的时域为定义域$[a,b]$周期函数$f(x)$（定义域在这里可以看做采样范围），在符合此前推导所用条件时，其频域函数为$F(x)=\frac{1}{b-a}\int_{a}^{b}f(t) \cos{(xt)} \mathrm{d}t$。这种从时域到频域的变换也就是傅里叶变换（Fourier transform）。

实际的傅里叶变换中一般不使用三角函数，而使用欧拉公式（Euler’s formula）中的$e^{ix}$，并将得到的结果取实部或者取模。

傅里叶变换得到的虚部也是有用的。虚部结合实部可以计算得到频率的相位谱。相位的傅里叶公式由于没那么常用，篇幅原因就不放了。

欧拉公式：$e^{ix}=\cos x + i \sin x$

连续的傅里叶变换公式为：