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

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

受时间所限,本文中的大部分内容作者并没有亲自代码实现过……所以基本都是基于理论的讲解,如有错误疏漏望谅解。本篇还会涉及一点点的量子力学。是的,我听个歌,混个音,还真能和量子力学扯上关系。

快速傅里叶变换(FFT)

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

也有用j替换i的,正负号在这里只是旋转方向的区别,不影响计算。

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

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

以下只考虑$2^N$个采样点的情况。首先为了简化公式,我们记$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$具有以下性质:

  1. 周期性:$W_N^k=W_N^{k+tN}$(k,t为整数),也就是三角函数的周期性;$\sin 2\pi x=\sin 2\pi (x+k)$么。
  2. 对称性:$W_N^{k+N/2}=-W_N^k$。从三角函数的角度想就是$\sin x=-\sin (x+\pi)$
  3. 缩放性:$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)的信号分别施加矩形窗(红色)、汉宁窗(绿色)和平顶窗(蓝色),用对数显示幅值,加窗后的结果如下图所示。

矩形窗,Hann窗和平顶窗

不确定性原理

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

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

@HDTT咸鱼

例如,对于如图的正弦波信号:

来源:知乎@J Pan

如果用较小的窗对其进行时频分析,可以看到时域分辨率高,频域分辨率低。

来源:知乎@J Pan

如果用较大的窗对其进行时频分析,频域分辨率倒是变高了,但是时域分辨率却变差了:

来源:知乎@J Pan

还记得量子力学中的德布罗意波吗?根据量子力学的基本假设,波函数描述了粒子出现在空间各处的概率,即粒子位置是一个空间上的波。同时,粒子的动量正比于波的频率。

当我们用傅里叶变换来从一个求另一个的时候,问题出现了:动量波函数和位置波函数可以互为傅里叶变换得到(共轭),就像我们的时域信号和频域信号一样。提高一个的分辨率,必然会导致另一个的分辨率降低。

海森堡不确定性原理,其实就是傅里叶变换中不确定性原理的一个例子。

(写到这里查资料的时候,发现3b1b对此也出过一期视频。那就把链接也附上吧。

有时,为了在时间分辨率和频率分辨率间合理取舍,变换时会对不同频率施加不同的窗大小,来做到尽量好的效果。此外还有小波变换(wavelet)等改进方法,这里就不详细讨论了。

时频分析在音频处理中非常有用。后面我们会用时频处理相关知识制作一个简单的虚拟歌手引擎

引用文献

[1] https://www.zhihu.com/question/264560305/answer/302819915

[2] https://zhuanlan.zhihu.com/p/24318554

[3] https://zhuanlan.zhihu.com/p/41548628

暂无评论

发送评论 编辑评论


|´・ω・)ノ
ヾ(≧∇≦*)ゝ
(☆ω☆)
(╯‵□′)╯︵┴─┴
 ̄﹃ ̄
(/ω\)
∠( ᐛ 」∠)_
(๑•̀ㅁ•́ฅ)
→_→
୧(๑•̀⌄•́๑)૭
٩(ˊᗜˋ*)و
(ノ°ο°)ノ
(´இ皿இ`)
⌇●﹏●⌇
(ฅ´ω`ฅ)
(╯°A°)╯︵○○○
φ( ̄∇ ̄o)
ヾ(´・ ・`。)ノ"
( ง ᵒ̌皿ᵒ̌)ง⁼³₌₃
(ó﹏ò。)
Σ(っ °Д °;)っ
( ,,´・ω・)ノ"(´っω・`。)
╮(╯▽╰)╭
o(*////▽////*)q
>﹏<
( ๑´•ω•) "(ㆆᴗㆆ)
😂
😀
😅
😊
🙂
🙃
😌
😍
😘
😜
😝
😏
😒
🙄
😳
😡
😔
😫
😱
😭
💩
👻
🙌
🖕
👍
👫
👬
👭
🌚
🌝
🙈
💊
😶
🙏
🍦
🍉
😣
Source: github.com/k4yt3x/flowerhd
洛天依九周年
颜文字
Emoji
小恐龙
花!
上一篇
下一篇