整段频谱为什么找不到声音发生的时间:短时傅里叶变换
两段声音都先后出现一个低音和一个高音,只把顺序交换,整段频率结果就可能看起来完全一样。要让程序找回“哪个声音先出现”,需要把录音分成互相重叠的小段,每到一个位置就重新计算一次。
导读: 两段声音都先后出现一个低音和一个高音,只把顺序交换,整段频率结果就可能看起来完全一样。要让程序找回“哪个声音先出现”,需要把录音分成互相重叠的小段,每到一个位置就重新计算一次。
第 14 课已经把一整段录音换成了一条曲线:横轴是频率,也就是一秒重复多少次,纵轴是每个频率有多强。这条曲线能看出小提琴、钢琴、萨克斯和噪声(也就是沙沙作响、没有固定高低的那种声音)各自由哪些成分组成,却指不出钢琴是在哪一刻敲下的。
本课要做的实验是:先故意造出两段“整段结果相同、声音顺序相反”的信号,再检查移动小段以后,程序能不能找回先后顺序。完整代码在第 15 课脚本,运行环境仍沿用课程总纲中的一次性配置。
实验:交换低音和高音的顺序,整段频谱还能分辨吗
先生成两段 1 秒声音。两段都包含每秒重复 300 次的低音和每秒重复 900 次的高音:
- 第一段先播放 300 Hz,再播放 900 Hz。
- 第二段正好反过来,先播放 900 Hz,再播放 300 Hz。
这里的 Hz 读作赫兹,就是一秒重复多少次。第 02 课解释过它与听到的高低有什么关系。
脚本先对两段完整声音各算一次幅度谱,也就是第 12 课那条「每个频率有多强」的曲线。两条曲线的最大差只有 4.547e-13。这个数已经接近浮点计算留下的舍入误差,意味着只看整段幅度谱,程序无法判断谁先谁后。
接着,脚本每次只看 800 个样本,每隔 400 个样本向右移动一次。第一段开头最强的是 300 Hz,结尾最强的是 900 Hz;第二段则从 900 Hz 变成 300 Hz。顺序重新出现了。
图 1:中间的整段幅度谱只有相同的 300 Hz 和 900 Hz 峰。下方逐帧结果才显示:一段从低到高,另一段从高到低。
这里要补一个边界,免得把话说过头:信息并没有凭空消失。只要保留完整的复数结果——也就是第 11 课那种一个数同时记住「有多强」和「从哪里起步」——DFT 仍然能把原声音一个数不差地重建回来。
真正的问题有两层。第一,每个 DFT 系数概括的都是整段录音,它身上没有一个属于自己的时间坐标。第二,我们最常看的幅度谱只取了模,把相位(也就是那个「从哪里起步」)整个丢掉了。两件事加起来,整段幅度谱就答不出「这个成分在第几秒」。
一次只看一小段,时间位置就留下来了
最直接的办法,是把录音切成许多短帧,对每一帧分别计算 FFT:
- 从录音里选出一帧。
- 给这一帧乘上一个两端低、中间高的窗,让切口平缓下来。
- 对这一帧计算 FFT,得到一列频率结果。
- 向右移动一小段距离,重复计算。
- 把所有列按先后顺序排在一起。
这个过程叫短时傅里叶变换,英文是 Short-Time Fourier Transform,简称 STFT。这里的“短时”不是换了一种傅里叶变换,而是让同一个计算不断处理不同时间位置的小段。
每次移动的距离叫帧移,下面记作 H。帧移比帧长小时,相邻两帧会共享一部分样本,这就是重叠。
把录音切成短帧、每帧乘一个窗,这两步第 06 课都做过。那一课解释了为什么它们要放在一起,本课要做的是把那套操作接到 FFT 后面。
每一帧会得到一列数字。把这些列按先后顺序并排成一张二维数字表,这张表就叫矩阵。
图 2:四步走一遍——取一帧、给它乘上窗、算一次 FFT、把得到的这一列排进矩阵。帧移 400 只有帧长 800 的一半,所以下一帧会和这一帧共用一半样本。
窗必须和帧一样长。原因很具体:程序要让帧里的每个样本都乘上一个窗值;两边数量不同,就无法逐项配对。本课的公共函数因此让窗长、帧长和 FFT 长度使用同一个 frame_length。
DFT 公式怎样多出一个时间位置
只处理一段长度为 N 的数据时,DFT 可以写成:
STFT 只改了三处:给输出增加“第几帧”,让输入从这一帧的起点开始,再给这一帧乘窗:
逐项读就不复杂了:
m表示现在计算第几帧。mH表示这帧从原录音的第几个样本开始。n表示帧内部的第几个样本。w[n]是这个样本位置对应的窗值。k仍然表示正在检查第几个频率位置。
当 m 从 0 变成 1、2、3 时,窗依次向右移动;每个位置都得到一整列 k。公式通常写成 S[m,k],本课程代码则把频率放在行、时间放在列,所以数组形状按“频率格数,帧数”排列。
10000 个样本为什么得到 501 × 19
现在给程序三个数:录音共有 10000 个样本,每帧取 1000 个样本,每次前进 500 个样本。这里只保留完整帧,不在两端补值。
1000 点实数帧只需要保留非负频率,因此行数是:
第一帧从 0 开始,最后一帧不能伸出第 10000 个样本,所以列数是:
最终形状就是 (501, 19):501 行频率,19 列时间。
图 3:上半部分沿录音摆出 19 个窗口;下半部分把每个窗口的 501 个非负频率结果竖成一列。两种计数在右侧合成 501 × 19。
公共函数的核心计算只有三行:
frames = frame(y, frame_length, hop_length)
windowed = frames * np.hanning(frame_length)[None, :]
S = np.fft.rfft(windowed, axis=1).T
脚本还用完全相同的 Hann 窗和 center=False 调用 librosa。两边形状都是 (501, 19),把 9519 个复数逐个相减,最大差是 0.000e+00。所以呢?这不是“看起来差不多”,而是每一个复数系数逐位相同——手写的这三行和库做的是同一件事。
center=False 表示只从原录音中取完整帧。有些库默认会在两端补值,因此列数可能不同;第 16 课真正调用绘图接口时会专门核对这个默认行为。
帧长决定时间和频率哪边看得更细
每个频率格之间相差“采样率除以帧长”。本课程使用 22050 Hz 采样率时,脚本得到:
| 帧长 | 一帧覆盖时间 | 相邻频率格间隔 |
|---|---|---|
| 512 | 23.22 ms | 43.07 Hz |
| 1024 | 46.44 ms | 21.53 Hz |
| 2048 | 92.88 ms | 10.77 Hz |
| 4096 | 185.76 ms | 5.38 Hz |
| 8192 | 371.52 ms | 2.69 Hz |
帧越长,一次看到的振动圈数越多,相邻频率越容易分开;代价是这一列结果覆盖了更长时间,很难指出短促声音究竟落在哪个瞬间。帧越短则相反:时间范围更窄,但频率格更疏。
下面两张图分析同一段从低到高的声音,并在 0.75 秒加入一次短促敲击。16 ms 窗的频率间隔是 62.50 Hz,64 ms 窗缩小到 15.62 Hz。
图 4:看 0.75 秒附近的竖亮线——短窗让它集中在更窄的时间范围;再看斜向亮带——长窗让它在频率方向更细。两边无法同时达到最细。
这种选择叫时间与频率的取舍。定位点击声、爆破音等短事件时,通常更在意时间;区分两个持续而接近的音高时,通常更在意频率。帧长没有脱离任务的“最佳值”。
帧移只决定隔多远再看一次
保持帧长 1000 不变,脚本只缩小帧移:
| 帧移 | 相邻帧重叠 | 完整帧数 |
|---|---|---|
| 500 | 50% | 19 |
| 250 | 75% | 37 |
| 125 | 87.5% | 73 |
帧移越小,程序在时间轴上取点越密,图看起来也更连续;代价是列数和计算量增加。它不会改变每一帧覆盖多长时间,因此不能消除帧长造成的取舍。把帧移减半,不等于突然获得了短一倍的观察窗口。
Hann 窗把每一帧的两端收下去
本课使用的 Hann 窗中间接近 1,两端是 0。1000 点窗的实测值是:左端 0.0,中心 0.999998,右端 0.0。
逐样本相乘以后,帧中央几乎保持原样,两端逐渐压低。这样做能减轻硬切口造成的频率泄漏,但不会免费消除所有误差;窗也会改变幅度和频率主瓣的宽度。这里沿用第 06 课的窗函数解释,不再重讲一次。
复数矩阵怎样变成一张声谱图
STFT 的每个位置仍然是复数。若只想用颜色显示强弱,可以先取模再平方:
得到的 Y 是功率矩阵。把时间帧放在横向、频率放在纵向、功率大小换成颜色,就得到功率声谱图。图 4 正是用这套数值画出来的:亮处表示该时间和频率附近更强,暗处表示更弱。
这张图暂时只负责解释“矩阵怎样变成颜色”。下一课才会把真实音阶和三段音乐依次交给 librosa,比较线性频率、相对 dB 和对数频率显示,并检查矩阵方向与补边方式。
这一课得到了什么,下一课接着讲什么
整段幅度谱把整段录音汇总成一条曲线,所以能回答“有哪些频率”,却不能直接回答“它们什么时候出现”。STFT 用移动短窗反复计算 FFT:每个窗口产生一列频率结果,所有列按时间并排成复数矩阵。
矩阵有多少行,由每帧包含多少个非负频率格决定;有多少列,由录音长度、帧长、帧移和补边方式决定。帧长控制时间与频率的取舍,帧移控制时间取点密度,窗函数控制每一帧的切口。
下一课要把这张复数矩阵真正画成可比较的声谱图。那时最容易出错的不再是公式,而是“画的是模还是功率”“颜色是否换成 dB”“纵轴按普通频率还是对数频率”以及“库是否在两端补了值”。