librosa | 梅尔谱图最幽默的解释

小白进来！

---

写在前面

Medium的一篇文章，特别幽默地介绍了梅尔谱图，快来一起轻松学习吧！

正文开始

原文：《Getting to Know the Mel Spectrogram》
作者：Dalya Gartzman
时间：2019.08.20

了解梅尔频谱图

阅读这篇短文，如果你想像Neo（电影《黑客帝国》男主角，即上图中的人物）一样，并了解所有关于Mel Spectrogram的信息！（可能不会涵盖到梅尔谱图的全部内容，但至少有一点）

有关 tl;dr 和完整代码，请转到此处：达利亚格/CodeSnippetsForPosterity。

几天前发生在我脑海中的真实对话

• 我：嗨，梅尔频谱图，我可以叫你梅尔吗？

• 梅尔：当然。

• 我：谢谢。所以梅尔，当我们第一次见面时，你对我来说是一个谜。

• 梅尔：真的吗？怎么样？

• 我：你由两个概念组成，它们的全部目的是让人类能够理解抽象的概念——梅尔尺度和频谱图——但你自己对我这个人类来说很难理解。

• 梅尔：这种片面的演讲有意义吗？

• 我：你知道是什么让我更加困扰吗？我通过小道消息听说你是DSP（数字信号处理）中的嗡嗡声，但我在网上发现关于你的直观信息很少。

• 梅尔：我应该为你难过吗？

• 我：不管怎样，我不想让你被误解，所以我决定写你。

• 梅尔：哎呀！这个主意不错。希望现在有更多的人能得到我。

• 我：很高兴，我的朋友。我认为我们可以讨论你的核心元素是什么，然后在python上使用librosa包展示一些不错的技巧。

• 梅尔：噢噢，太好了！我爱librosa！它可以用一行代码生成我！

• 我：太棒了！让我们在这篇文章中使用这首美丽的鲸鱼歌曲作为我们的玩具示例！你觉得怎么样？

  鲸鱼歌声链接：https://www.youtube.com/watch?v=W5Trznre92c

• 梅尔：你知道你在自言自语吗？

频谱图

可视化声音是一个迷幻的概念。有一些令人着迷的方法可以做到这一点，还有更多的数学方法，我们将在这篇文章中探讨。

当我们谈论声音时，我们通常谈论不同压力强度的一系列振动，因此将声音可视化意味着可视化电波。

import librosa
import librosa.display

filename = 'Haunting_song_of_humpback_whales.wav'
y, sr = librosa.load(filename)
# trim silent edges
whale_song, _ = librosa.effects.trim(y)
librosa.display.waveplot(whale_song, sr=sr);

但这只是这首复杂而丰富的鲸鱼之歌的二维表现！

声音的另一个数学表示是傅里叶变换。傅里叶变换无需过多细节（观看此教育视频以获得可理解的解释），傅里叶变换是一种函数，可将时域中的信号作为输入，并将其分解输出为频率。

让我们以一个短时间窗口为例，看看我们从应用傅里叶变换中得到了什么。

import numpy as np
import matplotlib.pyplot as plt

n_fft = 2048
D = np.abs(librosa.stft(whale_song[:n_fft], n_fft=n_fft, hop_length=n_fft+1))
plt.plot(D);

现在，让我们获取完整的鲸鱼歌曲，将其分离到时间窗口，并在每个时间窗口上应用傅里叶变换。

hop_length = 512
D = np.abs(librosa.stft(whale_song, n_fft=n_fft,  hop_length=hop_length))
librosa.display.specshow(D, sr=sr, x_axis='time', y_axis='linear');
plt.colorbar();

哇，这里看不出来，不是吗？这是因为人类听到的大多数声音都集中在非常小的频率和幅度范围内。

注：作者的这张插图意在说明：狗狗对50,000Hz最敏感。
狗狗的听阈范围为60~65000Hz
人的听阈范围为20~20kHz

让我们再做一个小的调整 - 将y轴（频率）转换为对数刻度，并将“颜色”轴（振幅）转换为分贝，这有点像振幅的对数刻度。

DB = librosa.amplitude_to_db(D, ref=np.max)
librosa.display.specshow(DB, sr=sr, hop_length=hop_length, x_axis='time', y_axis='log');
plt.colorbar(format='%+2.0f dB');

这就是我们所说的频谱图！

