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

上篇讲完了傅里叶变换的基本原理后，本篇将会讲一些傅里叶变换在实际应用中的相关算法。

注：本文写作时间比较早，其中关于快速傅里叶变换的部分可能不是完全准确/易于理解。如果读者有兴趣进行代码实现更推荐阅读这篇文章。

快速傅里叶变换（FFT）

当在计算机中实现离散傅里叶变换DFT时，算法的时间复杂度就不得不纳入考量了。

根据上一篇中的DFT公式，我们不难发现，针对每一个频率求权重（振幅）的过程都需要$N$次复数乘法和$(N-1)$次复数加法。当求解$N$个频率时，合起来就是$O(n^2)$的时间复杂度。

当采样率来到40kHz以上后，可不是每台计算机都撑得住……我们需要对算法进行一些优化。

以下只考虑序列长度$N$为2的幂的情况。首先为了简化公式，我们记$W_N=e^{-j2\pi /N}$，此时DFT公式可写为

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

如果不理解的话，$W_N^k$的实部等于$\cos 2\pi kn/N$，虚部等于$i\sin 2\pi kn/N$。$W_N$具有以下性质：

周期性：$W_N^k=W_N^{k+tN}$（k,t为整数），也就是三角函数的周期性；$\sin 2\pi x=\sin 2\pi (x+k)$么。
对称性：$W_N^{k+N/2}=-W_N^k$。从三角函数的角度想就是$\sin x=-\sin (x+\pi)$
缩放性：$W_N^k=W_{N/2}^{k/2}$。因为代进去化简一下压根就是一回事么。

然后按照n的值把奇偶项分开。$X[k]=G[k]+H[k]$，偶数项$G[k]=\sum_{n为偶数}^{N-2}x[n]W_N^{kn}$，奇数项$H[k]=\sum_{n为奇数}^{N-1}x[n]W_N^{kn}$。

$X[k+N/2]$和$X[k]$有什么关系呢？

$X[k+N/2]=\sum_{n为偶数}^{N-2}x[n]W_N^{nk+nN/2}$。注意，n为偶数，所以$n/2$必然为整数。

利用周期性可得

$\sum_{n为偶数}^{N-2}x[n]W_N^{nk+nN/2}=\sum_{n为偶数}^{N-2}x[n]W_N^nk=G[k]$

奇数项同理：

$X[k+N/2]=\sum_{n为奇数}^{N-1}x[n]W_N^{nk+nN/2}$

$=\sum_{n为奇数}^{N-1}x[n]W_N^{nk+(n-1)N/2+N/2}$

因为n为奇数，$(n-1)/2$为整数。利用周期性得到：

$=\sum_{n为奇数}^{N-1}x[n]W_N^{nk+N/2}$

利用对称性得到：

$=-\sum_{n为奇数}^{N-1}x[n]W_N^{nk}=-H[k]$

因此：

$X[k+N/2]=G[k+N/2]+H[k+N/2]=G[k]-H[k]$

也就是说，我们只需要计算前半个序列的$G[k]$和$H[k]$就足够了，后半个序列可以复用前半个序列的计算结果。好耶！计算量减半了。

利用缩放性，我们还可以把$G[k]$和$H[k]$进行以下化简：

偶数项：

$G[k]=\sum_{n=0}^{N/2-1}x[2n]W_N^{2kn}=\sum_{n=0}^{N/2-1}x[2n]W_{N/2}^{kn}$

奇数项：

$G[k]=\sum_{n=0}^{N/2-1}x[2n+1]W_N^{(2n+1)k}$

$=W_N^k\sum_{n=0}^{N/2-1}x[2n+1]W_N^{2nk}$

$=W_N^k\sum_{n=0}^{N/2-1}x[2n+1]W_{N/2}^{nk}$

注意到了吗？偶数项和奇数项各自又是一个标准的DFT，可以再用上述方法再次减半计算量。现在算法的时间复杂度降低到了$o(n\log n)$。FFT在代码上用递归就可以实现。当然，DFT是一个对性能相当敏感的应用，而递归的大量函数调用比较影响速度，所以代码实现中经常会展开成蝶形计算等非递归方法。这样的DFT被称作快速傅里叶变换（Fast Fourier Transform, FFT）

短时傅里叶变换（STFT）

傅里叶变换有一个小问题：它是针对一整段采样的，求出的频域是一整段采样的平均值。现实中的音频显然不是这样：每一秒都有不同的音发声和停止发声，声音在不断变大变小。傅里叶变换只能求出一张固定的整体频域图，但真实的声音有着不断变化的频域。

短时傅里叶变换（Short-time Fourier Transform, STFT）假设音频的频域在很短的时间内是不变的。STFT将音频切成一片片帧（frame）后对每个帧分别进行FFT，这样就能得到以帧长度为时间分辨率的频域-时间变化图。时频图大概长这样：

由于截取的帧大小不一定是周期的整数倍，帧两端会出现频谱泄露（spectrum leakage）问题，导致求出来的频谱包含不该有的频率。解决方式一般是给帧加窗，即乘以一个窗函数（window function）。常用的有矩形窗，Hann窗和Hamming窗等。例如，信号可以长这样：

加如下窗：

加窗后的信号如下：

由于加窗后靠近帧两端的信号幅度变得比较小，处理中容易被忽略，所以帧与帧之间一般会有一定程度的重叠来保证所有频率都被均匀分析。

不同的窗有不同的频谱特性，一般使用时根据频率准确度和振幅准确度的需求来选择。对50.5Hz（频率分辨率为1Hz）的信号分别施加矩形窗（红色）、汉宁窗（绿色）和平顶窗（蓝色），用对数显示幅值，加窗后的结果如下图所示。

不确定性原理

傅里叶变换在分析周期函数时，考虑到截断效应等影响，分析的采样时间越长，分析出的频谱精确度越高。但这和我们的时域分析是矛盾的；分析的单个采样越长，时间分辨率越低。

信号的不确定性原理，也就是量子力学里的海森堡测不准原理。它指出：信号的时宽和频宽不可能同时任意地窄，这也就意味着，不可能通过任何方法上的改进同时得到完全精确的信号时间、频率信息。当使用短时窗进行傅里叶变换，信号的时间信息变得精确，同时频率信息变得不精确；当使用长时窗进行傅里叶变换，信号的频率信息变得精确，同时时间分布信息变得不精确，因为你得到的频率信息是分布在整个时窗中的。
@HDTT咸鱼

例如，对于如图的正弦波信号：