梅尔量表

让我们暂时忘记所有这些可爱的可视化和数学。从数学上讲，梅尔标度是频率尺度的一些非线性变换的结果。这个Mel音阶的构造使得在Mel音阶上彼此之间距离相等的声音对人类的听觉来说，距离也是相等的。

与Hz标度相反，500和1000 Hz之间的差异是显而易见的，而7500和8000 Hz之间的差异几乎不明显。

幸运的是，有人为我们计算了这个非线性变换，我们需要做的就是使用librosa的适当命令来应用它。

n_mels = 128
mel = librosa.filters.mel(sr=sr, n_fft=n_fft, n_mels=n_mels)

是的。就是这样。
但这给了我们什么呢？
它将 Hz 标度划分为条柱，并使用重叠的三角形滤波器将每个箱转换为 Mel 尺度中的相应箱。

plt.figure(figsize=(15, 4));

plt.subplot(1, 3, 1);
librosa.display.specshow(mel, sr=sr, hop_length=hop_length, x_axis='linear');
plt.ylabel('Mel filter');
plt.colorbar();
plt.title('1. Our filter bank for converting from Hz to mels.');

plt.subplot(1, 3, 2);
mel_10 = librosa.filters.mel(sr=sr, n_fft=n_fft, n_mels=10)
librosa.display.specshow(mel_10, sr=sr, hop_length=hop_length, x_axis='linear');
plt.ylabel('Mel filter');
plt.colorbar();
plt.title('2. Easier to see what is happening with only 10 mels.');

plt.subplot(1, 3, 3);
idxs_to_plot = [0, 9, 49, 99, 127]
for i in idxs_to_plot:
    plt.plot(mel[i]);
plt.legend(labels=[f'{i+1}' for i in idxs_to_plot]);
plt.title('3. Plotting some triangular filters separately.');

plt.tight_layout();

现在这给了我们什么？
现在，我们可以采用一个时间窗口的振幅，计算点积来执行变换，并在这个新的`mel频率尺度下获得声音的可视化。

plt.plot(D[:, 1]);
plt.plot(mel.dot(D[:, 1]));
plt.legend(labels=['Hz', 'mel']);
plt.title('One sampled window for example, before and after converting to mel.');

嗯……好像还没有多大意义？
好吧，让我们把它全部总结起来，看看我们得到了什么。

梅尔频谱图

我们现在知道什么是三维频谱图，也知道什么是梅尔尺度，所以，令人惊讶的是，梅尔频谱图是一个以梅尔尺度为y轴的频谱图。

这就是你用一行代码生成Mel Spectrogram的方法，并且只用3行代码就可以很好地显示它：

S = librosa.feature.melspectrogram(whale_song, sr=sr, n_fft=n_fft, hop_length=hop_length, n_mels=n_mels)
S_DB = librosa.power_to_db(S, ref=np.max)
librosa.display.specshow(S_DB, sr=sr, hop_length=hop_length, x_axis='time', y_axis='mel');
plt.colorbar(format='%+2.0f dB');

回顾

Mel Spectrogram 是以下步骤的结果：

1. 与窗口分开：使用一定大小的窗口对输入进行采样，每次制作一定大小的帧移以对下一个窗口进行采样。n_fft=2048 hop_length=512
2. 为每个窗口计算FFT（快速傅里叶变换），实现从时域变换到频域。
3. 生成 Mel 标度：获取整个频谱，并将其分离成均匀分布的频率。
   我们所说的均匀间隔是什么意思？不是通过频率维度上的距离，而是人耳听到的距离。n_mels=128
4. 生成频谱图：对于每个窗口，将信号的大小分解为其分量，对应于mel尺度中的频率。

# Sanity check that indeed we understood the underlying pipeline
S = librosa.feature.melspectrogram(whale_song, sr=sr, n_fft=n_fft, hop_length=hop_length, n_mels=n_mels)
fft_windows = librosa.stft(whale_song, n_fft=n_fft, hop_length=hop_length)
magnitude = np.abs(fft_windows)**2
mel = librosa.filters.mel(sr=sr, n_fft=n_fft, n_mels=n_mels)

assert (mel.dot(magnitude) == S).all()

未完待续…

现在我们知道了Mel Spectrogram，我们打算用它做什么呢？

好吧，那是另一篇文章…

---

（完）