Valerio-机器学习与音频信号处理笔记-全-

Valerio 机器学习与音频信号处理笔记（全）

001：课程概述 🎵

在本节课中，我们将对《面向机器学习的音频信号处理》系列课程进行概述。我们将了解开设本课程的原因、将要涵盖的核心内容、目标受众以及学习本课程所需的预备知识。

许多深度学习工程师在开发图像相关的应用时，拥有丰富的资源来指导如何处理图像数据以适配模型。然而，对于音频数据，情况则不同。音频数据及其在机器学习应用中的使用方法，常常笼罩着一层迷雾。这正是本系列课程旨在解决的问题。

当我们讨论音频AI应用时，可以将其分为两个主要阶段。第一个阶段是模型的开发与评估，这在我的另一个系列课程《Python音频深度学习》中已有涵盖。第二个阶段则是准备原始的音频数据，使其能够作为模型的输入。虽然我在之前的课程中已经涉及了一些音频预处理和特征提取的内容，但显然这还不够深入，因此有了本系列课程。

课程应用场景 🎯

那么，音频数字信号处理具体在哪些机器学习（特别是深度学习）场景中应用呢？以下是一些主要的应用领域：

  • 音频分类问题：例如语音识别、说话人验证、说话人分离。
  • 音频增强：例如音频降噪、音频上采样。
  • 音乐信息检索：这是一个结合数字信号处理与机器学习的领域，用于解决如乐器识别、音乐情绪与流派分类等问题。

课程内容大纲 📚

本系列课程将涵盖大量内容，并且会根据大家的反馈进行调整。目前确定的核心主题包括：

  • 声波基础：理解音频的物理本质。
  • 数模/模数转换：了解数字音频是如何产生和记录的。
  • 音频特征：深入探讨时域和频域特征，例如RMS、频谱质心、MFCCs。
  • 音频变换：学习关键的数学变换，包括傅里叶变换、短时傅里叶变换（其结果是频谱图），并对比其他变换如常数Q变换、梅尔频谱图和色谱图。
  • 听觉感知：探讨如何利用人类听觉感知的知识来预处理音频数据，使其更贴合待解决的问题。

课程形式与资源 💻

如果你熟悉The Sound of AI频道，就会知道我喜欢同时涵盖理论和实践。本系列也不例外：

  • 理论环节：深入探讨所涉及主题背后的理论思想。
  • 编程环节：使用Python实现讨论过的理论内容。

课程将主要使用Python编程语言，并依赖一个强大的开源音频处理库——Librosa，来方便地提取各种音频特征。

所有课程资料，包括代码示例和幻灯片，都可以在下方描述区链接的GitHub页面找到，方便大家复习。

学习目标 🎓

从实践角度，完成本系列课程后，你将能够：

  • 深入理解音频数据，掌握其本质及预处理方法。
  • 熟悉时域和频域音频特征，并能从音频中提取这些特征。
  • 为不同的音频ML应用选择合适的音频特征。
  • 预处理音频数据，使其适用于深度学习模型。
  • 理解核心音频变换背后的数学原理，从而能够根据具体问题调整特征提取参数。
  • 高效使用Librosa库，为你的音频ML项目提取所需特征。

最重要的是，本课程成功的衡量标准是：当你再次看到类似下图的频谱图时，不再感到困惑，而是能够理解其含义并解读其中蕴含的音频信息。

目标受众 👥

本系列课程适合以下人群：

  • 涉足音频领域的机器学习/深度学习工程师。
  • 寻求音频数据预处理方法的计算机科学学生。
  • 对音频和音乐感兴趣的软件工程师。
  • 希望深入了解音频与计算的音乐技术专家或技术导向的音乐人。

预备知识与社区 🚀

本课程并非面向Python初学者，你需要具备中级的Python技能以跟上编程环节。

最后，我邀请你加入The Sound of AI的Slack社区。在那里，你可以找到一个由对AI音乐、AI音频、音频与音乐处理感兴趣的同好组成的成长型社区。你可以提问、深化理解，并与许多酷炫且博学的人交流网络。社区链接请在下方描述区查看。

本节课中，我们一起了解了《面向机器学习的音频信号处理》系列课程的开设背景、核心内容、学习目标以及适合的学习者。我期待与你一同开启这段学习旅程。

002：声音与波形 🎵

在本节课中，我们将学习声音的基本概念，包括声音的产生原理、波形的表示方法，以及频率、振幅、音高等核心参数的定义和关系。这些基础知识是理解后续音频处理技术的关键。

声音是什么？

声音是由振动物体产生的。物体振动导致空气分子振荡并相互碰撞，从而改变局部区域的空气压力状态，并在此过程中形成波。因此，我们可以将声音视为一种通过空气分子将能量从一点传递到另一点的波。

具体来说，声音是一种机械波。机械波是一种在空间中传播的振荡，能量随之传递。机械波的特点是它需要一个介质来传播。对于声音而言，这个介质通常是空气。当声波传播时，介质（空气）会发生形变：空气分子碰撞时形成高压区（压缩），分子相互远离时形成低压区（稀疏）。

我们可以用压力图来可视化声波。下图展示了一个简单的正弦波，中心线代表平均大气压。随时间变化，压缩点（空气分子密集）和稀疏点（空气分子稀疏）交替出现。

波形：声音的可视化

我们可以使用波形来表示复杂的声音。波形本质上是一个压力图，它绘制了空气压力相对于零基准线的偏差随时间的变化。

通过波形，我们可以直观地表示一段音乐、噪音或任何声音。它不仅是可视化声音的好方法，也是后续提取音频信号重要特征时进行各种变换的基础。

波形之所以重要，是因为它提供了丰富的信息，不仅包括频率信息，还包括强度、音色、持续时间等时间信息。例如，在一段复杂音乐的波形中，我们可以识别音符的起始点和持续时间。

声音的分类

根据其特性，声音可以分为以下几类：

  • 周期性声音：在这种声音中，压缩和稀疏有规律地重复。最简单的周期性声音是单一正弦波。更复杂的声音（如管弦乐）是由多个正弦波叠加组合而成的复合音。
  • 非周期性声音：顾名思义，这类声音没有周期性。它又可以分为两种：
    • 连续非周期性声音：例如噪音，其空气压力变化不遵循任何模式，是随机的。
    • 瞬态非周期性声音：例如爆裂声或咔嗒声，它们是能量的突然爆发，导致空气压力瞬间改变，类似于脉冲。

正弦波与核心参数

让我们从最简单的正弦波开始分析。下图展示了一个正弦波的波形及其数学表达式。

正弦波的方程由几个参数决定：
y(t) = A * sin(2πft + φ)
其中：

  • A 代表振幅。
  • f 代表频率。
  • t 代表时间。
  • φ 代表相位。

接下来，我们逐一分析这些参数如何影响波形。

频率与周期

周期是指波形中两个连续波峰（或波谷）之间的时间间隔。频率是周期的倒数，表示每秒完成的周期数，单位为赫兹（Hz）。
f = 1 / T
其中 T 是周期，f 是频率。

振幅

振幅衡量的是空气压力扰动的最大幅度。振幅越大，表示能量传递越强，我们感知到的声音就越响。

相位

相位决定了波形在时间零点（t=0）时的位置。通过改变相位，我们可以将波形向左或向右平移。

通过振幅、频率和相位，我们可以完全确定一个正弦波。正如我们将在后续课程中看到的，复杂声音可以被视为许多正弦波的组合或叠加，因此理解这些参数至关重要。

频率、振幅与听觉感知

频率和振幅与我们的听觉感知直接相关。

  • 频率与音高：频率越高，我们感知到的音调就越高。下图展示了不同频率（周期）的波形。
  • 振幅与响度：振幅越大，我们感知到的声音就越响。下图展示了不同振幅的波形。

听觉范围与音高

人类可听到的频率范围（听觉范围）大约在20 Hz到20,000 Hz之间。不同动物的听觉范围差异很大，例如猫、狗可以听到更高频率的超声波。

频率是一个客观的物理量，而我们主观感知到的频率称为音高。我们对音高的感知不是线性的，而是对数的。一个有趣的现象是，如果两个频率相差一倍（即2的幂次关系），我们感知到的音高差异是相同的，这引出了八度的概念。

MIDI音符、音高与频率的映射

在音乐中，我们常用MIDI音符来表示音高。每个琴键对应一个MIDI编号（例如中央C是60），也可以映射为音符名（如C4）。

一个八度包含12个半音。如果我们将频率提高一倍，音高就上升一个八度。例如，标准音A4（MIDI音符69）的频率是440 Hz，那么高一个八度的A5频率就是880 Hz。

频率 f 和音高 p（以MIDI编号表示）之间的转换公式如下：
f(p) = 440 * 2^((p - 69)/12)
例如，计算MIDI音符60（C4）的频率：
f(60) = 440 * 2^((60 - 69)/12) ≈ 261.6 Hz

相邻两个半音（即MIDI编号相差1）之间的频率比是恒定的：
f(p+1) / f(p) = 2^(1/12) ≈ 1.059

音分

我们的听觉分辨率很高，能够分辨小于一个半音的音高差异。为了精确度量这种微小差异，引入了音分的概念。一个八度被分为1200音分，因此每个半音对应100音分。人耳可察觉的最小音高差异阈值大约在10到25音分之间，音乐家通常具有更好的分辨能力。

总结

本节课我们一起学习了声音与波形的基础知识。我们了解了声音作为机械波的产生原理，以及如何使用波形来可视化声音。我们深入探讨了描述声音的核心参数：频率（及相关的周期）、振幅和相位，并明确了它们与听觉感知（音高、响度）之间的关系。此外，我们还介绍了人类听觉范围、音乐中音高与频率的对数映射关系（通过MIDI音符和八度概念），以及用于度量微小音高差异的单位——音分。

下一节课，我们将继续深入探讨声音的其他描述参数，包括与振幅相关的强度、功率和响度，以及声音中非常有趣且主观的一个属性——音色。

003：强度、响度与音色

在本节课中，我们将继续探索声音的特征，重点学习强度、响度与音色这三个核心概念。我们将了解声音的客观物理量（如强度）与主观感知量（如响度）之间的区别与联系，并初步探讨音色这一复杂而迷人的属性。

声音强度与功率

上一节我们介绍了声音的基本特征和频率。本节中，我们来看看声音的强度与功率。

声音功率是一个物理量，表示声源在单位时间内向所有方向辐射的总能量。其公式为：
[
P = \frac{E}{t}
]
其中，( P ) 是功率，单位是瓦特 (W)；( E ) 是能量；( t ) 是时间。

与声音功率紧密相关的是声音强度。声音强度定义为通过单位面积的声功率。其公式为：
[
I = \frac{P}{A}
]
其中，( I ) 是强度，单位是瓦特每平方米 (W/m²)；( P ) 是功率；( A ) 是面积。

一个有趣的事实是，一个大型交响乐团或一次强烈的雷鸣，其声音功率大约仅为1瓦特，这仅相当于一个旧式白炽灯泡功率的百分之一。

听觉阈值与分贝

人类能够感知的声音强度范围极其宽广。以下是两个关键的听觉阈值：

  • 听阈：这是人耳能感知到的最小声音强度，约为 (10^{-12}) W/m²。
  • 痛阈：这是开始引起听觉疼痛的强度阈值，约为 (10) W/m²。

由于人耳可感知的强度范围跨越了13个数量级，我们通常使用对数标度——分贝 (dB) 来描述声音的强度级。分贝的计算公式为：
[
L = 10 \log_{10}\left(\frac{I}{I_0}\right)
]
其中，( L ) 是强度级（分贝），( I ) 是当前声音强度，( I_0 ) 是参考强度（通常取听阈强度 (10^{-12}) W/m²）。当 ( I = I_0 ) 时，强度级为0 dB。

分贝每增加3 dB，声音强度大约翻倍。以下是一些常见声音的强度级参考：

  • 耳语：约 20 dB
  • 正常交谈：约 60 dB
  • 痛阈：约 130 dB
  • 喷气发动机起飞：约 140 dB

响度：主观的强度感知

虽然强度和功率是客观的物理量，但响度是我们对声音大小的主观感受。它与强度相关，但并非线性对应。

响度感知受以下因素影响：

  • 持续时间：在相同强度下，持续时间较长的声音听起来比短促的声音更响。
  • 频率：人耳对不同频率声音的敏感度不同。在中等频率范围（约几百赫兹到几千赫兹）内，我们对响度的感知最敏感。
  • 年龄：不同年龄的人对相同强度声音的响度感知也不同。

研究人员通过心理声学实验，引入了“方 (phon)”作为响度的度量单位。等响曲线图揭示了在相同响度水平下，所需声压级（分贝）随频率变化的规律。例如，要让人耳在20 Hz的低频处感知到与1 kHz处相同的响度，需要更高的声压级。

音色：声音的“色彩”

音色是区分具有相同音高、响度和时长的不同声音的属性。例如，同样演奏中央C，小提琴和钢琴的声音截然不同，这种差异就源于音色。

音色是一个多维度的复杂特征，尚无单一、全面的定义，但通常用以下词语描述：明亮、暗淡、柔和、刺耳、温暖等。影响音色的主要因素包括以下三个方面。

1. 声音包络

声音包络（或称振幅包络）描述了一个声音振幅随时间变化的形状。在音乐中，常用ADSR模型来描述：

  • 起音 (Attack)：声音振幅从静音上升到峰值的阶段。
  • 衰减 (Decay)：峰值过后，振幅下降到某一稳定水平的阶段。
  • 延音 (Sustain)：振幅保持相对稳定的阶段。
  • 释音 (Release)：声音振幅从稳定水平衰减到零的阶段。

不同乐器的包络形状差异显著。例如，钢琴的起音尖锐、衰减快；而小提琴的起音则相对平缓。

2. 谐波成分

复杂的声音可以看作是许多正弦波（称为“分音”或“谐波分音”）的叠加。

  • 基频：最低的分音，决定了我们感知到的音高。
  • 谐波：频率为基频整数倍的分音。例如，如果基频是440 Hz，那么二次谐波就是880 Hz，三次谐波是1320 Hz，依此类推。

音色在很大程度上取决于这些谐波分音上的能量分布。频谱图是一种可视化工具，可以展示声音中不同频率成分的能量随时间的变化。通过对比小提琴和钢琴演奏同一音符的频谱图，可以清晰地看到它们谐波能量分布的差异。

3. 振幅与频率调制

调制是指信号的某些参数（如振幅或频率）按照另一个信号（调制信号）的规律发生变化。

  • 频率调制 (FM) / 颤音 (Vibrato)：声音频率发生周期性的微小波动。在小提琴等乐器中，颤音常用于情感表达。
  • 振幅调制 (AM) / 震音 (Tremolo)：声音振幅发生周期性的波动。这种效果在音乐中也常能听到。

这两种调制都会显著改变声音的音色特征。

总结

本节课中，我们一起学习了声音的三个重要特征：

  1. 强度与功率：客观的物理量，描述了声音的能量属性，常用分贝度量其级别。
  2. 响度：人对声音强度的主观感受，受频率、持续时间等因素影响，使用“方”作为单位。
  3. 音色：区分声音“色彩”的多维度属性，主要受声音包络、谐波成分以及振幅/频率调制的影响。

理解这些基础概念，是后续深入学习音频数字信号处理、频谱分析以及机器学习应用的关键。下一节，我们将正式引入音频信号的概念，并探讨模拟信号与数字信号之间转换的核心过程——模数转换与数模转换。

004：理解音频信号

在本节课中，我们将学习如何将声音转换为数字化的音频信号，以便我们能够处理它、提取特征，并用于机器学习任务。首先，我们需要理解什么是音频信号。

音频信号是声音的一种可能表示形式，它包含了重建该声音所需的全部信息。然而，这里存在一个核心问题：声音本质上是模拟的机械波，而我们希望使用计算机等数字技术来处理它。那么，如何将模拟信号转换为数字信号呢？这正是本节课的主题。

在深入探讨之前，我们先简要概述一下模拟信号和数字信号。

模拟信号与数字信号

模拟信号的特点是，在表示时间的X轴和表示振幅（例如声压）的Y轴上，其值都是连续的实数。

如上图所示，模拟信号的曲线是连续的。如果我们想以数字格式存储模拟信号，就会遇到问题：它在时间和振幅轴上都具有无限的分辨率。这意味着无论我们如何提高观察精度，总能找到更细微的时间点或振幅值。显然，这需要无限的内存来存储，而我们无法负担。因此，我们需要转向数字信号。

数字信号本质上是一系列离散值的序列。可以将其理解为在连续信号的不同时间点拍摄快照。这些数据点只能取有限数量的值，而不是所有可能的实数。

模数转换

从模拟信号转换到数字信号的过程称为模数转换，其英文缩写是 ADC。这个过程包含两个子步骤：采样 和 量化。ADC的结果就是数字音频信号，我们通常也称之为脉冲编码调制信号。

上一节我们介绍了ADC的基本概念，本节中我们来看看它的两个核心步骤。

1. 采样

采样的过程，就是在特定时间点获取模拟波形上的数据点。下图中的黑点就是采样点。

我们通常以固定的采样周期（用大写 T 表示）在等间隔的时间点上进行采样。每个周期获取一个数据点。

要定位第 n 个采样点在时间轴上的位置，可以使用以下公式：

时间(T_n) = n × T

采样的一个关键特征是采样率，用 SR 表示。它是采样周期 T 的倒数：

SR = 1 / T

采样率的单位是赫兹，它表示每秒钟数字信号中包含的样本数量。

采样率的高低会影响采样质量。采样率越低，采样周期越长，采样点之间的间隔越大，这会导致采样误差增大——即连续曲线下的面积与采样点构建的阶梯状面积之间的差异。反之，采样率越高，时间分辨率越高，采样误差就越小。

那么，我们应该使用多高的采样率呢？这引出了一个重要概念：奈奎斯特频率。

奈奎斯特频率与混叠

根据奈奎斯特定理，奈奎斯特频率 F_N 等于采样率的一半：

F_N = SR / 2

奈奎斯特频率定义了一个数字信号中能够被无失真重建的最高频率上限。如果信号中包含高于奈奎斯特频率的成分，就会产生名为混叠的失真现象。

混叠会使高于奈奎斯特频率的信号成分“折叠”到较低的频率上，从而在重建信号中引入原本不存在的低频成分。下图展示了这一现象：红色高频原始信号经过采样（黑点）后，重建出的蓝色信号频率远低于原始信号。

为了直观感受混叠的效果，我们可以听一段音乐。在44.1 kHz的原始采样率下，音乐听起来是正常的。但如果将其重采样到极低的1 kHz，奈奎斯特频率变为500 Hz，所有高于500 Hz的频率成分都会发生混叠，导致声音严重失真、变得低沉模糊。

因此，在实际的ADC过程中，通常会加入一个抗混叠滤波器（低通滤波器），在采样前滤除信号中高于奈奎斯特频率的成分，以防止混叠发生。

理解了采样之后，我们来看看ADC的第二个步骤：量化。

2. 量化

量化处理的是信号的振幅轴。其核心思想是：将振幅轴划分为有限个离散的层级，然后将每个采样点的振幅值近似到最接近的可用层级。

与采样会产生采样误差类似，量化也会产生量化误差，即原始振幅值与量化后近似值之间的差异。直观上，量化分辨率越高（即振幅轴上可用的层级越多），量化误差就越小。

在数字系统中，这些振幅层级用二进制数值表示。量化的分辨率通常用比特深度来衡量。例如，CD的比特深度是16位。

比特深度为 Q 时，可表示的振幅层级数量为 2^Q。对于16位深度，就是 2^16 = 65,536 个不同的振幅值。

量化过程关联着动态范围的概念，它指的是系统能够记录的最大信号与最小信号之间的差异，可以理解为数字声音的响度范围。比特深度越高，动态范围通常也越大。

动态范围又与信噪比相关。信噪比衡量的是最大信号强度与量化误差（噪声）之间的关系。对于量化过程，信噪比 SQNR 可以近似估算为：

SQNR ≈ 6.02 × Q dB

其中 Q 是比特深度。对于16位深度，信噪比约为 96 dB，这也大致对应了其动态范围。

音频的录制与回放流程

现在，我们可以将整个过程串联起来，理解声音从录制到回放的完整数字链路：

  1. 录制：

    • 声音（机械波）冲击麦克风的振膜，使其振动。
    • 振动产生模拟电信号。
    • 模拟电信号进入声卡（充当ADC设备）。
    • 声卡进行抗混叠滤波、采样和量化，输出数字信号。
    • 数字信号被存储到计算机中。

  2. 回放：

    • 数字信号从计算机送入声卡。
    • 声卡进行数模转换，将数字信号还原为模拟电信号。
    • 模拟电信号驱动扬声器的振膜振动。
    • 振膜振动产生机械波（声音），传入我们的耳朵。

总结

本节课中，我们一起学习了音频数字化的核心知识：

  • 我们了解了模拟信号（连续值）与数字信号（离散值）的根本区别。
  • 我们深入探讨了模数转换的两个关键步骤：
    • 采样：在时间轴上以固定周期（采样率）获取离散点。采样率需满足奈奎斯特定理，以避免混叠失真。
    • 量化：在振幅轴上将连续值近似到有限的离散层级（由比特深度决定）。量化会引入误差，并决定了信号的动态范围和信噪比。
  • 我们梳理了声音从物理振动，经过ADC变为数字信号，再通过DAC还原为声音的完整流程。

掌握了音频信号的基础后，在接下来的课程中，我们将开始深入探讨如何从这些数字音频信号中提取各种特征，例如时域特征和频域特征。这些特征正是我们训练机器学习或深度学习模型所使用的核心数据。

欢迎加入Sound of AI Slack社区，与更多对AI、音频、音乐和数字信号处理感兴趣的同好交流。链接请在视频描述中查找。如果你对本视频内容有任何疑问，请在下方评论区提出。

005：机器学习中的音频特征类型 🎵

在本节课中，我们将学习什么是音频特征，以及如何从不同角度对它们进行分类。理解这些分类策略，有助于我们为机器学习任务选择合适的音频描述符。

音频特征是声音的描述符。不同的音频特征能为我们提供声音不同方面的信息。我们可以利用这些特征来训练智能音频系统。具体方法是，识别出我们感兴趣的若干音频特征，然后将它们输入到机器学习系统中，以期让系统学习到模式并解决我们的任务或问题。

现在，让我们回顾几种可用于对音频特征进行分类的策略。这里列出了五种策略。这些策略的设计初衷是尽可能通用，能够处理任何类型的声音。当然，其中一些策略（如“音乐层面”和“抽象层次”）更侧重于音乐信号。

以下是五种分类策略：

  • 抽象层次
  • 时间范围
  • 音乐层面
  • 信号域
  • 机器学习方法

在接下来的内容中，我们将深入探讨每一种分类策略。

抽象层次 🎼

这种分类策略主要适用于音乐信号，而非一般的声音或音频。我们可以将其分为三个层次：低层音频特征、中层音频特征和高层音频特征。

以下是三个层次的具体说明：

  • 低层特征：这些特征对机器有意义，但对我们人类来说不太容易理解。它们通常是直接从音频中提取的统计特征。未经音频处理训练的人很难理解这些特征的含义。
  • 中层特征：这些特征开始从感知角度变得有意义。它们涉及音高和节拍等属性。例如，音符起始点、波动模式和梅尔频率倒谱系数都属于此类。
  • 高层特征：这些是非常抽象的特征，倾向于映射到我们在听音乐时能够理解和感知的音乐结构。例如，乐器、调性、和弦、旋律、节奏和速度。基本思想是，层次越高，音频特征就越抽象。

时间范围 ⏱️

这种策略适用于任何类型的声音，无论是音乐还是非音乐。我们可以将音频特征大致分为三类。

以下是基于时间范围的分类：

  • 瞬时级特征：这些特征提供音频信号的瞬时信息。它们通常考虑非常短的音频片段，大约在20到100毫秒之间。需要记住，人类能够感知的最小时间分辨率大约为10毫秒，低于这个阈值在感知层面上就没有意义了。
  • 片段级特征：这些是可以在音频片段上计算的特征。这里我们谈论的是秒级的时间，从2秒到5秒、10秒、15秒、20秒甚至更长。以音乐为例，这些特征可以提供关于一个小节或一个乐句的信息。它们倾向于从更长的时段中聚合瞬时信息。
  • 全局特征：这些特征提供关于整个声音的信息，基本上是聚合特征。我们可以聚合来自较低时间分辨率特征（如瞬时或片段级特征）的结果，然后使用某种平均值或更复杂的聚合方法，最终得到一个描述整个声音或整个信号的数字或特征向量。

音乐层面 🎹

这种分类策略显然只专注于音乐。它基于音频特征所揭示的音乐层面进行分类。

例如，音符起始点与节拍和旋律有关。而调性则与和声和音高有关。我们可以根据音频特征所反映的音乐层面来对它们进行分类。由于显而易见的原因，这种分类策略主要用于音乐。

信号域 📊

这是对音频特征进行分类的最重要策略之一，它基于我们所在的信号域。通过几个例子，这个概念会变得非常清晰。

以下是基于信号域的分类：

  • 时域特征：这些特征直接从波形（原始音频）中提取。在波形中，X轴是时间，Y轴是振幅。时域音频特征从这个表示中提取信息。例如：振幅包络、均方根能量、过零率。
  • 频域特征：声音的特征很大程度上由频率决定，而时域表示不包含频率信息。频域特征专注于声音的频率成分。我们通过对时域信号应用傅里叶变换，将其转换到频域，从而得到频谱。在频谱中，X轴是频率，Y轴是幅度。例如：频带能量比、频谱质心、频谱通量、频谱扩展。
  • 时频域特征：时域和频域表示无法同时提供时间和频率信息。时频域特征则能同时提供这两方面的信息，使用时频表示。例如：频谱图、梅尔频谱图、常数Q变换。其中，频谱图是最著名的时频表示。通过对时域信号应用短时傅里叶变换获得。在频谱图中，X轴是时间，Y轴是频率，颜色亮度表示特定时间点上特定频带的贡献大小。

机器学习方法 🤖

最后一种分类音频特征的方法是基于所使用的机器学习方法。这里可以在传统机器学习算法和深度学习架构之间做出重要区分。

以下是基于机器学习方法的分类：

  • 传统机器学习：对于传统机器学习算法（如支持向量机、逻辑回归），我们通常会手动挑选我们认为对解决问题最有效的音频特征（来自时域和频域）。例如，对于音频分类任务（区分汽车引擎、飞机或枪声），我们可能选择振幅包络、过零率和频谱通量这三个特征，将它们从音频文件中提取出来，然后输入给算法进行训练和推理。
  • 深度学习：在深度学习中，我们倾向于使用非结构化的数据表示。类似于图像处理中直接输入像素，在音频处理中，我们可以直接输入原始音频（时域表示），或者输入非结构化的音频表示，如频谱图、梅尔频谱图等。深度学习的承诺是算法能够自动从这些原始表示中提取特征，而无需人工手动挑选。

总结与回顾 📝

在本节课中，我们一起学习了音频特征及其多种分类策略。我们回顾了五种主要的分类方式：抽象层次、时间范围、音乐层面、信号域和机器学习方法。

其中，基于信号域的分类（时域、频域、时频域）是最核心和通用的策略，这将是我们未来课程中重点探讨的内容。

理解这些分类后，我们现在对将要使用的“食材”——即各种音频特征——有了清晰的图景。下一步是学习如何直接从音频中提取这些特征。在下一节课中，我们将探讨针对时域和频域特征的特征提取流程。

希望本节课对你有所帮助。我们下节课再见！

006：音频特征提取流程

在本节课中，我们将学习如何构建提取时域和频域音频特征的完整流程。我们将从模拟声音开始，逐步了解数字化、分帧、加窗等关键步骤，并解释为何需要这些处理。

概述

在上一节视频中，我们介绍了几种不同类型的音频特征：存在于时域的时域特征、存在于频域的频域特征，以及同时提供时间和频率信息的时频域特征。

本节我们将重点探讨用于提取时域和频域特征的提取流程。让我们先从时域特征开始。

时域特征提取流程

时域特征的提取流程始于原始声音。

首先，我们需要处理一个模拟声音，例如小提琴声或噪音。为了能在计算机中编辑和处理这个声音，我们必须将其转换为数字信号。这个过程称为模数转换，包括采样和量化两个步骤。

如果你对这个过程不熟悉，建议回顾本系列中详细介绍此过程的视频。

获得声音的数字化版本后，下一步是分帧。分帧意味着将一系列样本捆绑在一起。例如，第一帧可能包含第1到第128个样本，第二帧包含第64到第192个样本，第三帧包含第128到第256个样本，依此类推。

这里有一个有趣的现象：这些帧是重叠的。我们暂时不解释原因，在本视频结束时你会明白。

为何需要分帧？

在提取声学特征之前，我们需要理解为何使用分帧。可以将帧视为可感知的音频块。

问题在于，如果我们以CD标准的44.1 kHz采样率来看，一个单一样本的持续时间仅为约0.0227毫秒。这个时间非常短，远低于人耳听觉的时间分辨率阈值（约10毫秒）。这意味着持续时间低于10毫秒的事件，我们几乎无法将其感知为独立的声学事件。

通过分帧，我们获得一段有足够持续时间的音频信号，使其从声学角度变得有意义，并能被人类感知。我们想要提取的特征与人类的听觉体验相关，因此分帧是必要的。

帧的大小与持续时间

关于帧的另一个要点是，帧大小（即每帧包含的样本数）通常是2的幂。这听起来可能很奇怪，原因在于当我们进入频域进行处理时，通常会应用傅里叶变换。而应用快速傅里叶变换（FFT）时，样本数为2的幂可以极大地加速计算过程。

典型的帧大小在256到8192之间。帧的持续时间可以通过以下公式计算：

帧持续时间 = (1 / 采样率) × 帧大小

这个公式中，(1 / 采样率) 计算的是单个样本的持续时间，乘以帧大小 K 就得到了整帧的持续时间。

让我们代入一些常用数值。假设采样率为44.1 kHz，帧大小为512，计算可得帧持续时间约为11.6毫秒。这个值略高于人耳约10毫秒的时间分辨率。

特征计算与聚合

回到时域特征提取流程。在完成分帧后，下一步是在每一帧上计算时域特征。

但之后，我们需要聚合这些结果，为整个声音得到一个单一的特征描述。这可以通过使用统计方法来实现，例如计算均值、中位数、总和，或使用更复杂的方法如高斯混合模型。

通过这个过程，我们最终得到一个值、一个特征向量，甚至一个特征矩阵。这基本上是整个被分析音频信号的一个快照。

以上就是时域特征的提取流程。

频域特征提取流程

现在让我们转向频域特征。你会发现，许多步骤与上述时域流程基本相同。

我们同样从模拟声音开始，进行模数转换（采样和量化），得到数字化音频信号，然后对其进行分帧。

现在我们有了一系列帧。下一步是使用一个“魔法工具”——傅里叶变换——将信号从时域转换到频域。目前，我们可以将其视为一个黑箱，它神奇地将时间表示转换为频率表示。我们将在后续视频中更详细地介绍傅里叶变换。

简单回顾一下时域和频域：

  • 时域：声音被可视化为随时间变化的函数，我们看到声音在时间轴上的所有事件。
  • 频域：我们观察声音的所有频率成分，看它们对整体声音的贡献有多大。X轴是频率（Hz），Y轴是幅度，表示各频段对整体声音的贡献程度。

频谱泄漏问题

你可能会认为，提取流程中下一步自然就是应用傅里叶变换。但不幸的是，这里存在一个主要问题，称为频谱泄漏。

频谱泄漏发生在我们处理的信号周期数不是整数时。这在实际中经常发生，因为一段音频的时间长度很少恰好包含整数个周期。

通常，信号的端点会存在不连续性，因为它们不是整数个周期。

这些不连续性在通过傅里叶变换转换到频谱（频域）时，会表现为高频成分。但这些高频成分在原始信号中并不真实存在，它们只是由于端点不连续性而产生的伪影。换句话说，这些不连续频率“泄漏”到了其他更高频率中，因此得名“频谱泄漏”。

在频谱中，你可能会看到一些本不该有实质性贡献的高频成分（例如下图中红框部分），这就是频谱泄漏的例子。

加窗：解决频谱泄漏

幸运的是，我们有一种方法可以最小化频谱泄漏，那就是加窗。

加窗的思想是：在将帧送入傅里叶变换之前，对每一帧应用一个窗函数。这样做基本上消除了帧端点的样本信息。

这会产生一个周期性信号，从而最小化频谱泄漏。

一个著名且最常用的窗函数是汉宁窗。其函数形式如下（其中 k 代表样本索引）：

w(k) = 0.5 * (1 - cos(2πk / (K-1)))， 其中 0 ≤ k ≤ K-1

下图展示了应用于50个样本的汉宁窗，它形成了一个钟形曲线，端点趋向于0。

将汉宁窗应用于原始信号 s 的数学操作是逐点相乘：s_windowed = s * w。得到的新信号端点被平滑，不连续性消失，因为信号在窗函数的作用下自然趋于零。

重叠帧：解决信息丢失

但加窗带来了一个新问题：我们丢失了信号。因为在帧的起点和终点，我们通过加窗移除了信息。我们当然不希望在进行傅里叶变换时丢失信号。

解决方案就是使用重叠帧。这解开了视频开头提到的“帧为何重叠”的谜团。

在非重叠帧的情况下，我们选择一个帧大小并将其应用于信号，帧与帧之间没有重叠。但如果加窗后这样做，就会丢失端点信息。

通过使用重叠帧，我们可以弥补在端点丢失的信息。我们让连续的帧之间有一部分样本重叠。

关键概念：帧大小与跃迁长度

这里有两个重要概念：

  1. 帧大小：每帧包含的样本数量。
  2. 跃迁长度（有时称为跃迁大小）：每次我们向右移动以获取新一帧时，跳过的样本数。

在频域特征提取流程中，经过加窗后，我们现在可以应用傅里叶变换了。我们希望得到的是一个频谱泄漏最小化的频谱。

之后，我们回到与时域特征提取相似的步骤：在每一帧上计算频域特征，然后使用统计方法聚合整个音频信号的所有结果，最终得到一个频域特征值、向量或矩阵。

总结

本节课我们一起学习了提取时域和频域音频特征的完整流程。我们了解了从模拟信号到数字信号的转换、分帧的必要性、加窗如何解决频谱泄漏问题，以及为何需要使用重叠帧来避免信息丢失。

你现在应该对处理时域和频域特征所需的流程有了清晰的认识。接下来，我们将开始深入探讨具体的时域特征，了解它们是什么以及如何在机器学习的不同应用中使用它们。敬请期待下一节内容。

007：理解时域音频特征

在本节课中，我们将要学习三种重要的时域音频特征：振幅包络、均方根能量和过零率。我们将逐一探讨它们的定义、计算方法、应用场景以及各自的优缺点。

概述

上一节我们介绍了时域和频域特征的提取流程。本节中，我们来看看三种核心的时域音频特征：振幅包络、均方根能量和过零率。这些特征是低级的声学特征，并且是瞬时性的，即我们在音频信号的每一帧上计算它们。理解这些特征是进行更高级音频分析的基础。

振幅包络

振幅包络是我们在一个音频帧中所有样本里取出的最大振幅值。它为我们提供了信号响度的粗略概念，因为振幅与强度和响度相关。

计算方法

振幅包络在帧 t 的计算公式如下：

AE(t) = max(|s(k)|)  对于 k = t*K, ..., (t+1)*K - 1

其中：

  • s(k) 是样本 k 的振幅值。
  • K 是帧大小，即一帧中包含的样本数量。
  • t 是当前帧的索引。

这个公式的含义是：找出从第 t 帧的第一个样本（索引为 t*K）到最后一个样本（索引为 (t+1)*K - 1）之间，振幅绝对值的最大值。

可视化理解

想象我们将一个音频信号分割成连续的帧。对于每一帧，我们找出该帧波形中振幅绝对值最大的那个点。连接所有帧的这些最大振幅点，就形成了振幅包络线。它勾勒出了信号振幅变化的轮廓。

特点与应用

振幅包络的主要问题是对离群值非常敏感。因为它是从一帧中只取一个样本值（最大值），如果该帧内恰好有一个异常的振幅尖峰，那么这个尖峰就会成为代表整帧的包络值，可能无法真实反映整帧的平均能量水平。

以下是振幅包络的一些应用场景：

  • 起始点检测：在音符、单词或音素开始时，通常会出现能量的爆发（振幅尖峰），振幅包络可用于检测这些起始点。
  • 高级分类问题：例如音乐流派分类，不同风格的音乐可能具有不同的振幅包络形态。

均方根能量

均方根能量计算的是音频帧内所有样本的均方根值。它是衡量信号响度的另一个指标，并且比振幅包络更常用。

计算方法

均方根能量在帧 t 的计算公式如下：

RMS(t) = sqrt( (1/K) * Σ (s(k)^2) )  对于 k = t*K, ..., (t+1)*K - 1

计算步骤分解：

  1. 计算能量：对帧内每个样本的振幅 s(k) 进行平方，得到该样本的能量 s(k)^2。
  2. 求和：将帧内所有样本的能量值相加 Σ (s(k)^2)。
  3. 求平均：将总能量除以帧大小 K，得到平均能量 (1/K) * Σ (s(k)^2)。
  4. 开方：对平均能量开平方根，得到均方根能量 sqrt(平均能量)。

特点与应用

均方根能量对离群值的敏感度低于振幅包络。因为它考虑了帧内所有样本的贡献，并通过平均过程平滑了单个样本的极端值影响，因此能更稳定地表示整帧的能量水平。

以下是均方根能量的一些应用场景：

  • 音频分割：当音频中出现新的段落或事件（如说话人切换、音乐段落改变）时，RMS值通常会发生显著变化，可用于识别这些边界。
  • 音乐流派分类：与振幅包络类似，RMS能量特征也可用于区分不同音乐风格。

过零率

过零率是指信号在单位时间内穿过零轴（即从正振幅变为负振幅，或反之）的次数。它是一个直观的特征，广泛应用于语音识别和音乐信息检索。

计算方法

过零率在帧 t 的计算公式如下：

ZCR(t) = (1/(2K)) * Σ |sgn(s(k)) - sgn(s(k+1))|  对于 k = t*K, ..., (t+1)*K - 2

其中 sgn() 是符号函数：

  • 如果 s(k) > 0，则 sgn(s(k)) = 1
  • 如果 s(k) < 0，则 sgn(s(k)) = -1
  • 如果 s(k) = 0，则 sgn(s(k)) = 0

计算逻辑：

  1. 比较相邻样本：对于帧内每一对连续的样本 s(k) 和 s(k+1)，计算它们符号的差值 sgn(s(k)) - sgn(s(k+1))。
  2. 判断过零：
    • 如果两个样本符号相同（同为正或同为负），差值的绝对值为0，表示没有发生过零。
    • 如果两个样本符号相反（一正一负），差值的绝对值为2。
  3. 计数与归一化：将所有差值的绝对值求和，然后除以 2K（因为一次过零会产生值为2的贡献，且最多有 K-1 个样本对），得到该帧的过零率。

特点与应用

过零率提供了关于信号频率内容的粗略信息。高频信号通常比低频信号过零更频繁。

以下是过零率的一些应用场景：

  • 区分打击乐与音高乐音：打击乐声音的过零率通常变化随机且剧烈，而具有稳定音高的声音其过零率则相对稳定。
  • 简单的单音高估计：对于单音信号，过零率与基频存在一定的近似关系，过零率越高，通常音高也越高。这是一种基础但非绝对可靠的方法。
  • 语音识别中清/浊音判别：浊音（如元音）部分的过零率通常较低，而清音（如/s/、/f/等辅音）部分由于类似噪声的特性，过零率通常较高。

总结

本节课中我们一起学习了三种基础的时域音频特征：

  1. 振幅包络：提取每帧的最大振幅，反映响度轮廓，但对噪声敏感。
  2. 均方根能量：计算每帧样本的均方根值，是更稳健的响度度量。
  3. 过零率：统计每帧内信号穿过零轴的次数，与信号频率特性相关。

这些特征是构建更复杂音频分析系统的基石。在接下来的课程中，我们将动手实践，使用代码实现这些特征，并将它们应用于具体的音频分析任务中。

008：使用Python从零开始提取振幅包络特征 🎵

在本节课中，我们将学习如何使用Python和Librosa库从零开始提取音频信号的振幅包络特征。我们将加载音频文件、绘制波形图，并实现计算振幅包络的算法。

概述

我们将使用Librosa库加载音频文件，并实现一个函数来计算振幅包络。振幅包络是音频信号在每个时间帧内振幅的最大值，它能有效反映信号的整体能量变化。我们还将学习如何可视化波形和振幅包络。

导入必要的库

首先，我们需要导入处理音频和绘图所需的库。

import librosa
import librosa.display
import IPython.display as ipd
import matplotlib.pyplot as plt
import numpy as np

加载音频文件

我们将加载三个不同风格的音乐片段：古典乐、摇滚乐和爵士乐。

# 定义音频文件路径
bach_file = ‘audio/bach.wav’
redhot_file = ‘audio/red.wav’
duke_file = ‘audio/duke.wav’

# 使用IPython播放音频
ipd.Audio(bach_file)
ipd.Audio(redhot_file)
ipd.Audio(duke_file)

使用Librosa加载音频

Librosa的load函数可以将音频文件加载为单声道信号，并返回信号数据及其采样率。

# 加载音频文件，默认采样率为22050 Hz，转换为单声道
bach, sr = librosa.load(bach_file)
redhot, _ = librosa.load(redhot_file)
duke, _ = librosa.load(duke_file)

计算信号时长

我们可以通过采样率和样本总数来计算音频信号的时长。

# 计算单个样本的时长
sample_duration = 1 / sr
print(f“Duration of one sample: {sample_duration:.6f} seconds”)

# 计算整个信号的时长
duration = sample_duration * len(bach)
print(f“Duration of signal: {duration:.2f} seconds”)

可视化波形

为了直观比较不同音乐的波形，我们将使用Matplotlib绘制它们的波形图。

# 创建图形，设置大小
plt.figure(figsize=(15, 17))

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/391cbaece27f64c284238d800f2bc1c4_15.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/391cbaece27f64c284238d800f2bc1c4_16.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/391cbaece27f64c284238d800f2bc1c4_17.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/391cbaece27f64c284238d800f2bc1c4_18.png)

# 绘制巴赫音频的波形
plt.subplot(3, 1, 1)
librosa.display.waveplot(bach, alpha=0.5)
plt.title(‘Bach’)
plt.ylim(-1, 1)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/391cbaece27f64c284238d800f2bc1c4_19.png)

# 绘制红辣椒乐队音频的波形
plt.subplot(3, 1, 2)
librosa.display.waveplot(redhot, alpha=0.5)
plt.title(‘Red Hot Chili Peppers’)
plt.ylim(-1, 1)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/391cbaece27f64c284238d800f2bc1c4_21.png)

# 绘制艾灵顿公爵音频的波形
plt.subplot(3, 1, 3)
librosa.display.waveplot(duke, alpha=0.5)
plt.title(‘Duke Ellington’)
plt.ylim(-1, 1)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/391cbaece27f64c284238d800f2bc1c4_22.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/391cbaece27f64c284238d800f2bc1c4_23.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/391cbaece27f64c284238d800f2bc1c4_24.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/391cbaece27f64c284238d800f2bc1c4_26.png)

plt.show()

通过观察波形图，我们可以发现古典音乐的振幅变化范围较大，而流行音乐的波形则相对稳定。

计算振幅包络

振幅包络反映了信号在每个时间帧内的最大振幅。以下是两种实现方式。

基础实现方法

这种方法通过循环遍历每个帧来计算振幅包络，易于理解。

def amplitude_envelope(signal, frame_size, hop_length):
    “”“计算信号的振幅包络。”“”
    ae = []
    for i in range(0, len(signal), hop_length):
        current_frame = signal[i:i + frame_size]
        ae.append(max(current_frame))
    return np.array(ae)

# 设置帧大小和跳跃长度
FRAME_SIZE = 1024
HOP_LENGTH = 512

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/391cbaece27f64c284238d800f2bc1c4_28.png)

# 计算巴赫音频的振幅包络
ae_bach = amplitude_envelope(bach, FRAME_SIZE, HOP_LENGTH)
print(f“Amplitude envelope size: {ae_bach.size}”)

进阶实现方法

这种方法使用列表推导式，代码更简洁，但算法逻辑相同。

def fancy_amplitude_envelope(signal, frame_size, hop_length):
    “”“使用列表推导式计算振幅包络。”“”
    return np.array([max(signal[i:i + frame_size]) for i in range(0, len(signal), hop_length)])

# 计算并验证两种方法的结果是否一致
fancy_ae_bach = fancy_amplitude_envelope(bach, FRAME_SIZE, HOP_LENGTH)
print(np.array_equal(ae_bach, fancy_ae_bach))  # 应输出 True

可视化振幅包络

现在，我们可以在波形图上叠加绘制振幅包络，以更直观地观察信号的能量变化。

# 计算其他音频的振幅包络
ae_redhot = amplitude_envelope(redhot, FRAME_SIZE, HOP_LENGTH)
ae_duke = amplitude_envelope(duke, FRAME_SIZE, HOP_LENGTH)

# 将帧索引转换为时间
frames = range(ae_bach.size)
time = librosa.frames_to_time(frames, hop_length=HOP_LENGTH)

# 绘制波形和振幅包络
plt.figure(figsize=(15, 17))

# 巴赫音频
plt.subplot(3, 1, 1)
librosa.display.waveplot(bach, alpha=0.5)
plt.plot(time, ae_bach, color=‘r’)
plt.title(‘Bach with Amplitude Envelope’)
plt.ylim(-1, 1)

# 红辣椒乐队音频
plt.subplot(3, 1, 2)
librosa.display.waveplot(redhot, alpha=0.5)
plt.plot(time, ae_redhot, color=‘r’)
plt.title(‘Red Hot Chili Peppers with Amplitude Envelope’)
plt.ylim(-1, 1)

# 艾灵顿公爵音频
plt.subplot(3, 1, 3)
librosa.display.waveplot(duke, alpha=0.5)
plt.plot(time, ae_duke, color=‘r’)
plt.title(‘Duke Ellington with Amplitude Envelope’)
plt.ylim(-1, 1)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/391cbaece27f64c284238d800f2bc1c4_32.png)

plt.show()

从图中可以看出，振幅包络（红色曲线）紧密跟随波形的轮廓。在摇滚乐中，振幅包络的峰值通常对应鼓点，而古典音乐的振幅包络则表现出更大幅度的起伏。

总结

本节课中，我们一起学习了如何使用Python和Librosa库从零开始提取振幅包络特征。我们掌握了加载音频文件、计算信号时长、绘制波形图以及实现振幅包络算法的方法。通过比较不同音乐风格的振幅包络，我们了解到这一特征如何反映音频信号的动态变化。

009：如何从音频中提取均方根能量与过零率 🎵

在本节课中，我们将学习如何从音频信号中提取两种重要的时域特征：均方根能量 和 过零率。我们将使用 librosa 库来实现这些功能，并会从零开始编写代码来理解其背后的原理。

概述与准备工作

首先，我们需要导入必要的库并加载一些音频文件用于演示。我们将使用三首不同风格的音乐片段：一首德彪西的管弦乐、一首艾灵顿公爵的爵士乐和一首红辣椒乐队的摇滚乐。

import matplotlib.pyplot as plt
import numpy as np
import librosa
import librosa.display
import IPython.display as ipd

接下来，我们加载音频文件。为了节省时间，我们直接复用之前课程中的代码。

# 定义音频文件路径
debussy_path = ‘audio/debussy.wav’
redhot_path = ‘audio/redhot.wav’
duke_path = ‘audio/duke.wav’

# 加载音频信号
debussy, _ = librosa.load(debussy_path)
redhot, _ = librosa.load(redhot_path)
duke, _ = librosa.load(duke_path)

现在我们已经准备好了音频数据，可以开始提取特征了。

使用Librosa提取均方根能量

均方根能量是衡量音频信号在一段时间内平均能量的指标。它比振幅包络更平滑，因为它考虑了帧内所有样本的贡献。

以下是使用 librosa 提取均方根能量的步骤。首先，我们定义帧长和跳数。

FRAME_LENGTH = 1024
HOP_LENGTH = 512

然后，我们使用 librosa.feature.rms 函数来计算每首音乐的均方根能量。

rms_debussy = librosa.feature.rms(y=debussy, frame_length=FRAME_LENGTH, hop_length=HOP_LENGTH)
rms_redhot = librosa.feature.rms(y=redhot, frame_length=FRAME_LENGTH, hop_length=HOP_LENGTH)
rms_duke = librosa.feature.rms(y=duke, frame_length=FRAME_LENGTH, hop_length=HOP_LENGTH)

该函数返回一个NumPy数组。为了绘图方便，我们通常需要提取其第一行数据。

rms_debussy = rms_debussy[0]
rms_redhot = rms_redhot[0]
rms_duke = rms_duke[0]

现在，我们可以将均方根能量与原始波形一起可视化。

frames = range(len(rms_debussy))
t = librosa.frames_to_time(frames, hop_length=HOP_LENGTH)

plt.figure(figsize=(15, 17))

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_22.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_23.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_24.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_25.png)

# 绘制德彪西的波形和RMS
ax = plt.subplot(3, 1, 1)
librosa.display.waveshow(debussy, alpha=0.5)
plt.plot(t, rms_debussy, color=‘r’)
plt.title(‘Debussy’)

# 绘制红辣椒乐队的波形和RMS
plt.subplot(3, 1, 2)
librosa.display.waveshow(redhot, alpha=0.5)
plt.plot(t, rms_redhot, color=‘r’)
plt.title(‘Red Hot Chili Peppers’)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_27.png)

# 绘制艾灵顿公爵的波形和RMS
plt.subplot(3, 1, 3)
librosa.display.waveshow(duke, alpha=0.5)
plt.plot(t, rms_duke, color=‘r’)
plt.title(‘Duke Ellington’)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_29.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_31.png)

plt.show()

从图中可以看出，均方根能量曲线比之前学习的振幅包络曲线平滑得多，因为它平均了帧内所有样本的能量，而不是只取最大值。

从零实现均方根能量计算

为了深入理解均方根能量的计算原理，我们现在尝试自己编写一个函数来实现它。

均方根能量的计算公式如下：

RMS(t) = sqrt( (1/N) * Σ_{n=0}^{N-1} x(t, n)^2 )

其中，N 是帧长，x(t, n) 是第 t 帧中的第 n 个样本。

以下是该函数的实现代码：

def rms_manual(signal, frame_length, hop_length):
    rms = []
    for i in range(0, len(signal), hop_length):
        current_frame = signal[i:i+frame_length]
        frame_energy = np.sum(current_frame ** 2)
        frame_rms = np.sqrt(frame_energy / frame_length)
        rms.append(frame_rms)
    return np.array(rms)

现在，我们用这个自定义函数来计算均方根能量，并与 librosa 的结果进行比较。

rms_debussy_manual = rms_manual(debussy, FRAME_LENGTH, HOP_LENGTH)
rms_redhot_manual = rms_manual(redhot, FRAME_LENGTH, HOP_LENGTH)
rms_duke_manual = rms_manual(duke, FRAME_LENGTH, HOP_LENGTH)

我们将自定义函数的结果（黄色）与 librosa 的结果（红色）绘制在同一张图上进行对比。

plt.figure(figsize=(15, 17))

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_47.png)

ax = plt.subplot(3, 1, 1)
librosa.display.waveshow(debussy, alpha=0.5)
plt.plot(t, rms_debussy, color=‘r’)
plt.plot(t, rms_debussy_manual, color=‘y’)
plt.title(‘Debussy’)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_49.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_50.png)

plt.subplot(3, 1, 2)
librosa.display.waveshow(redhot, alpha=0.5)
plt.plot(t, rms_redhot, color=‘r’)
plt.plot(t, rms_redhot_manual, color=‘y’)
plt.title(‘Red Hot Chili Peppers’)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_52.png)

plt.subplot(3, 1, 3)
librosa.display.waveshow(duke, alpha=0.5)
plt.plot(t, rms_duke, color=‘r’)
plt.plot(t, rms_duke_manual, color=‘y’)
plt.title(‘Duke Ellington’)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_54.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_56.png)

plt.show()

可以看到，两条曲线几乎完全重合，这验证了我们自定义函数的正确性。

提取过零率

过零率是指音频信号在单位时间内穿过零点的次数。它是一个简单的时域特征，常用于语音/音乐分类和打击乐检测。

使用 librosa 提取过零率非常简单。

zcr_debussy = librosa.feature.zero_crossing_rate(debussy, frame_length=FRAME_LENGTH, hop_length=HOP_LENGTH)[0]
zcr_redhot = librosa.feature.zero_crossing_rate(redhot, frame_length=FRAME_LENGTH, hop_length=HOP_LENGTH)[0]
zcr_duke = librosa.feature.zero_crossing_rate(duke, frame_length=FRAME_LENGTH, hop_length=HOP_LENGTH)[0]

librosa 返回的过零率是归一化的值（除以帧长）。为了在同一图中比较三首音乐，我们进行可视化。

frames = range(len(zcr_debussy))
t = librosa.frames_to_time(frames, hop_length=HOP_LENGTH)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_59.png)

plt.figure(figsize=(10, 5))
plt.plot(t, zcr_debussy, color=‘y’, label=‘Debussy’)
plt.plot(t, zcr_redhot, color=‘r’, label=‘Red Hot Chili Peppers’)
plt.plot(t, zcr_duke, color=‘b’, label=‘Duke Ellington’)
plt.ylim(0, 1) # 归一化值范围在0到1之间
plt.legend()
plt.show()

从图中可以看出，红辣椒乐队的摇滚乐过零率最高，这是因为摇滚乐中打击乐成分较多。古典乐（德彪西）的过零率较低，而爵士乐（艾灵顿公爵）则介于两者之间。

如果我们想要得到实际的过零次数（而非归一化值），只需将结果乘以帧长即可。

zcr_debussy_abs = zcr_debussy * FRAME_LENGTH
zcr_redhot_abs = zcr_redhot * FRAME_LENGTH
zcr_duke_abs = zcr_duke * FRAME_LENGTH

语音与噪声的过零率对比

最后，我们通过对比语音和纯噪声的过零率，来进一步理解这个特征。我们加载一段15秒的语音和一段白噪声。

voice_path = ‘audio/voice.wav’
noise_path = ‘audio/noise.wav’

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/b3729772e16d0e187e40453400f2ad8b_80.png)

voice, _ = librosa.load(voice_path, duration=15)
noise, _ = librosa.load(noise_path, duration=15)

计算它们的过零率。

zcr_voice = librosa.feature.zero_crossing_rate(voice, frame_length=FRAME_LENGTH, hop_length=HOP_LENGTH)[0]
zcr_noise = librosa.feature.zero_crossing_rate(noise, frame_length=FRAME_LENGTH, hop_length=HOP_LENGTH)[0]

然后进行可视化比较。

frames_voice = range(len(zcr_voice))
t_voice = librosa.frames_to_time(frames_voice, hop_length=HOP_LENGTH)

plt.figure(figsize=(10, 5))
plt.plot(t_voice, zcr_voice, color=‘y’, label=‘Voice’)
plt.plot(t_voice, zcr_noise, color=‘r’, label=‘White Noise’)
plt.ylim(0, 1)
plt.legend()
plt.show()

正如预期，白噪声的过零率明显高于语音信号，因为白噪声的随机性导致其更频繁地穿过零点。

总结

本节课中，我们一起学习了两种重要的音频时域特征：

  1. 均方根能量：它提供了信号能量的平滑估计，计算公式为 RMS = sqrt( (1/N) * Σ x(n)^2 )。我们使用 librosa.feature.rms 进行提取，并成功从零实现了该计算。
  2. 过零率：它衡量信号符号变化的速率，是区分语音、音乐和噪声的简单有效特征。我们使用 librosa.feature.zero_crossing_rate 进行提取，并观察了它在不同音乐风格以及语音与噪声之间的差异。

至此，关于时域音频特征的学习就告一段落了。下一节课，我们将进入频域，开始学习信号处理中一个极其重要的工具——傅里叶变换。

010：直观理解

在本节课中，我们将要学习傅里叶变换的核心概念。傅里叶变换是信号处理的基石，它能够将复杂的时域信号分解为其组成的频率分量。我们将通过直观的类比和简单的代码演示来理解其工作原理，为后续深入学习打下基础。

傅里叶变换的直观类比

为了引入傅里叶变换背后的概念，我们使用一个类比。这个类比借鉴了平克·弗洛伊德乐队《月之暗面》专辑的经典封面。

这个封面与傅里叶变换有何关系？关键在于类比。这里有一束复杂的光波（复合波形），它穿过一个棱镜（算法或机器），然后从棱镜中我们得到了不同颜色的光谱带（分解后的不同频率成分）。

傅里叶变换做的事情与此非常相似。其高层直觉是：我们有一个复杂的声音，然后使用傅里叶变换（类比中的棱镜），将这个复杂的声音分解成其组成的频率。

因为正如我们在早期视频中所学，复杂的声音是由许多不同的纯音叠加在一起构成的。这是你需要掌握的第一个直观理解。

从时域到频域的旅程

当我们使用傅里叶变换时，我们实际上是在进行一次从时域到频域的旅程。

下图展示了一个声音的波形图，其X轴代表时间。

然后我们应用傅里叶变换，经过一番“魔法”，我们得到了同一个声音在频域中的表示，即对声音进行频率分析。此时，X轴代表频率。频谱中的尖峰意味着该特定赫兹处的频率是原始声音的重要组成部分。

那么，我们是如何做到这一点的呢？这就是今天视频要探讨的核心问题。

深入理解傅里叶变换的过程

现在，让我们进一步深入，更好地理解傅里叶变换内部发生了什么。

当我们应用傅里叶变换时，我们进行了一系列操作。

以下是其核心步骤的概述：

  1. 比较信号与正弦波：我们将原始信号与一系列不同频率的正弦波进行比较。
  2. 计算幅度和相位：通过比较，我们得到两个关键值：幅度和相位。
  3. 幅度表示相似度：幅度告诉我们原始信号与特定频率的正弦波有多么相似。幅度越高，相似度越高。
  4. 相位表示对齐：相位则告诉我们正弦波在时间上的偏移量，我们稍后会详细讨论。

为了具体化这些直觉，让我们通过一些代码和可视化来观察这个过程。

代码演示：观察傅里叶变换

我们将分析一个C5音符的钢琴声。首先，我们加载音频文件并绘制其时域波形。

如上图所示，我们在时域中看到了声音的波形，X轴是时间，Y轴是振幅。声音在约1.4秒处结束。

接下来是关键步骤：使用SciPy库中的fft函数计算傅里叶变换。我们取结果的绝对值来获得幅度谱，并对应地设置频率轴（从0 Hz到采样率）。然后绘制频谱图。

在频谱中，我们可以看到在523 Hz附近有一个明显的峰值，这正好是C5音符的基频。此外，在1040 Hz、1560 Hz和2080 Hz附近还有几个较小的尖峰，这些是基频的谐波（泛音），分别是基频的2倍、3倍和4倍。这个轮廓是谐波声音的典型特征。

可视化比较：信号与正弦波

为了理解傅里叶变换如何工作，我们需要放大波形，观察其周期，并将其与正弦波进行比较。

我们首先创建一个频率为523 Hz（基频）的正弦波，其相位初始设为0。正弦波的公式为：
s(t) = sin(2π * f * t - p)
其中 f 是频率，t 是时间，p 是相位。

上图展示了原始信号（蓝色）和初始正弦波（红色）。可以看到它们的频率相似，但并未对齐。

相位的作用：相位 p 控制正弦波在时间轴上的左右移动。当相位从0变化到1时，正弦波会经历一个完整的周期偏移。例如，相位为0.5时，两个正弦波完全反相（一个波峰对应另一个波谷）。

为了找到最佳对齐，我们调整正弦波的相位。通过尝试，我们发现当相位约为0.555时，正弦波与原始信号对齐得非常好。

计算相似度：点乘与面积

如何量化两个信号的相似度？一个直观简单的方法是：将两个信号逐点相乘（点乘），然后计算乘积曲线下的净面积（正面积减去负面积）。

我们将对齐后的正弦波与原始信号相乘，得到黄色曲线。

直觉：当两个信号同号（都为正或都为负）时，乘积为正，贡献正面积；当它们异号时，乘积为负，贡献负面积。净面积越大，说明两个信号在形状和相位上越相似。

如果我们改变相位（例如回到0），乘积曲线下几乎全是负面积，净面积很小，表明相似度很低。

因此，对于一个给定频率，傅里叶变换的核心步骤可以总结为：

  1. 选择该频率并生成一个正弦波。
  2. 优化相位，使得该正弦波与原始信号点乘后的净面积最大。
  3. 这个最大的净面积就是该频率对应的幅度（相似度度量）。

扩展到所有频率

傅里叶变换不仅仅对一个频率进行上述操作，而是对所有可能的频率（在实数范围内）重复这一过程。

对于每一个频率，我们都找到使其与原始信号最匹配的相位，并计算最大相似度（幅度）。最终，我们将所有频率及其对应的幅度绘制出来，就得到了频谱图。频谱中的高峰就对应了原始信号中能量最强的频率成分。

回到我们的频谱图，523 Hz处的峰值正反映了我们之前观察到的现象：频率为523 Hz的正弦波与原始信号具有很高的相似性。

数学形式化

上述直觉可以用数学公式进行形式化描述。对于连续信号，我们寻找对于给定频率 f 的最优相位 φ，以最大化点乘的积分（净面积）。

寻找最优相位的公式可以表示为：
φ_opt(f) = argmax_{φ ∈ [0,1)} ∫ s(t) * sin(2πft - φ) dt
其中 s(t) 是原始信号，积分计算的是点乘后的净面积，argmax 表示寻找使该积分最大的 φ 值。

一旦找到最优相位 φ_opt(f)，该频率的幅度 M(f) 就可以通过将最优相位代入来计算：
M(f) = ∫ s(t) * sin(2πft - φ_opt(f)) dt

这本质上是将我们手动调整相位、计算面积的过程用数学语言表达了出来。需要注意的是，这里使用的是连续时间信号的积分表示。在实际的数字信号处理中，我们处理的是离散采样信号，使用的是离散傅里叶变换，其数学形式会有所不同，但核心思想一致。

逆傅里叶变换：从频域回到时域

我们不仅可以从时域变换到频域，还可以反过来，从频域重建时域信号。这个过程称为逆傅里叶变换。

其思想是：将我们从原始信号中提取出的所有频率的正弦波，按照它们各自的幅度（作为权重）和相位叠加（相加）起来，就能重构出原始信号。

这意味着时域信号和其傅里叶变换（频域表示）包含完全相同的信息，我们可以在两者之间自由转换。这类似于加法合成，在合成器中，通过叠加不同频率、幅度和相位的正弦波来创造复杂的声音。

下图展示了一个加法合成器的例子，通过激活基频（C4， 261 Hz）及其谐波（C5， 523 Hz 等），可以合成出丰富的音色。

总结与展望

本节课中我们一起学习了傅里叶变换的直观理解。我们通过棱镜的类比引入了频率分解的概念，详细讲解了傅里叶变换通过将信号与不同频率的正弦波进行比较、计算相似度（幅度）和相位偏移来工作的过程，并通过代码可视化加深了理解。我们还简要了解了其数学形式以及逆变换的可能性。

傅里叶变换是音频和信号处理的基石。为了深入其数学原理和实现，我们需要先掌握另一个重要工具：复数。复数可以非常优雅地表示我们刚才讨论的幅度和相位。在接下来的视频中，我们将学习复数，为最终攻克傅里叶变换的数学核心做好准备。

011：音频信号处理中的复数

在本章中，我们将学习复数的基本概念。这些知识对于深入理解傅里叶变换至关重要。我们将探讨复数的起源、表示方法及其与音频信号处理的联系。

上一节我们介绍了傅里叶变换的直观理解，并提到需要复数知识才能更进一步。本节中，我们将深入探讨复数的核心概念。

为什么需要复数？🤔

在处理傅里叶变换时，每个纯音频率分量会得到两个参数：幅度和相位。幅度是实数。如果只关心幅度谱，则无需复数。但如果想完整理解傅里叶变换并处理相位信息，就需要一种能同时处理这两个参数的数学工具。复数正是这样的工具。

复数的起源

长期以来，数学家们被一个简单问题所困扰：负数的平方根。在实数范围内，无法求解如 √(-1) 这样的问题。为此，数学家引入了虚数单位 i，其定义为：

i² = -1

这使得处理负数的平方根成为可能。虽然这看起来像是“编造”的，但复数在信号处理等领域被广泛应用，是真实有效的数学工具。

复数的定义与表示

一个复数 C 通常定义为：

C = a + i * b

其中 a 和 b 都是实数。a 称为实部，b 称为虚部。

为了可视化复数，我们使用复平面。复平面的横轴（x轴）是实轴，纵轴（y轴）是虚轴。复数 a + i*b 可以表示为复平面上的一个点，其坐标为 (a, b)。

例如，复数 3 + 2i 对应点 (3, 2)；复数 -1 + 4i 对应点 (-1, 4)。

复数的极坐标表示

在音频信号处理中，极坐标表示法更为方便。它将复数表示为距离和角度，这与傅里叶变换中的幅度和相位概念直接对应。

从笛卡尔坐标 (a, b) 转换到极坐标需要两个参数：

  1. 模（绝对值） |C|：复数到原点的距离。
  2. 辐角 γ：正实轴与连接原点和复数的线段之间的夹角。

以下是计算这两个参数的公式。

计算模 |C|：
根据勾股定理，在复平面上，a 和 b 构成直角三角形的两条直角边，|C| 是斜边。因此：

|C| = √(a² + b²)

计算辐角 γ：
利用三角函数关系：

cos(γ) = a / |C|
sin(γ) = b / |C|

两式相除可得：

tan(γ) = sin(γ) / cos(γ) = b / a

因此，辐角 γ 可以通过反正切函数求得：

γ = arctan(b / a)

得到模和辐角后，我们可以将复数用极坐标形式重新表示。利用三角恒等式 a = |C| * cos(γ) 和 b = |C| * sin(γ)，代入复数定义：

C = a + i*b = |C| * cos(γ) + i * |C| * sin(γ) = |C| * [cos(γ) + i * sin(γ)]

这就是复数的极坐标形式。

欧拉公式与复指数表示

欧拉公式建立了复指数与三角函数的美妙联系：

e^(i * γ) = cos(γ) + i * sin(γ)

其中 e 是自然对数的底数。这个公式意味着，e^(i * γ) 在复平面上代表单位圆（半径为1的圆）上的一个点。随着角度 γ 增加，这个点沿单位圆逆时针移动。

利用欧拉公式，我们可以将复数的极坐标形式写得极其简洁：

C = |C| * e^(i * γ)

这个表示法非常强大：|C| 决定了复数的“长度”（距离原点的远近），而 e^(i * γ) 决定了它的“方向”（在复平面上的角度）。

一个著名的特例：欧拉恒等式
当 γ = π 时，代入欧拉公式：

e^(i * π) = cos(π) + i * sin(π) = -1 + i * 0 = -1

由此得到被誉为“数学中最优美公式”的欧拉恒等式：

e^(i * π) + 1 = 0

它巧妙地将数学中五个最重要的常数（e, i, π, 1, 0）联系在了一起。

可视化复数的两个组成部分

理解公式 C = |C| * e^(i * γ) 的直观含义非常重要：

  • 指数部分 e^(i * γ)：它像一个方向指针，指向单位圆上的某个角度 γ。它决定了复数在复平面上的朝向。
  • 模 |C|：它像一个缩放因子，将单位圆上的点拉近或推远原点。如果 |C| > 1，则点被推向圆外；如果 |C| < 1，则点被拉向圆内。

这种理解方式将抽象的公式转化为清晰的几何图像。

复数与傅里叶变换的联系🔗

现在，我们可以清晰地看到复数表示法与傅里叶变换参数之间的对应关系：

  • 复数的模 |C| 直接对应傅里叶变换中纯音分量的幅度（Magnitude）。
  • 复数的辐角 γ 直接对应傅里叶变换中纯音分量的相位（Phase）。

因此，复数 C = |C| * e^(i * γ) 以一种紧凑而优雅的方式，同时封装了音频信号中某个频率成分的强度和时间偏移信息。这为我们在复数域中理解和操作傅里叶变换奠定了坚实的基础。

总结

本节课中我们一起学习了复数的核心知识。我们了解了复数因解决负数开方问题而被引入，并掌握了其两种主要表示法：笛卡尔坐标形式 a + i*b 和极坐标形式 |C| * e^(i * γ)。我们推导了两种形式间的转换公式，并借助欧拉公式得到了极其简洁的复指数表示。最重要的是，我们建立了复数与傅里叶变换参数（幅度和相位）之间的直观联系，这将是后续深入理解傅里叶变换复数形式的钥匙。

012：用复数定义傅里叶变换 🎵

在本节课中，我们将学习如何使用复数来定义傅里叶变换。我们将看到，这种表示方法非常紧凑和优雅，它将傅里叶变换输出的幅度和相位两个参数，巧妙地编码进一个复数中。

上一节我们介绍了复数，并提到它将用于定义音频信号处理的许多方面。本节中，我们将在此基础上，使用复数来定义傅里叶变换。

回顾：傅里叶变换的参数与复数表示

首先，我们回顾一下之前介绍过的几个关键公式和概念。

以下是傅里叶变换中提取相位和幅度的公式：

相位公式：
φ = arctan2(b, a)

幅度公式：
M = sqrt(a² + b²)

这些公式在介绍傅里叶变换直观理解的视频中已经详细讲解过。在高层次上，相位和幅度是当我们用傅里叶变换分解一个复杂声音时，为每个频率提取的两个参数。幅度告诉我们原始信号中有多少该频率的成分。

接下来，我们回顾复数的极坐标表示法：

复数极坐标公式：
c = |c| * e^(iγ)

其中 i 是虚数单位，|c| 是复数的模（绝对值），γ 是角度。这是用指数形式表示的复数极坐标定义。

核心思想：将傅里叶系数编码为复数

现在，让我们进入正题。用复数定义傅里叶变换的直观思想是：每次应用傅里叶变换，我们为每个频率得到一对参数（幅度和相位）。我们可以将这对参数作为极坐标，来定义一个复数。换句话说，我们可以将这两个系数编码进一个复数中。

让我们看看这在数学上是如何实现的。

我们从之前提到的数学公式开始。我们有相位 φ 和幅度 M，以及复数的极坐标定义。

那么，我们如何得到傅里叶变换系数呢？公式如下：

复数傅里叶变换系数公式：
C_f = (M / sqrt(2)) * e^(-i * 2πφ)

这个复数被称为复数傅里叶变换系数。我们为分解原始信号得到的每个频率都有一个这样的系数。

仔细观察这个公式，我们可以将这个复数系数映射到复数的极坐标表示上：

  • 复数的模 |C_f| 对应 M / sqrt(2)。这里的 sqrt(2) 是由于归一化原因，它是一个常数因子。
  • 复数的角度 γ 对应 -2πφ。其中的负号 - 意味着当相位 φ 增加时，我们在复平面上顺时针旋转。

可视化：复平面上的傅里叶系数

现在，让我们回到复数傅里叶系数的可视化表示。

这个系数 C_f 不过是复平面上的一个点。在复平面中，水平轴是复数的实部，垂直轴是复数的虚部。

C_f 这个点具有以下几何意义：

  • 该点到原点的距离由系数的模 |C_f| 给出，即 M / sqrt(2)。
  • 该点与原点连线和正实轴之间的夹角是 γ，即 -2πφ。由于相位 φ 介于 0 和 1 之间，我们无法画出整个圆。

傅里叶变换的数学定义与视觉解读

现在，我想展示如何实际得到傅里叶变换，以及其数学定义。为此，我们从连续音频信号开始。

我们的连续音频信号是 g(t)。这是一个函数，输入一个实数（时间），输出另一个实数（振幅或强度）。这是一个波形。

这与复数傅里叶变换形成鲜明对比。复数傅里叶变换我们用 ĝ(f) 表示。它的输出是我们几分钟前看到的系数，即傅里叶变换系数。我们为每个频率都有一个这样的复数。

这个函数 ĝ(f) 有所不同：它输入一个实数（频率，单位赫兹），但输出一个复数，即傅里叶变换系数。

在复平面上看，对于频率 f1，复数傅里叶变换的输出可以可视化为复平面上的一个点。这非常有趣，你现在可能开始理解为什么我们使用复数表示傅里叶变换了。它的伟大之处在于，傅里叶变换产生两个不同的参数（幅度和相位），我们可以将它们打包成一个方便的工具——复数，并轻松地在复平面上可视化。

但你可能会问，我们如何从数学上得到这个系数 C_f 呢？

让我们再次从之前视频中看过的原始公式开始。我们有幅度和相位从傅里叶变换中产生的公式，以及复数傅里叶变换的定义。

我知道下面这个积分公式看起来有点吓人，但请耐心听我解释，因为我们将查看它的视觉解读，你会发现它其实很容易理解。

复数傅里叶变换积分公式：
ĝ(f) = ∫ g(t) * e^(-i*2πft) dt

如前所述，ĝ(f) 在特定频率 f 的输出不过是复平面上的一个复数，即我们的傅里叶变换系数。但我们如何从这个包含指数和复数的“怪兽”积分中得到它呢？

首先，我们来看看这个公式的不同部分。

第一部分是指数项 e^(-i*2πft)。这个指数在复平面上描绘单位圆。随着 t 增加，我们沿着整个单位圆移动。由于负号 - 的存在，我们是顺时针描绘这个圆。需要记住的另一点是，我们完成一个圆的速度取决于频率值 f。例如，如果 f = 1 Hz，则需要 1 秒完成一个整圆；如果 f = 2 Hz，则只需 0.5 秒（因为周期 T = 1/f）。这本质上就是用复数表示正弦波（纯音），我们将用它来分解原始信号。

第二部分是我们的信号 g(t)。为了展示视觉解读，我创建了一个自定义信号（由基频及其二次、三次谐波的正弦波叠加而成），并在 Jupyter Notebook 中进行了绘制。

第三部分是信号 g(t) 与纯音 e^(-i*2πft) 的乘法。这在视觉上是什么样子呢？

当我们绘制这个乘法结果时，我们看到一个形状出现在复平面上。正在发生的是，我们通过将原始信号乘以一个纯音，将原始信号“缠绕”在了复平面上。 这非常酷。

如果我稍微改变纯音频率，比如从 1 Hz 改为 1.1 Hz，输出形状会发生巨大变化。为了理解这一点，我们需要查看时域。如果我们只观察 1 秒内的数据，会看到一个形状；在下一秒，会绘制另一个形状。继续下去，直到 10 秒，我们基本上旋转了一整圈。

另一个有趣的现象是，如果我用频率 2 Hz 或 3 Hz（这些频率存在于原始信号中）的纯音去分解，会得到相当稳定的形状。但如果用 3.1 Hz（不存在于信号中），又会得到那种混乱的形状。直观上，这意味着每当原始声音中存在的频率成分被用来分解时，我们就会得到一个更稳定的输出。

理解积分：求“重心”与求和

现在我们有了信号与纯音相乘的视觉解读。但这个积分呢？它看起来像什么？

直观上，积分所做的是：取所有不同的点（不同时间点对应的复数）并对它们求和。在物理系统中，这类似于求该系统的重心（即所有点的平均值）。

当我绘制这个重心（用红点表示）时，例如用 3 Hz 的纯音分解，重心在某个位置。用 1 Hz（基频）分解，重心也在某个位置。但如果用 1.1 Hz 分解，重心就在原点。这并不意外，因为图形完全对称，所有点相互抵消，平均值为 0。

这与我们对傅里叶变换的直观理解产生了共鸣：当原始信号中存在某个频率成分时，其幅度不为零；当处理的频率成分不存在时，幅度为零。

现在，关键点来了：我们可以将复数傅里叶变换系数看作是复平面上的这个点（重心）。但这里我有点“作弊”，因为我们并不是在求平均值，而是在求和。积分告诉我们对所有值求和。

那么，如果我绘制所有这些点的和（用绿点表示），会发生什么？例如，在 1000 个时间步长下，和可能是一个很大的数（如 -500）。如果只取 100 步，和大约是 -50。你能看出这与重心的联系吗？重心大约在虚轴 -0.5 的位置。和（-50）大约是重心（-0.5）乘以步数（100）的结果。所以，概念上，我们取重心，然后乘以所考虑的时间步数。

对于 1.1 Hz 的情况，由于图形对称，所有点相互抵消，和就在原点。这就是我们如何在原始声音中实际存在的频率上得到那些尖峰的方法。

现在，你应该对如何得到傅里叶变换系数有了视觉上的理解：

  1. 我们从一个信号开始。
  2. 通过将信号乘以一个纯音，将信号缠绕在复平面上。
  3. 然后，我们对每个时间步的结果求和。

但你可能会问，为什么这里用积分符号 ∫，而不是求和符号 ∑？这取决于我们使用的是连续表示。如果我们考虑的时间步长无限多（即连续时间），那么求和就变成了积分。

幅度、相位与逆傅里叶变换

回到复系数 C_f 的定义，我们可以写出其实部和虚部：

ĝ(f) = a + bi

我们可以取傅里叶变换的幅度，即取该函数的绝对值：

|ĝ(f)|

并且我们知道，|ĝ(f)| 等于 M / sqrt(2)。同样，sqrt(2) 是归一化常数。

对于相位 φ，我们可以通过取 ĝ(f) 的角度 γ 来定义：

φ = -γ / (2π)

现在，你已经对傅里叶变换有了很好的深入理解。你知道幅度和相位是什么，如何将它们编码成一个称为傅里叶变换系数的复数，也知道了获取傅里叶变换的数学公式，并且理解了其视觉解读。

还缺少什么呢？缺少的是逆傅里叶变换。我们在之前的视频中已经见过，可以从频域回到时域。逆傅里叶变换的数学表示如下：

逆傅里叶变换公式：
g(t) = ∫ ĝ(f) * e^(i*2πft) df

高层次上的直观理解是：我们取所有不同的频率成分，用幅度对它们进行加权，加上相位，然后把它们全部加起来，从而重建原始信号 g(t)。关键在于，原始信号和它的傅里叶变换包含相同的数据，我们可以在两者之间任意转换而不丢失任何信息。

让我们快速看一下视觉解读：

  1. e^(i*2πft) 是频率为 f 的纯音（正弦波）。
  2. 我们将复数傅里叶变换系数 ĝ(f) 乘以这个纯音。这实际上是用幅度对纯音进行加权，并加上相位，从而得到一个可能具有不同相位和振幅的新信号。
  3. 最后，我们积分（即求和）。这里不是对时间积分，而是对频率积分。换句话说，我们将所有加权的正弦波相加，得到原始声音。

傅里叶变换的“往返之旅”

现在，我们准备好看看傅里叶变换的“往返之旅”了。

傅里叶变换：
ĝ(f) = ∫ g(t) * e^(-i*2πft) dt

逆傅里叶变换：
g(t) = ∫ ĝ(f) * e^(i*2πft) df

它们基本上是同一个函数，唯一真正的区别是：在傅里叶变换中，我们对时间积分；在逆傅里叶变换中，我们对频率积分。这非常优美。所有的复杂性都被一个非常直接和极简的公式优雅地捕捉到了，并且在时域和频域之间来回转换的想法也非常美妙。

我希望你现在能理解为什么我们要不厌其烦地引入复数，因为通过依赖这个数学构造，我们现在有了傅里叶变换及其逆变换的超级紧凑的定义。

总结与展望

本节课内容很长，但我希望对你来说非常值得。现在，我们终于破解了将复杂声音分解成频率成分以及反向合成的“魔法”。我们了解了傅里叶变换。

然而，我们还有一个问题。到目前为止，我们处理傅里叶变换的所有理论都是在连续数学的框架内，即我们假设处理的是模拟或连续信号。但在现实中，当我们实际应用傅里叶变换时，我们处理的是数字信号、离散信号。我们从模拟信号开始，然后对其进行采样和数字化。

因此，关于傅里叶变换，剩下的工作就是理解如何将我们迄今为止所学的理论适应到离散世界，即我们的数字计算机世界，也是我们实际从信号中提取信息所要使用的世界。下一次，我们将学习离散傅里叶变换。

这就是今天的全部内容。我知道信息量很大，但希望你觉得值得。现在你应该对傅里叶变换有了非常深入的理解。如果你有任何问题，欢迎在下方评论区留言。我们下次再见！

013：离散傅里叶变换详解 🎵

在本节课中，我们将要学习连续傅里叶变换如何演变为其离散版本——离散傅里叶变换。我们将理解为何需要这种转变，以及如何通过数学调整来处理数字信号。课程将涵盖从连续到离散的转换过程、DFT的最终公式、其可视化解释以及一个关键的高效计算算法。

从连续到离散的转变

上一节我们介绍了使用复数表达的连续傅里叶变换。本节中我们来看看如何处理数字世界中的离散信号。

连续傅里叶变换公式适用于时间上连续的模拟信号。然而，在实际应用中，我们处理的是由计算机等数字设备处理的离散信号。

从连续信号 G(t) 转换到离散信号 x[n] 的过程称为数字化。这包括采样和量化。采样是指在特定时间点（采样点）记录信号的值，采样周期 T 定义了采样点之间的时间间隔。

离散信号 x[n] 表示在第 n 个采样点获取的值，其对应的实际时间为 t = n * T。

构建离散傅里叶变换

现在，让我们尝试从连续傅里叶变换出发，构建离散傅里叶变换。我们将看到公式中的每个元素如何映射到其离散对应物。

我们从傅里叶变换的定义开始。在连续情况下，我们称之为 Ĝ(f)；在离散情况下，我们称之为 x̂(f)，以表明我们处理的是样本。

以下是核心元素的转变：

  1. 积分变为求和：连续时间意味着无限的时间点，因此使用积分。离散情况下，我们只有有限的样本点，因此使用求和符号 Σ，对样本索引 n 进行求和。
  2. 连续时间 t 变为离散索引 n：在复指数项中，我们用离散的样本索引 n 替换连续时间变量 t。

由此，我们得到离散傅里叶变换的初步版本：

x̂(f) = Σ (x[n] * e^(-i2πfn)) （对所有 n 求和）

DFT的可视化解释

为了直观理解DFT，我们可以参考其可视化解释。

在连续傅里叶变换中，我们计算信号曲线下的面积。在离散版本中，我们近似这个面积。我们将每个样本点视为一个矩形的顶部，矩形的宽度由采样周期决定。DFT的结果近似于所有这些矩形面积的总和，它代表了原始连续信号下面积的离散近似。

采样周期越小（采样率越高），矩形数量越多，近似误差就越小。

完善DFT公式：解决两个问题

我们初步的DFT公式仍存在两个问题，使其无法直接用于数字计算。

以下是这两个问题及其解决方案：

  1. 问题一：无限数量的样本：理论上，傅里叶变换需要对无限时间（即无限样本）进行积分。
    • 解决方案（技巧一）：我们只考虑信号在有限时间区间内的样本。例如，对于一首3分钟的歌曲，我们只分析这3分钟内的样本。这假设该时间段包含了我们所需的所有频率信息。因此，求和限制在有限的 N 个样本上。
  2. 问题二：连续的频率：频率 f 在理论上是连续变量，意味着有无限多个频率分量。
    • 解决方案（技巧二）：我们只计算有限个频率点上的变换。一个巧妙且高效的做法是：考虑的频率数量 M 等于信号样本的数量 N，即 M = N。这样做有两个主要原因：一是保证了从时域到频域再通过逆变换返回时域的完整“往返”能力；二是计算上更高效。

应用这两个技巧后，我们得到了最终广泛使用的离散傅里叶变换公式。

离散傅里叶变换的最终公式

应用上述调整后，我们得到离散傅里叶变换的标准定义：

X[k] = Σ (x[n] * e^(-i2πkn/N)) （对 n=0 到 N-1 求和）

其中：

  • N 是信号的总样本数。
  • n 是当前样本的索引（从0到N-1）。
  • k 是当前频率分量的索引（从0到N-1）。
  • X[k] 是第 k 个频率分量的复数傅里叶系数。

这里，我们用 k/N 替代了连续的频率 f。k 的取值范围是 0 到 N-1，对应着 N 个离散的频率点。

第 k 个索引对应的实际频率（以赫兹Hz为单位）由以下公式给出：

频率 = k / (N * T) = k * Fs / N

其中 T 是采样周期，Fs 是采样率（Fs = 1/T）。

这意味着DFT输出的频率范围是从 0 Hz（当 k=0）到采样率 Fs（当 k=N 时，但实际 k 最大为 N-1，接近 Fs）。这个范围被均匀地划分为 N 个点。

DFT幅度谱与奈奎斯特频率

DFT的输出 X[k] 是复数，包含每个频率分量的幅度和相位信息。对于音频分析，幅度谱尤为重要，它显示了每个频率在信号中的强度。

观察DFT的幅度谱时，会发现一个有趣现象：频谱在频率 Fs/2 处呈现对称性。Fs/2 这个频率称为奈奎斯特频率。

  • 频率轴从 0 到 Fs。
  • 从 0 到 Fs/2 的部分是“真实”的频率信息。
  • 从 Fs/2 到 Fs 的部分是前半部分的镜像复制，是冗余的。

因此，在分析时，我们通常只关心从 0 到 Fs/2 的频率范围。奈奎斯特频率是数字系统能够无失真重建的最高频率，高于此频率会发生混叠现象。

快速傅里叶变换

直接计算DFT需要进行大约 N² 次运算，当样本数 N 很大时计算量非常昂贵。

快速傅里叶变换 是一种高效计算DFT的算法。它通过利用正弦和余弦计算中的对称性和周期性冗余，显著减少了运算量。

FFT将计算复杂度从 O(N²) 降低到 O(N log₂ N)。为了充分发挥FFT的效率，通常要求样本数 N 是2的幂（如256, 512, 1024等）。现代所有的傅里叶变换实现（如在NumPy、SciPy或音频处理软件中）都使用FFT算法。

本节课中我们一起学习了离散傅里叶变换。我们从连续傅里叶变换出发，理解了为何及如何将其调整为处理离散数字信号的版本。我们探讨了DFT的构建过程、最终公式、其输出的频率含义、幅度谱的对称性（奈奎斯特频率的关键作用），以及用于高效计算的快速傅里叶变换算法。下一节，我们将暂时抛开理论，进行实际编程，从声音中提取幅度谱并分析其语义含义。

014：使用Python提取傅里叶变换 🎵

概述

在本节课中，我们将学习如何使用Python提取音频信号的傅里叶变换系数，并绘制其幅度谱。我们将通过对比不同乐器演奏同一音符的频谱，来初步理解音色差异的成因。

上一节我们介绍了离散傅里叶变换和快速傅里叶变换的理论基础。本节中，我们将动手实践，使用Python代码来实现这些概念。

准备工作

首先，我们需要导入必要的Python库。

以下是所需的库及其用途：

  • numpy：用于数值计算，特别是执行快速傅里叶变换。
  • matplotlib.pyplot：用于绘制频谱图。
  • librosa：用于加载音频文件。
  • IPython.display：用于在Jupyter笔记本中直接播放音频。

import numpy as np
import matplotlib.pyplot as plt
import librosa
import IPython.display as ipd

接下来，我们定义音频文件所在的目录路径，并加载四个音频文件：小提琴、萨克斯、钢琴演奏的C音，以及一段噪声。

base_dir = ‘/home/valerio/pycharm_projects/audio_processing_for_ML/14/audio/‘
violin_file = ‘violin_C4.wav‘
sax_file = ‘sax_C4.wav‘
piano_file = ‘piano_C5.wav‘
noise_file = ‘noise.wav‘

现在，我们可以在笔记本中直接播放这些音频片段，以便了解它们的声音。

ipd.Audio(base_dir + violin_file)
ipd.Audio(base_dir + sax_file)
ipd.Audio(base_dir + piano_file)
ipd.Audio(base_dir + noise_file)

加载音频数据

为了进行信号处理，我们需要将音频文件加载为NumPy数组。我们将使用librosa.load函数，它会返回音频波形数据及其采样率。

violin_C4, sr = librosa.load(base_dir + violin_file)
sax_C4, _ = librosa.load(base_dir + sax_file)
piano_C5, _ = librosa.load(base_dir + piano_file)
noise, _ = librosa.load(base_dir + noise_file)

加载后，我们可以查看小提琴信号的形状，它代表时域中的样本数量。

print(violin_C4.shape)

提取傅里叶变换系数

现在，我们使用NumPy的fft.fft函数来计算信号的快速傅里叶变换。这将把信号从时域转换到频域。

violin_ft = np.fft.fft(violin_C4)

变换后得到的数组与原始信号具有相同的长度，每个元素对应一个频率仓。每个元素都是一个复数，包含该频率分量的幅度和相位信息。

print(violin_ft[0]) # 查看第一个频率仓的复数系数

计算幅度谱

在音频AI应用中，我们通常更关心频率的幅度而非相位。因此，我们通过取复数系数的绝对值来获得幅度谱。

violin_magnitude_spectrum = np.abs(violin_ft)

此时，幅度谱中的值变为实数，代表了每个频率分量的能量大小。

print(violin_magnitude_spectrum[0]) # 查看第一个频率仓的幅度

可视化幅度谱

为了直观比较不同声音的频谱，我们创建一个函数来绘制幅度谱。

以下是plot_magnitude_spectrum函数的关键步骤：

  1. 计算信号的傅里叶变换。
  2. 计算幅度谱。
  3. 生成对应的频率轴（0 Hz 到采样率）。
  4. 使用matplotlib绘制频率与幅度的关系图。
def plot_magnitude_spectrum(signal, title, sample_rate, freq_ratio=1):
    ft = np.fft.fft(signal)
    magnitude_spectrum = np.abs(ft)
    
    frequency = np.linspace(0, sample_rate, len(magnitude_spectrum))
    num_freq_bins = int(len(frequency) * freq_ratio)
    
    plt.figure(figsize=(10, 5))
    plt.plot(frequency[:num_freq_bins], magnitude_spectrum[:num_freq_bins])
    plt.xlabel(‘Frequency (Hz)‘)
    plt.title(title)
    plt.show()

参数freq_ratio允许我们只查看频谱的一部分。例如，设置为0.5表示只显示到奈奎斯特频率（采样率的一半），设置为0.1可以放大观察低频部分。

对比不同乐器的频谱

现在，我们使用这个函数来绘制并对比四个音频文件的幅度谱。

首先绘制小提琴演奏C4的频谱，并放大观察前10%的频率范围。

plot_magnitude_spectrum(violin_C4, ‘Violin C4‘, sr, freq_ratio=0.1)

在频谱图中，我们可以在约260 Hz处看到一个明显的峰值，这是C4音符的基频。此外，在520 Hz（第一谐波）等处也有峰值。有趣的是，小提琴的第一谐波能量可能比基频更强。

接下来绘制萨克斯演奏C4的频谱。

plot_magnitude_spectrum(sax_C4, ‘Saxophone C4‘, sr, freq_ratio=0.1)

小提琴和萨克斯都演奏C4，因此它们的频谱在基频和谐波位置相似。然而，各频率分量能量的相对比例不同，这正是我们感知到它们音色不同的主要原因之一。

然后绘制钢琴演奏C5的频谱。

plot_magnitude_spectrum(piano_C5, ‘Piano C5‘, sr, freq_ratio=0.1)

钢琴演奏的是C5，其基频约为520 Hz（是C4的两倍）。钢琴的频谱特征通常是基频能量最强，谐波能量依次衰减，这与弦乐器有所不同。

最后，我们观察噪声信号的频谱。

plot_magnitude_spectrum(noise, ‘Noise‘, sr, freq_ratio=0.1)

噪声的频谱显示，能量广泛分布在几乎所有频率仓中，没有像乐音那样清晰的峰值结构。

总结

本节课中我们一起学习了如何使用Python提取和可视化音频信号的傅里叶变换。我们掌握了以下核心操作：

  • 使用np.fft.fft计算快速傅里叶变换。
  • 使用np.abs从复数系数中获取幅度谱。
  • 创建函数绘制并对比不同声音的幅度谱。

通过对比小提琴、萨克斯、钢琴和噪声的频谱，我们看到了谐波结构如何影响音色，以及噪声在频域中的能量分布特性。

需要指出的是，我们目前绘制的幅度谱是整个音频片段在时间上的平均。它丢失了频率成分随时间变化的信息。在下一节课中，我们将介绍短时傅里叶变换的理论，它将使我们能够提取随时间变化的频谱图，这将是音频AI和音乐信息检索中分析数字音频信号的主要工具。

015：短时傅里叶变换详解 🎵

在本节课中，我们将学习音频信号处理中的一个核心概念——短时傅里叶变换。我们将了解它的工作原理、数学表达、关键参数，以及它如何生成对机器学习至关重要的特征——频谱图。

概述

上一节我们使用Python和Librosa从音频文件中提取了离散傅里叶变换。本节我们将回到理论层面，探讨AI音频中的一个关键主题：短时傅里叶变换。短时傅里叶变换之所以重要，是因为它使我们能够提取频谱图，而频谱图可能是可以馈送给深度学习音频模型的最重要的特征。

从离散傅里叶变换到短时傅里叶变换

首先，让我们回顾一下离散傅里叶变换。其数学公式如下：

X[k] = Σ_{n=0}^{N-1} x[n] * e^{-j2πkn/N}

其高层次直觉是：我们从时域信号（如波形）开始，应用离散傅里叶变换，得到的结果是原始信号中不同频率成分存在情况的一幅“图片”，即幅度谱。然而，这是一个“静态图片”，因为它只提供了一个平均了整个信号持续时间内频率成分存在情况的单一图像。这里存在一个问题：我们知道信号中存在哪些频率成分，但我们不知道它们何时出现得更多或更少，因为所有信息都被平均了。

音频数据是关于频率成分随时间演变的动态过程。我们希望了解不同频率成分如何随时间变化，这正是短时傅里叶变换的核心目的——从静态图像转向能提供跨时间频率成分信息的“视频”。

短时傅里叶变换的核心思想

短时傅里叶变换的高层次思想是：我们不对整个信号持续时间进行傅里叶变换，而是考虑信号的小片段或块，技术上我们称之为“帧”，然后对每一帧应用离散傅里叶变换。

为了更直观地理解，让我们将其可视化。我们从一个音频信号开始，然后只考虑第一块（第一帧）。此时，仅对该帧内的样本应用离散傅里叶变换，得到幅度谱。接着，我们滑动到下一帧，再次对该帧应用离散傅里叶变换，如此重复，直到遍历整个信号。

分帧与加窗

一种用于推导信号片段的方法是通过加窗。这意味着我们取原始信号，然后逐样本地乘以一个窗函数，从而得到一个加窗信号。

例如，我们从一个音频信号开始，然后应用一个矩形窗函数。矩形窗函数是一个红色曲线，除了一个片段内值为1外，其他地方均为0。将信号与矩形窗相乘，我们就得到了下方的加窗信号。

现在，介绍两个对今天讨论非常重要的参数：窗长 和 帧长。它们都以样本数衡量，但指代略有不同。

  • 窗长：我们应用加窗的样本数量。
  • 帧长：当我们分割信号并传递给短时傅里叶变换以计算每帧的傅里叶变换时，我们在信号的每个块中考虑的样本数量。

通常，窗长和帧长是一致的，具有相同的样本数。但有时帧长可能大于窗长。不过，在大多数应用场景中，窗长和帧长是相同的。例如，在Librosa中提取短时傅里叶变换时，我们不必强制传递窗长，其默认值就是帧长。

如果窗长小于帧长会发生什么？我们应用傅里叶变换的块仍然是整个帧，但加窗只发生在窗长数量的样本上，剩余的样本（帧长与窗长之差）将进行零填充。为了简化，本视频中我们假设窗长等于帧长。

此时，我们已经有了一个帧，并对其应用了窗函数。接下来就是应用离散傅里叶变换，从该帧中获取频率成分。然后我们向右滑动，得到第二帧，同样滑动窗函数并应用离散傅里叶变换，直到信号结束。

但这是一个简化版本。在实际的短时傅里叶变换中，帧通常是重叠的，第二帧与第一帧重叠。这就需要引入另一个参数：跳数。

跳数与帧重叠

跳数（通常用大写H表示）告诉我们当获取新帧时向右滑动多少样本。帧的重叠很重要，原因涉及频谱泄漏等问题。如果你想深入了解为什么跳数如此重要以及为什么需要重叠帧，建议观看关于音频特征提取管道的视频。

短时傅里叶变换的数学公式

现在，让我们从视觉直觉转向数学公式。我们将对比离散傅里叶变换和短时傅里叶变换的数学表达。

离散傅里叶变换的公式为：

X[k] = Σ_{n=0}^{N-1} x[n] * e^{-j2πkn/N}

短时傅里叶变换的公式为：

S[m, k] = Σ_{n=0}^{N-1} (x[n + mH] * w[n]) * e^{-j2πkn/N}

让我们逐项比较它们如何对应：

  1. 定义与变量：

    • 离散傅里叶变换 X[k] 是频率 k 的函数。每个公式为我们提供一个复数傅里叶系数，包含相位和幅度信息。
    • 短时傅里叶变换 S[m, k] 不仅依赖于频率 k，还依赖于 m。m 是时间的代理，代表当前所在的帧号。因此，短时傅里叶变换的结果是第 m 个时间仓（帧）中第 k 个频率的傅里叶系数，它仍然是一个包含相位和幅度信息的复数。

  2. 求和范围：

    • 在离散傅里叶变换中，我们对信号中的所有样本（整个持续时间）求和，N 等于信号中的总样本数。
    • 在短时傅里叶变换中，我们只对一个帧内的样本求和，N 等于帧长，因为我们考虑的不是整个信号，而只是一帧。

  3. 信号部分：

    • 在离散傅里叶变换中，x[n] 是整个信号。
    • 在短时傅里叶变换中，我们只考虑当前帧 m 中存在的信号。m * H 是当前帧的起始样本（H 是跳数），加上 n（从0到帧长-1），覆盖了该帧内的所有样本。然后，我们将该信号乘以窗函数 w[n]，得到该帧的加窗信号。

  4. 复指数项：

    • 最后一步对于两者是相同的：乘以一个频率为 k/N 的复指数（纯音）。这样做的本质是将信号分解并投影到该频率的纯音上。

短时傅里叶变换的输出

理解了数学公式，你可能会问：我们从离散傅里叶变换和短时傅里叶变换中具体得到什么输出？

  • 离散傅里叶变换：我们提取一个频谱向量，它是一个一维数组，包含一定数量的频率仓的傅里叶系数。这里没有时间维度，因为所有信息都是整个信号持续时间的平均。
  • 短时傅里叶变换：我们得到一个频谱矩阵，它是一个二维数组，包含频率仓数量和帧数量。换句话说，我们为考虑的每个频率仓和每一帧得到一个复数傅里叶系数。这样，我们不仅保留了频率信息，还通过不同的帧（时间的代理）获得了时间信息。

计算输出形状

我们可以计算短时傅里叶变换输出的频率仓数和帧数。

频率仓数：
公式为：帧长 / 2 + 1
原因在于离散傅里叶变换具有围绕奈奎斯特频率的镜像对称性。第一半部分包含信息，第二半部分是其镜像。因此，在短时傅里叶变换中，我们只取第一半部分的信息以避免冗余。

帧数：
公式为：(总样本数 - 帧长) / 跳数 + 1
建议你作为练习，通过可视化来理解这个公式。

示例：
假设信号有10000个样本，帧长为1000，跳数为500。

  • 频率仓数 = 1000 / 2 + 1 = 501
  • 帧数 = (10000 - 1000) / 500 + 1 = 19
    因此，短时傅里叶变换的输出形状为 (501, 19)，是一个二维数组。

短时傅里叶变换的关键参数

短时傅里叶变换的输出取决于我们传递的一系列参数。不同的参数会导致不同的输出。

1. 帧长
帧长决定了我们将原始信号分割成的块的大小，以样本数衡量。常用值是2的幂次方，如512、1024、2048。选择2的幂次方很重要，因为这样我们可以使用快速傅里叶变换来计算离散傅里叶变换，这是一种非常快速且计算高效的方法。

选择帧长时涉及一个重要的概念：时频权衡。

  • 如果使用较大的帧长，频率分辨率会增加，但时间分辨率会下降。因为更多的样本意味着更多的频率仓（更好的频率分辨率），但也意味着考虑的时间块更大（更差的时间分辨率）。
  • 反之，如果使用较小的帧长，频率分辨率会下降（更少的频率仓），但时间分辨率会提高（考虑更小的时间块）。

这个权衡无法完全解决，通常需要根据具体应用找到一个合适的折中值。例如，某些问题可能更需要高频率分辨率，而像起始点检测这样的应用可能更需要高时间分辨率。

2. 跳数
跳数是我们获取新帧时向右滑动的样本数。常用值也是2的幂次方，如256、512。也可以将其定义为帧长的一部分，例如帧长的一半、四分之一或八分之一。

3. 窗函数
短时傅里叶变换的结果也取决于所选的窗函数，因为不同的窗函数会以不同的方式调制原始信号。

虽然我们介绍了矩形窗，但在数字信号处理中它并不常用，因为它会在边缘产生不连续性。为了避免这些，通常使用钟形曲线窗，其中最重要的是汉宁窗。大约90%的情况下，你在执行短时傅里叶变换时可能都在使用汉宁窗。其公式如下：

w[n] = 0.5 * (1 - cos(2πn / (N-1))) , 0 ≤ n ≤ N-1

当我们将汉宁窗应用于信号时，信号被调制，样本值在边缘趋向于0，从而避免了不连续性。这有助于减少频谱泄漏。

从短时傅里叶变换到频谱图

最后，我们来到本视频的最终主题，也是你可能最关心的部分：频谱图。通过频谱图，我们可以可视化声音。

我们如何得到频谱图？到目前为止，我们知道短时傅里叶变换是一个矩阵，其中每个元素都是复数傅里叶系数。我们取短时傅里叶变换的幅度平方，得到一个与原始短时傅里叶变换形状相同的矩阵，但区别在于现在所有元素都是实数。然后，我们可以使用热图将其可视化，这种可视化就称为频谱图。

频谱图对AI音频的所有应用都至关重要，因为很多时候我们都会使用频谱图作为输入算法的特征。

观察一个频谱图：

  • X轴 代表时间，显示离散的时间点（帧/时间仓）。
  • Y轴 代表频率，显示不同的频率仓。
  • 图中显示的是不同频率成分如何随时间（跨帧）演变。

现在，我们梦想成真了：我们不仅拥有频率成分的信息（这是幅度谱已有的），还拥有了成分随时间演变的信息（这是时域通常提供的信息）。因此，频谱图被称为时频表示，这也是频谱图在AI音频中如此重要的原因。

总结

本节课我们一起学习了短时傅里叶变换。我们从回顾离散傅里叶变换的局限性开始，引出了需要分析频率随时间变化的需求。接着，详细讲解了短时傅里叶变换通过分帧、加窗、重叠的方式，将一维时域信号转换为二维时频表示的核心思想。我们深入探讨了其数学公式，并与离散傅里叶变换进行了对比。然后，我们学习了如何计算短时傅里叶变换输出的形状（频率仓数和帧数），并分析了影响其结果的关键参数：帧长（涉及时频权衡）、跳数和窗函数。最后，我们了解到通过对短时傅里叶变换结果取幅度平方，可以得到能够可视化声音动态频率内容的频谱图，这是机器学习音频应用中最常用的特征之一。

在下一讲中，我们将使用Python和Librosa具体提取频谱图，研究不同种类的频谱图，并比较不同音频样本（例如不同音乐流派）的频谱图。

016：使用Python从音频中提取频谱图

在本节课中，我们将学习如何使用Python和Librosa库从音频文件中提取并可视化频谱图。我们将从加载音频开始，逐步进行短时傅里叶变换，计算频谱图，并最终生成对数幅度和对数频率的频谱图以进行更好的可视化分析。

导入必要的库

首先，我们需要导入一些必要的Python库来处理音频、进行数学运算和可视化。

import os
import librosa
import librosa.display
import numpy as np
import matplotlib.pyplot as plt
import IPython.display as ipd

加载音频文件

我们将使用四个不同的音频文件进行演示：一个钢琴音阶、一段德彪西的古典音乐、一段红辣椒乐队的摇滚乐以及一段艾灵顿公爵的爵士乐。

以下是加载音频文件的代码。librosa.load函数返回音频信号（一个NumPy数组）和采样率。

# 定义音频文件路径
scale_path = ‘path/to/scale.wav’
debussy_path = ‘path/to/debussy.wav’
redhot_path = ‘path/to/redhot.wav’
duke_path = ‘path/to/duke.wav’

# 加载音频文件
scale_signal, sr_scale = librosa.load(scale_path)
debussy_signal, sr_debussy = librosa.load(debussy_path)
redhot_signal, sr_redhot = librosa.load(redhot_path)
duke_signal, sr_duke = librosa.load(duke_path)

为了对处理的音频有一个直观印象，我们可以在Jupyter笔记本中直接播放它们。

# 播放音阶音频
ipd.Audio(scale_signal, rate=sr_scale)

提取短时傅里叶变换

上一节我们介绍了短时傅里叶变换的理论。本节中，我们来看看如何使用Librosa库实际计算它。

首先，我们需要设置两个关键参数：帧大小和跳跃长度。

# 设置参数
frame_size = 2048  # 帧大小（样本数）
hop_size = 512     # 跳跃长度（样本数）

现在，我们可以使用librosa.stft函数来计算短时傅里叶变换。

# 计算音阶音频的短时傅里叶变换
scale_stft = librosa.stft(scale_signal, n_fft=frame_size, hop_length=hop_size)

让我们检查一下结果的形状和数据类型。

# 检查STFT结果的形状
print(scale_stft.shape)  # 输出: (1025, 342)

# 检查数据类型（应为复数）
print(type(scale_stft[0, 0]))  # 输出: <class ‘numpy.complex128’>

结果矩阵的第一维（1025）是频率仓的数量，等于frame_size/2 + 1。第二维（342）是时间帧的数量。

计算频谱图

短时傅里叶变换的结果是复数。为了得到可以可视化的频谱图，我们需要计算其幅度平方。

以下是计算频谱图的公式：
频谱图 = |STFT|²

在代码中，我们使用NumPy的绝对值函数和平方运算来实现。

# 计算幅度平方频谱图
scale_spectrogram = np.abs(scale_stft) ** 2

检查新矩阵的形状和数据类型。

print(scale_spectrogram.shape)  # 输出: (1025, 342)
print(type(scale_spectrogram[0, 0]))  # 输出: <class ‘numpy.float64’>

现在，我们有了一个实数矩阵，代表每个频率和时间点上的能量。

可视化频谱图（线性表示）

有了频谱图数据，我们现在可以尝试将其可视化。Librosa提供了一个方便的函数librosa.display.specshow。

首先，我们创建一个辅助函数来绘制频谱图。

def plot_spectrogram(Y, sr, hop_length, y_axis=“linear”):
    plt.figure(figsize=(10, 6))
    librosa.display.specshow(Y, sr=sr, hop_length=hop_length, x_axis=“time”, y_axis=y_axis)
    plt.colorbar(format=“%+2.0f dB”)
    plt.title(“Spectrogram”)
    plt.show()

现在，用这个函数绘制音阶的线性频谱图。

plot_spectrogram(scale_spectrogram, sr_scale, hop_size, y_axis=“linear”)

您可能会注意到图像非常暗，大部分区域是黑色的。这是因为人耳对声音强度的感知是对数式的，而非线性的。

转换为对数幅度频谱图

为了使频谱图更符合人耳的感知，我们需要将对数变换应用于幅度。这可以通过Librosa的librosa.power_to_db函数轻松完成。

# 将对数变换应用于幅度
scale_log_spectrogram = librosa.power_to_db(scale_spectrogram)

现在，再次绘制频谱图，这次使用对数幅度。

plot_spectrogram(scale_log_spectrogram, sr_scale, hop_size, y_axis=“linear”)

现在图像看起来好多了！我们可以看到音阶中每个音符的能量爆发。较低的部分是基频，较高的部分是谐波成分。

转换为对数频率频谱图

除了幅度，我们对频率的感知也是对数式的。因此，一个更符合感知的表示是同时使用对数幅度和对数频率。

在我们的绘图函数中，只需将y_axis参数从“linear”改为“log”即可。

plot_spectrogram(scale_log_spectrogram, sr_scale, hop_size, y_axis=“log”)

在对数频率表示下，音阶的上升过程在视觉上更加清晰和均匀，反映了我们实际的听觉体验。

比较不同音乐风格的频谱图

现在，让我们将同样的流程应用于其他三个音乐片段，并比较它们的频谱图。

以下是处理每个音频文件的步骤：

  1. 计算短时傅里叶变换。
  2. 计算幅度平方频谱图。
  3. 转换为对数幅度频谱图。
  4. 使用对数频率轴进行绘制。
# 为德彪西片段处理
debussy_stft = librosa.stft(debussy_signal, n_fft=frame_size, hop_length=hop_size)
debussy_spectrogram = np.abs(debussy_stft) ** 2
debussy_log_spectrogram = librosa.power_to_db(debussy_spectrogram)

# 为红辣椒乐队片段处理
redhot_stft = librosa.stft(redhot_signal, n_fft=frame_size, hop_length=hop_size)
redhot_spectrogram = np.abs(redhot_stft) ** 2
redhot_log_spectrogram = librosa.power_to_db(redhot_spectrogram)

# 为艾灵顿公爵片段处理
duke_stft = librosa.stft(duke_signal, n_fft=frame_size, hop_length=hop_size)
duke_spectrogram = np.abs(duke_stft) ** 2
duke_log_spectrogram = librosa.power_to_db(duke_spectrogram)

现在，绘制所有三个对数幅度、对数频率的频谱图。

plot_spectrogram(debussy_log_spectrogram, sr_debussy, hop_size, y_axis=“log”)
plot_spectrogram(redhot_log_spectrogram, sr_redhot, hop_size, y_axis=“log”)
plot_spectrogram(duke_log_spectrogram, sr_duke, hop_size, y_axis=“log”)

通过观察，我们可以发现不同音乐风格在频谱图上的差异：

  • 古典音乐（德彪西）：能量分布变化丰富，弦乐部分在高频有平滑的能量增长和衰减。
  • 摇滚音乐（红辣椒乐队）：在低频有大量重复的能量模式（如鼓和贝斯），节奏感强。
  • 爵士乐（艾灵顿公爵）：结合了规律的模式和更流畅、即兴的能量变化。

总结

本节课中我们一起学习了如何使用Python和Librosa库从音频中提取和可视化频谱图。我们从加载音频和计算短时傅里叶变换开始，然后通过计算幅度平方得到频谱图。为了获得更符合人类感知的可视化效果，我们应用了对数变换到幅度（librosa.power_to_db）和频率轴（设置y_axis=“log”），最终生成了对数幅度-对数频率频谱图。通过比较不同音乐风格的频谱图，我们看到这种表示方法能够揭示音频信号在时频域上的特征差异。在下一课中，我们将介绍另一种更贴近听觉感知的频谱图变体——梅尔频谱图。

017：梅尔频谱图详解 🎵

在本节课中，我们将要学习梅尔频谱图。这是一种在人工智能音频研究和实际应用中广泛使用的频谱图变体。

在上一节中，我们介绍了短时傅里叶变换以及如何提取和可视化频谱图。本节中，我们来看看梅尔频谱图，它解决了普通频谱图在频率感知上的一个关键问题。

从线性频率到感知频率

首先，我们回顾一下普通频谱图。在频谱图中，X轴代表时间，Y轴代表频率，每个点的颜色表示特定频率在特定时间点的能量强度。

这里的关键在于，普通频谱图的频率表示是线性的，单位是赫兹。这与人类感知音高的方式存在矛盾。人类对音高的感知是非线性的，我们在低频区的分辨率远高于高频区。

为了证明这一点，我们可以进行一个简单的听觉实验。比较两组音符：第一组是从C2到C4，第二组是从G6到A6。虽然它们在赫兹尺度上的频率差都约为200Hz，但我们的听觉会明显感觉到第二组音符的音高距离比第一组要近得多。这证实了我们的音高感知是非线性的。

理想的音频特征

现在，让我们思考一下，对于机器学习算法而言，理想的音频特征应该具备哪些特性。

以下是三个关键条件：

  1. 时频表示：能够展示信号中不同频率随时间的变化。
  2. 感知相关的幅度表示：幅度应以对数尺度表示，因为人耳对响度的感知也是对数式的。
  3. 感知相关的频率表示：频率表示应符合人类的听觉感知特性。

普通频谱图可以满足前两个条件（例如通过幅度谱），但无法满足第三个条件。而这正是梅尔频谱图的核心目标。

什么是梅尔尺度？🎼

“梅尔频谱图”由“频谱图”和“梅尔”两部分组成。“梅尔”这个概念与梅尔尺度相关。梅尔尺度是一种基于人耳听觉感知的音高标度。

梅尔尺度是对数性质的。这意味着，在梅尔尺度上相等的距离，对应着相同的感知音高距离。例如，500梅尔到510梅尔的感知差异，与1000梅尔到1010梅尔的感知差异是相同的。而这在赫兹线性尺度上是不成立的。

梅尔（Mel）一词是“旋律”的缩写，因为旋律主要由音高和节奏构成。

从赫兹到梅尔的转换公式如下：
m = 2595 * log10(1 + f/700)

反之，从梅尔转换回赫兹的公式为：
f = 700 * (10^(m/2595) - 1)

如何提取梅尔频谱图？🔧

提取梅尔频谱图可以概括为三个主要步骤。

以下是具体流程：

  1. 提取短时傅里叶变换：这与生成普通频谱图的第一步相同。
  2. 将幅度转换为分贝：即对幅度取对数，这也是普通频谱图的常见操作。
  3. 将频率转换为梅尔尺度：这是梅尔频谱图独有的、最关键的一步。

那么，如何将频率转换到梅尔尺度呢？这需要通过梅尔滤波器组来实现。

构建与应用梅尔滤波器组

频率转换过程包含以下几个子步骤。

以下是构建梅尔滤波器组的步骤：

  1. 选择梅尔带数量：这是一个超参数，通常在40到128之间选择，具体取决于任务。可以参考钢琴的88个键作为启发。
  2. 确定频率范围：设定STFT中要考虑的最低和最高频率（赫兹），并将其转换为梅尔值。
  3. 创建中心频率点：在梅尔尺度上，在最低和最高梅尔频率之间，生成等间距的N个点（N为梅尔带数量）。这些点就是各个梅尔带的中心频率。
  4. 将中心频率点转换回赫兹：使用逆转换公式，将梅尔中心频率点转换回赫兹。
  5. 创建三角滤波器：以每个中心频率点为中心，以前一个和后一个中心频率点为边界，构建一个三角滤波器。在中心点权重为1，在边界点权重为0。所有三角滤波器的集合就构成了梅尔滤波器组。

梅尔滤波器组在数学上可以表示为一个矩阵，其形状为 (num_mel_bands, n_fft/2 + 1)。

得到梅尔频谱图

最后一步是将梅尔滤波器组应用到频谱图上。

应用方法就是进行矩阵乘法：梅尔频谱图 = 梅尔滤波器组矩阵 × 频谱图矩阵。前提是滤波器组矩阵的列数等于频谱图矩阵的行数（即n_fft/2 + 1）。

相乘后得到的梅尔频谱图是一个新矩阵，其形状为 (num_mel_bands, num_frames)。在视觉上，它看起来和普通频谱图类似，但Y轴的单位从线性赫兹变成了感知相关的梅尔带。

为什么使用梅尔频谱图？💡

梅尔频谱图在AI音频和音乐信息检索领域被广泛使用。

以下是其主要应用场景：

  • 自动情绪识别
  • 音乐流派分类
  • 乐器分类
  • 以及其他众多音频分类任务

它之所以受欢迎，是因为它为机器学习模型提供了更符合人类听觉感知的输入特征，这通常能带来更好的模型性能。

总结

本节课中我们一起学习了梅尔频谱图。我们首先指出了普通频谱图在线性频率表示上的局限性。接着，我们探讨了梅尔尺度如何以对数方式模拟人耳的音高感知。然后，我们详细拆解了提取梅尔频谱图的步骤，核心在于构建和应用梅尔滤波器组将线性频率转换为梅尔带。最后，我们了解了梅尔频谱图在各类AI音频任务中的广泛应用。

下一节，我们将使用Python和LibROSA库来实际提取和可视化梅尔频谱图及梅尔滤波器组。

018：使用Python提取梅尔频谱图

在本节课中，我们将学习如何利用Python和Librosa库，将上一节介绍的梅尔频谱图理论付诸实践，从音频文件中提取梅尔频谱图。

导入必要库

首先，我们需要导入一些必要的Python库。Librosa用于音频处理，Librosa.display和Matplotlib用于可视化。

import librosa
import librosa.display
import matplotlib.pyplot as plt

加载音频文件

接下来，我们将加载一个音频文件。这个文件是一个C大调音阶的钢琴演奏，我们在之前提取普通频谱图的视频中也使用过它。

以下是加载和播放音频的步骤：

  1. 获取音频文件的路径。
  2. 使用Librosa加载音频，得到信号数组和采样率。
# 获取音频文件路径
audio_path = ‘path/to/your/audio/file.wav’

# 加载音频文件
signal, sr = librosa.load(audio_path, sr=None) # sr=None 保持原始采样率

理解梅尔滤波器组

在从音频中提取梅尔频谱图之前，我们先来回顾并可视化梅尔滤波器组。梅尔滤波器组是获取梅尔频谱图的关键，其原理是将滤波器组矩阵与普通频谱图进行矩阵乘法运算。

在Librosa中，我们可以使用一个实用函数来生成梅尔滤波器组，而无需手动实现所有步骤。

# 定义参数
frame_size = 2048
n_mels = 10

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/ee9105527e3fdfb9472f7eeda67c2563_3.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/ee9105527e3fdfb9472f7eeda67c2563_4.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/ee9105527e3fdfb9472f7eeda67c2563_5.png)

# 生成梅尔滤波器组
mel_filter_banks = librosa.filters.mel(sr=sr, n_fft=frame_size, n_mels=n_mels)

# 查看滤波器组的形状
print(mel_filter_banks.shape) # 输出应为 (n_mels, frame_size // 2 + 1)

滤波器组是一个矩阵，其第一维的大小等于梅尔带的数量（本例中为10），第二维的大小等于帧大小的一半加一（即 2048 // 2 + 1 = 1025），这对应于频谱图中的频率单元数。

可视化梅尔滤波器组

我们可以用两种方式可视化梅尔滤波器组。第一种是绘制每个三角滤波器的权重随频率变化的曲线。第二种是使用频谱图显示函数，将整个滤波器组矩阵显示为图像。

以下是使用频谱图方式可视化的代码：

# 将梅尔滤波器组可视化为图像
plt.figure(figsize=(10, 4))
librosa.display.specshow(mel_filter_banks, sr=sr, x_axis=‘linear’)
plt.colorbar(format=‘%+2.0f dB’)
plt.title(‘Mel Filter Banks’)
plt.xlabel(‘Frequency (Hz)’)
plt.ylabel(‘Mel Band’)
plt.tight_layout()
plt.show()

在生成的图像中，X轴表示频率（Hz），Y轴表示不同的梅尔带（本例中为10个）。颜色亮度代表在特定频率下，某个梅尔带的权重值，最亮的点（权重为1）对应每个梅尔带的中心频率。

提取梅尔频谱图

理解了滤波器组后，我们现在可以提取梅尔频谱图了。在Librosa中，这可以通过一个函数轻松完成，该函数在内部自动执行了提取普通频谱图、创建梅尔滤波器组并应用滤波器组这三个步骤。

以下是提取梅尔频谱图的代码：

# 提取梅尔频谱图
mel_spectrogram = librosa.feature.melspectrogram(y=signal, sr=sr, n_fft=frame_size, hop_length=512, n_mels=90)

# 查看梅尔频谱图的形状
print(mel_spectrogram.shape) # 输出应为 (n_mels, 时间帧数)

输出的形状中，第一维是梅尔带的数量（90），第二维是时间帧的数量（例如342），这取决于音频长度和跳幅大小。

转换为对数梅尔频谱图

由于人耳对振幅的感知是对数关系而非线性关系，因此将功率谱转换为分贝标度的对数谱是重要的一步。这不会改变频谱图的形状，但会转换矩阵中每个元素的值。

以下是转换代码：

# 转换为对数梅尔频谱图（分贝标度）
log_mel_spectrogram = librosa.power_to_db(mel_spectrogram, ref=np.max)

# 验证形状不变
print(log_mel_spectrogram.shape) # 形状保持不变

可视化梅尔频谱图

最后，我们可以将提取的梅尔频谱图可视化出来，观察音频信号在梅尔频率尺度下的时频特征。

以下是可视化代码：

plt.figure(figsize=(10, 4))
librosa.display.specshow(log_mel_spectrogram, sr=sr, hop_length=512, x_axis=‘time’, y_axis=‘mel’)
plt.colorbar(format=‘%+2.0f dB’)
plt.title(‘Log-Mel Spectrogram’)
plt.tight_layout()
plt.show()

在生成的频谱图中，X轴是时间，Y轴是梅尔带。对于我们的钢琴音阶，可以看到能量在梅尔带上依次上升的模式。为了更清晰地看到离散的梅尔带，我们可以将 n_mels 参数设置为一个较小的值（如10）并重新运行，这样在Y轴上就能看到10个清晰的色块。

总结

本节课中，我们一起学习了如何使用Python的Librosa库从音频文件中提取梅尔频谱图。我们回顾了梅尔滤波器组的概念，并演示了如何生成和可视化它。接着，我们使用 librosa.feature.melspectrogram 函数一步完成了梅尔频谱图的提取，并将其转换为对数尺度以便于分析。最后，我们可视化了结果，观察了音频在梅尔频率尺度下的特征。掌握这些步骤，你就具备了为机器学习任务准备音频特征的基础能力。

019：梅尔频率倒谱系数详解 🎵

在本节课中，我们将要学习一个非常重要的音频特征：梅尔频率倒谱系数。我们将从基础概念开始，逐步深入其数学原理、计算过程、优缺点以及应用场景，确保初学者能够轻松理解。

概述

梅尔频率倒谱系数是一种广泛用于语音和音乐处理的音频特征。它能够有效捕捉声音的“音色”信息，同时过滤掉我们不关心的细节（如音高）。本节课我们将详细拆解其构成、计算步骤和背后的原理。

梅尔频率倒谱系数是什么？

上一节我们介绍了梅尔频谱图，它是理解MFCC的重要基础。本节中，我们来看看MFCC到底是什么。

MFCC是“梅尔频率倒谱系数”的缩写。这个名字包含了几个关键部分：

  • 梅尔频率：指我们使用了梅尔刻度，这是一种与人耳对音高感知相关的非线性频率刻度。
  • 倒谱：这个词与“频谱”相关，它本质上是一个频谱的频谱。我们稍后会详细解释。
  • 系数：最终我们会得到一系列数值（系数），这些数值描述了声音片段的特征。

理解“倒谱” 🧠

“倒谱”是理解MFCC的核心。为了理解它，我们需要先了解一些历史背景和数学形式化。

历史背景与应用

倒谱的概念最初由MIT的研究人员在20世纪60年代提出，用于研究地震信号中的回声。随后，语音处理领域的研究者发现这个概念非常适用于语音分析。从70年代起，倒谱及其衍生特征（如MFCC）成为了语音识别、说话人识别等任务的首选特征，直到深度学习兴起。在21世纪初，MFCC也开始被广泛应用于音乐信息检索领域。

数学形式化

倒谱的计算过程可以通过一个数学公式来清晰地描述。以下是计算倒谱 C(t) 的步骤：

  1. 时域信号：我们从时域信号 x(t)（即波形）开始。
  2. 傅里叶变换：对 x(t) 应用离散傅里叶变换 F，得到频谱 X(f)，即从时域转换到频域。
    X(f) = \mathcal{F}\{x(t)\}
  3. 取对数：计算频谱 X(f) 的幅度（或功率）的对数，得到对数幅度谱 log(|X(f)|)。
  4. 逆傅里叶变换：对对数幅度谱应用逆傅里叶变换，得到的结果就是倒谱 C(t)。
    C(t) = \mathcal{F}^{-1}\{\log(|\mathcal{F}\{x(t)\}|)\}

关键理解：这个过程可以看作是在计算一个频谱的频谱。因为我们对频谱（对数幅度谱）进行了傅里叶分析（逆变换）。研究者们玩了一个文字游戏，将“spectrum”（频谱）的前四个字母“spec”倒过来，创造了“cepstrum”（倒谱）这个词。

可视化理解

让我们通过可视化来巩固这个概念：

  1. 从波形到频谱：我们从一个短时音频波形开始，通过傅里叶变换得到其功率谱（频率 vs. 功率）。
  2. 从频谱到对数谱：对功率谱的幅度值取对数，得到对数功率谱（频率 vs. 分贝）。这个信号看起来是连续的，并且由于原始信号中的谐波结构，它呈现出一定的周期性。
  3. 从对数谱到倒谱：现在，我们将这个对数功率谱视为一个时域信号。对这个“信号”应用逆傅里叶变换，就得到了倒谱。此时，横轴不再是频率，而是一种被称为倒频率的伪频率轴，单位通常是毫秒。

在倒谱图中，峰值（称为倒频峰）对应了原始对数谱中的周期性结构。第一个显著的倒频峰通常与原始信号的基频相关，因此倒谱也曾被用于基频检测。

倒谱在语音处理中的意义 🗣️

为了理解为什么倒谱如此有用，我们需要了解人类语音的产生机制。

语音产生模型

人类语音的产生可以简化为一个管道模型：

  1. 声门脉冲：由声带振动产生的一个类似噪声的、高频的信号，它主要携带音高（基频） 信息。
  2. 声道：包括舌头、牙齿、鼻腔、喉咙等。声道像一个滤波器，对声门脉冲进行滤波。
  3. 输出：滤波后的结果就是我们听到的语音信号。声道的形状决定了我们发出什么音素（元音或辅音），即声音的音色。

频谱分解与倒谱的作用

根据这个模型，语音信号 x(t) 可以看作是声道的脉冲响应 v(t)（即滤波器特性）与声门脉冲 g(t) 的卷积：

x(t) = v(t) * g(t)

转换到频域，卷积变为乘法：

X(f) = V(f) \cdot G(f)

然后，我们对两边取对数：

\log(|X(f)|) = \log(|V(f)|) + \log(|G(f)|)

这就是关键所在！ 取对数后，乘积变成了加法。log(|V(f)|) 对应声道频率响应（决定音色），log(|G(f)|) 对应声门脉冲谱（决定音高）。它们在倒谱域中天然地分离了：

  • 倒频率轴的低频部分（对应缓慢变化的包络）包含了 log(|V(f)|) 的信息，即共振峰和音色。
  • 倒频率轴的高频部分（对应快速变化的细节）包含了 log(|G(f)|) 的信息，即音高细节。

对于大多数语音识别任务，我们只关心音色（共振峰），不关心音高。因此，我们可以使用一个低通“倒频滤波器” 来滤除倒谱中高频部分（与声门脉冲相关）的系数，只保留低频部分（与声道响应相关）的系数。这样就得到了我们想要的、描述音色的特征。

从倒谱到梅尔频率倒谱系数 🎚️

理解了倒谱之后，MFCC就很容易理解了。它是在倒谱的基础上，引入了梅尔刻度，以更好地匹配人耳听觉特性。

MFCC的计算步骤

以下是计算MFCC的多步流程，其中许多步骤与计算倒谱共享：

  1. 预加重、分帧、加窗：对原始音频波形进行预处理。
  2. 傅里叶变换：对每一帧应用FFT，得到幅度谱。
  3. 计算功率谱：求幅度谱的平方。
  4. 应用梅尔滤波器组：将功率谱通过一组梅尔三角滤波器。这一步将线性频率刻度映射到梅尔刻度，并求每个滤波器内的能量和，得到梅尔频谱。这是与普通倒谱计算的第一处不同。
  5. 取对数：计算梅尔频谱的对数。
  6. 离散余弦变换：对对数梅尔频谱应用离散余弦变换。这是与普通倒谱计算的第二处不同（普通倒谱使用逆傅里叶变换）。DCT的结果就是梅尔频率倒谱系数。

每一步都具有一定的感知合理性：对数运算模拟了人耳对响度的非线性感知；梅尔刻度模拟了人耳对音高的非线性感知；DCT则帮助我们解相关并压缩数据。

为什么使用离散余弦变换？

使用DCT而非逆傅里叶变换有几个原因：

  • 实数系数：DCT输出的是实数系数，而FFT输出复数系数。对于MFCC，实数系数足够且更易于处理。
  • 解相关：DCT能有效解相关梅尔滤波器组能量（因为相邻滤波器有重叠），这有利于后续的机器学习算法。
  • 降维：DCT可以看作一种降维技术，它能够用少数几个系数捕捉对数梅尔频谱的主要形状。

应该取多少个系数？

传统上，我们只取前12-13个MFCC系数。这是因为：

  • 前几个系数包含了关于频谱包络（共振峰） 的大部分信息，这正是我们关心的音色特征。
  • 高阶系数包含更多关于频谱细节（音高） 的信息，这些信息对于语音识别通常不重要，甚至可能成为噪声。

为了提升模型性能，一个常见的做法是计算MFCC的一阶差分和二阶差分（也称为Delta和Delta-Delta系数）。它们描述了MFCC系数随时间的变化（类似于速度、加速度），对于捕捉动态特征非常有效。通常会将13个静态MFCC、13个一阶差分和13个二阶差分拼接起来，形成一个39维的特征向量。

MFCC的可视化

MFCC通常以热图形式可视化，非常类似于频谱图。横轴是时间帧，纵轴是MFCC系数索引（从0开始，0对应第一个系数）。每个点的颜色深浅代表该系数在对应时间帧上的数值大小。

MFCC的优缺点与应用总结 📊

优点

  • 聚焦音色：能有效描述频谱的大体结构（包络、共振峰），忽略细微的频谱细节和音高信息。
  • 感知相关：计算步骤融合了人耳听觉特性（对数、梅尔刻度）。
  • 久经考验：在语音和音乐处理领域有长期的成功应用历史。

缺点

  • 对加性噪声敏感：在嘈杂环境中性能会下降。
  • 包含人为假设：梅尔刻度、三角滤波器形状等选择带有一定的主观性和局限性，可能不是最优的，也限制了数据驱动方法从原始数据中自行学习特征的能力。
  • 不可逆：主要用于分析，难以完美地从MFCC重构回原始音频，因此不适用于合成任务。

主要应用

  • 语音处理：语音识别、说话人识别、语种识别。
  • 音乐处理：音乐流派分类、情绪分类、自动打标签、乐器识别等与音色相关的任务。

总结

本节课我们一起深入学习了梅尔频率倒谱系数。我们从其名称解析开始，深入探讨了核心概念“倒谱”的数学定义、历史背景及其在语音产生模型中的重要意义。随后，我们详细讲解了MFCC的计算步骤，理解了为何引入梅尔刻度和离散余弦变换，并讨论了系数选取、差分计算以及可视化方法。最后，我们总结了MFCC的优缺点和典型应用场景。

MFCC是连接传统信号处理与机器学习的重要桥梁，理解它对于掌握音频特征工程至关重要。在接下来的课程中，我们将使用Python实际提取并可视化MFCC，将理论付诸实践。

020：使用Python提取梅尔频率倒谱系数 🎵

在本节课中，我们将学习如何使用Python和Librosa库来提取音频信号的梅尔频率倒谱系数。我们将从加载音频文件开始，逐步完成MFCC的提取、可视化，并计算其一阶和二阶导数，最后将它们组合成一个综合性的音频特征。

概述

上一节我们从理论角度探讨了梅尔频率倒谱系数。本节我们将通过实践，学习如何使用Python代码来提取和操作MFCCs。

导入必要的库

首先，我们需要导入本教程中将要用到的Python库。

import librosa
import librosa.display
import IPython.display as ipd
import matplotlib.pyplot as plt
import numpy as np

加载音频文件

接下来，我们将加载一个音频文件。我们使用一段熟悉的古典音乐片段。

audio_file = ‘debusy.wav’

为了确认音频内容，我们可以在Jupyter笔记本中播放它。

ipd.Audio(audio_file)

现在，我们使用Librosa加载音频文件，获取其波形信号和采样率。

signal, sr = librosa.load(audio_file)

我们可以查看信号的形状，了解它包含的样本数量。

signal.shape

提取MFCCs

提取MFCCs是整个过程的核心。Librosa库提供了一个非常简便的函数来完成此操作。

我们将提取13个MFCC系数，这是一个常用的数量。

mfccs = librosa.feature.mfcc(y=signal, n_mfcc=13, sr=sr)

提取完成后，我们可以查看MFCCs矩阵的形状。

mfccs.shape

可视化MFCCs

为了直观理解MFCCs，我们可以将其可视化为频谱图。

以下是可视化步骤：

plt.figure(figsize=(25, 10))
librosa.display.specshow(mfccs, x_axis=‘time’)
plt.colorbar(format=‘%+2.0f dB’)
plt.show()

计算MFCCs的导数

MFCCs的一阶和二阶导数（Delta和Delta-Delta）描述了这些系数随时间的变化情况，是重要的补充特征。

以下是计算步骤：

# 计算一阶导数 (Delta MFCCs)
delta_mfccs = librosa.feature.delta(mfccs)

# 计算二阶导数 (Delta-Delta MFCCs)
delta2_mfccs = librosa.feature.delta(mfccs, order=2)

我们可以检查导数的形状，确认它们与原始MFCCs一致。

delta_mfccs.shape

同样地，我们可以将导数可视化。

以下是可视化Delta MFCCs的代码：

plt.figure(figsize=(25, 10))
librosa.display.specshow(delta_mfccs, x_axis=‘time’)
plt.colorbar(format=‘%+2.0f dB’)
plt.show()

以下是可视化Delta-Delta MFCCs的代码：

plt.figure(figsize=(25, 10))
librosa.display.specshow(delta2_mfccs, x_axis=‘time’)
plt.colorbar(format=‘%+2.0f dB’)
plt.show()

组合特征

最后，我们可以将原始MFCCs及其一阶、二阶导数在特征维度上拼接起来，形成一个包含静态和动态信息的综合性特征向量。

以下是组合步骤：

comprehensive_mfccs = np.concatenate((mfccs, delta_mfccs, delta2_mfccs))

查看组合后特征的形状，可以看到特征维度从13增加到了39。

comprehensive_mfccs.shape

总结

本节课中，我们一起学习了使用Python和Librosa库提取梅尔频率倒谱系数的完整流程。我们掌握了如何加载音频、提取MFCCs、将其可视化、计算并可视化其一阶和二阶导数，以及如何将这些特征组合成一个更全面的表示。现在，你应该能够将这些技术应用到自己的音频处理项目中。在下一节视频中，我们将探讨频域音频特征。

021：频域音频特征 🎵

在本节课中，我们将学习三种重要的频域音频特征：频带能量比、频谱质心和带宽。我们将探讨它们的数学原理、直观理解以及应用场景。

在之前的几节视频中，我们学习了梅尔频率倒谱系数。现在，我们将利用对傅里叶变换和短时傅里叶变换的深入理解，转向频域音频特征的分析。本节将重点介绍三种核心特征。

频带能量比 📊

频带能量比用于衡量信号中低频能量与高频能量的相对关系，可以理解为低频成分主导程度的指标。

以下是计算频带能量比的公式：

BER(t) = ( Σ_{n=1}^{F-1} |M_t(n)|^2 ) / ( Σ_{n=F}^{N} |M_t(n)|^2 )

其中：

  • M_t(n) 代表在时间帧 t 和频率仓 n 处的信号幅度。
  • N 是频谱图中的总频率仓数。
  • F 是分割频率，它将频谱划分为低频和高频两部分。

核心概念解析：

  • 分子：计算从第1个频率仓到分割频率 F-1 的所有频率仓的功率（幅度平方）之和，即低频带总能量。
  • 分母：计算从分割频率 F 到最高频率仓 N 的所有频率仓的功率之和，即高频带总能量。
  • 分割频率 F：这是一个任意设定的阈值频率。通常设为2000赫兹，但可根据需要调整。低于 F 的为低频，高于 F 的为高频。

这个公式会应用于频谱图的每一时间帧，从而得到一个随时间变化的频带能量比序列。

应用场景：频带能量比在音频处理中用途广泛，尤其常用于区分音乐与语音，以及在音乐流派分类、情绪分类等问题中作为关键特征。

频谱质心 ⚖️

频谱质心提供了幅度谱的“重心”位置，即能量最集中的频率带。它与一个重要的听觉感知特征——明亮度——有很好的对应关系。声音越明亮，频谱质心通常越高。

频谱质心在数学上定义为频率仓的加权平均值。

其计算公式如下：

SC(t) = ( Σ_{n=1}^{N} n * |M_t(n)| ) / ( Σ_{n=1}^{N} |M_t(n)| )

核心概念解析：

  • 权重：每个频率仓 n 的权重是该仓在时间帧 t 的幅度值 |M_t(n)|。
  • 加权平均：公式计算的是所有频率仓编号 n 以其幅度为权重的加权平均。结果 SC(t) 就是一个代表能量重心的频率值。

应用场景：频谱质心是音频和音乐分类任务中最核心、最常用的频域特征之一，在各种应用中都有广泛的历史。

带宽 📏

带宽与频谱质心密切相关。我们可以将带宽理解为围绕频谱质心的、有意义的频谱范围，或者看作是能量分布相对于质心的方差。它同样与听觉感知特征（如音色）直接相关。

带宽在数学上定义为各频率仓到频谱质心距离的加权平均值。

其计算公式如下：

BW(t) = ( Σ_{n=1}^{N} |M_t(n)| * |n - SC(t)| ) / ( Σ_{n=1}^{N} |M_t(n)| )

核心概念解析：

  • 距离：|n - SC(t)| 计算了每个频率仓 n 与当前帧的频谱质心 SC(t) 的绝对距离。
  • 加权平均：同样以幅度 |M_t(n)| 为权重，对所有频率仓的这个距离值求加权平均。
  • 物理意义：如果能量广泛分布在多个频率带上，带宽值会增大；如果能量集中在质心附近的少数频率带上，带宽值会减小。因此，带宽有时也被称为频谱扩展。

应用场景：与频谱质心类似，带宽也广泛应用于音乐流派分类、情绪分类等音乐信息检索任务。

总结与展望 🚀

本节课我们一起学习了三种基础的频域音频特征：频带能量比、频谱质心和带宽。它们都源于对短时傅里叶变换得到的频谱图进行计算，并在传统的机器学习时代（依赖于特征工程）被大量使用。虽然深度学习时代更倾向于使用频谱图或波形本身作为输入，但理解这些基础特征对于掌握音频分析的核心思想至关重要。

掌握了这些基本概念后，你将能更容易地理解其他更复杂的频域特征。

下节预告：本节我们探讨了这些特征的理论基础。下一节，我们将使用Python从零开始实现频带能量比特征，并将其应用于不同流派的音乐片段，可视化并观察是否能仅凭此特征区分它们。

022：Python实现频带能量比

在本节课中，我们将学习如何从零开始，使用Python实现一个重要的频域音频特征——频带能量比。我们将通过分析两段不同风格的音乐（古典乐和摇滚乐），来直观地理解这个特征如何反映音频信号的频谱能量分布。

上一节我们介绍了几种频域音频特征的理论基础。本节中，我们来看看如何将其中一个特征——频带能量比——通过代码实现出来。

准备工作

首先，我们需要导入必要的库并加载音频文件。我们将使用两段音乐：一段是德彪西的古典乐，另一段是红辣椒乐队的摇滚乐。

import librosa
import numpy as np
import matplotlib.pyplot as plt

# 加载音频文件
beethoven_path = ‘path_to_beethoven.wav’
rhcp_path = ‘path_to_red_hot_chili_peppers.wav’

beethoven_waveform, sr = librosa.load(beethoven_path, sr=22050)
rhcp_waveform, _ = librosa.load(rhcp_path, sr=sr)

# 定义帧长和跳数
frame_size = 2048
hop_size = 512

# 计算频谱图
beethoven_spec = librosa.stft(beethoven_waveform, n_fft=frame_size, hop_length=hop_size)
rhcp_spec = librosa.stft(rhcp_waveform, n_fft=frame_size, hop_length=hop_size)

计算分割频率对应的频点索引

频带能量比需要一个分割频率，将频谱分为低频和高频两部分。由于频谱图是离散的，我们需要一个函数将连续的频率值映射到最接近的离散频点索引上。

以下是实现此功能的步骤：

  1. 计算频谱图覆盖的总频率范围（0 Hz 到奈奎斯特频率）。
  2. 计算相邻频点之间的频率间隔。
  3. 将给定的分割频率除以这个间隔，得到对应的浮点数索引。
  4. 使用向下取整函数，将其转换为整数索引。

def calculate_split_frequency_bin(spectrogram, split_frequency, sample_rate):
    """
    将连续的分割频率映射到频谱图最接近的离散频点索引。
    参数:
        spectrogram: 输入的频谱图（复数矩阵）
        split_frequency: 分割频率（Hz）
        sample_rate: 音频采样率
    返回:
        分割频率对应的频点索引（整数）
    """
    frequency_range = sample_rate / 2
    frequency_delta_per_bin = frequency_range / spectrogram.shape[0]
    split_frequency_bin = np.floor(split_frequency / frequency_delta_per_bin)
    return int(split_frequency_bin)

计算频带能量比

有了分割频点索引，我们就可以计算每一帧的频带能量比了。其定义为低频带能量之和与高频带能量之和的比值。

以下是计算频带能量比的核心步骤：

  1. 获取分割频点索引。
  2. 将复数频谱图转换为功率谱（取绝对值后平方）。
  3. 转置功率谱，使第一维是时间（帧），第二维是频率，便于按帧迭代。
  4. 遍历每一帧，分别计算低频部分和高频部分的能量和。
  5. 将低频能量和除以高频能量和，得到该帧的频带能量比值。

def calculate_band_energy_ratio(spectrogram, split_frequency, sample_rate):
    """
    计算频谱图每一帧的频带能量比。
    参数:
        spectrogram: 输入的频谱图（复数矩阵）
        split_frequency: 分割频率（Hz）
        sample_rate: 音频采样率
    返回:
        一个NumPy数组，包含每一帧的频带能量比值
    """
    split_frequency_bin = calculate_split_frequency_bin(spectrogram, split_frequency, sample_rate)
    
    # 转换为功率谱并转置，使时间维度在前
    power_spectrogram = np.abs(spectrogram) ** 2
    power_spectrogram = power_spectrogram.T
    
    band_energy_ratio = []
    
    # 遍历每一帧
    for frequencies_in_frame in power_spectrogram:
        sum_power_low_frequencies = np.sum(frequencies_in_frame[:split_frequency_bin])
        sum_power_high_frequencies = np.sum(frequencies_in_frame[split_frequency_bin:])
        
        ber_current_frame = sum_power_low_frequencies / sum_power_high_frequencies
        band_energy_ratio.append(ber_current_frame)
    
    return np.array(band_energy_ratio)

应用与可视化

现在，我们可以对两段音乐应用这个函数，并可视化结果以进行比较。

# 设置分割频率为2000 Hz
split_freq = 2000

# 计算两段音乐的频带能量比
ber_beethoven = calculate_band_energy_ratio(beethoven_spec, split_freq, sr)
ber_rhcp = calculate_band_energy_ratio(rhcp_spec, split_freq, sr)

# 将帧索引转换为时间（秒）
frames = range(len(ber_beethoven))
times = librosa.frames_to_time(frames, hop_length=hop_size)

# 绘制结果
plt.figure(figsize=(25, 10))
plt.plot(times, ber_beethoven, color=‘blue’, label=‘Beethoven (Classical)’)
plt.plot(times, ber_rhcp, color=‘red’, label=‘RHCP (Rock)’)
plt.xlabel(‘Time (s)’)
plt.ylabel(‘Band Energy Ratio’)
plt.title(‘Band Energy Ratio Comparison: Classical vs. Rock’)
plt.legend()
plt.show()

结果分析

从生成的图表中，我们可以清晰地看到两条曲线的显著差异：

  • 蓝色曲线（古典乐）的频带能量比值普遍远高于红色曲线（摇滚乐）。
  • 这表明在古典音乐片段中，大部分声学能量集中在频谱的低频部分。
  • 而在摇滚乐片段中，能量在频谱上的分布更为均衡，高频部分（如鼓、镲片等乐器产生的噪声）贡献了更多能量。

这种差异正是频带能量比特征能够有效区分不同音乐风格或音频事件的关键所在。

总结

本节课中我们一起学习了如何从零开始实现频带能量比。我们首先编写了将连续频率映射到离散频点的函数，然后实现了计算每一帧能量比的核心逻辑，最后通过可视化对比了该特征在古典乐和摇滚乐中的不同表现。频带能量比是一个简单但强大的特征，能够有效表征音频信号频谱能量的分布情况，在音乐信息检索和音频分类等任务中非常有用。下一节课，我们将学习如何使用librosa库来便捷地计算其他频域音频特征。

023：使用Python和Librosa提取频谱质心与带宽 🎵

在本节课中，我们将学习如何使用Librosa库来计算音频信号的频谱质心与频谱带宽。这是本系列的最后一课，我们将通过分析三首不同类型的音乐片段来实践这两个重要的频谱特征。

概述与准备工作

首先，我们需要导入必要的库并加载音频文件。我们将使用三首风格迥异的音乐：一首德彪西的管弦乐作品、一首红辣椒乐队的歌曲以及一首艾灵顿公爵的爵士乐。

import librosa
import librosa.display
import matplotlib.pyplot as plt
import numpy as np

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/29e74a766bda69812908c49e36510eae_4.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/29e74a766bda69812908c49e36510eae_5.png)

# 加载音频文件
debussy_path = ‘audio/debussy.wav’
redhot_path = ‘audio/redhot.wav’
duke_path = ‘audio/duke.wav’

debussy, sr_debussy = librosa.load(debussy_path)
redhot, sr_redhot = librosa.load(redhot_path)
duke, sr_duke = librosa.load(duke_path)

接下来，我们需要设置用于特征提取的帧长和跳跃长度参数。

FRAME_SIZE = 2048
HOP_LENGTH = 512

提取频谱质心

上一节我们完成了准备工作，本节中我们来看看如何提取频谱质心。频谱质心可以被视为频谱的“重心”或平均频率，它反映了声音的明亮度。

以下是使用Librosa提取频谱质心的步骤：

# 为德彪西作品计算频谱质心
spectral_centroid_debussy = librosa.feature.spectral_centroid(y=debussy,
                                                              sr=sr_debussy,
                                                              n_fft=FRAME_SIZE,
                                                              hop_length=HOP_LENGTH)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/29e74a766bda69812908c49e36510eae_12.png)

# 为红辣椒乐队歌曲计算频谱质心
spectral_centroid_redhot = librosa.feature.spectral_centroid(y=redhot,
                                                             sr=sr_redhot,
                                                             n_fft=FRAME_SIZE,
                                                             hop_length=HOP_LENGTH)

# 为艾灵顿公爵爵士乐计算频谱质心
spectral_centroid_duke = librosa.feature.spectral_centroid(y=duke,
                                                           sr=sr_duke,
                                                           n_fft=FRAME_SIZE,
                                                           hop_length=HOP_LENGTH)

Librosa返回的是一个二维数组，其中第一维是1，第二维是帧的数量。我们通常只需要时间序列上的值。

# 获取一维的时间序列数据
sc_debussy = spectral_centroid_debussy[0]
sc_redhot = spectral_centroid_redhot[0]
sc_duke = spectral_centroid_duke[0]

可视化频谱质心

计算完成后，我们可以将三首音乐的频谱质心随时间变化的情况绘制在同一张图中进行比较。

以下是创建可视化图表的代码：

# 计算时间轴
frames = range(len(sc_debussy))
t = librosa.frames_to_time(frames, hop_length=HOP_LENGTH)

# 绘制图形
plt.figure(figsize=(25, 10))

plt.plot(t, sc_debussy, color=‘b’, label=‘Debussy (Orchestral)’)
plt.plot(t, sc_redhot, color=‘r’, label=‘Red Hot Chili Peppers (Rock)’)
plt.plot(t, sc_duke, color=‘y’, label=‘Duke Ellington (Jazz)’)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/29e74a766bda69812908c49e36510eae_20.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/29e74a766bda69812908c49e36510eae_21.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/29e74a766bda69812908c49e36510eae_23.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/29e74a766bda69812908c49e36510eae_25.png)

plt.xlabel(‘Time (seconds)’)
plt.ylabel(‘Spectral Centroid (Hz)’)
plt.legend()
plt.show()

观察图表可以发现，红辣椒乐队摇滚乐的频谱质心整体上高于德彪西的管弦乐和艾灵顿公爵的爵士乐。这通常是因为摇滚乐或电子音乐中高频成分更丰富。

提取频谱带宽

了解了频谱质心后，我们进一步看看频谱带宽。频谱带宽衡量的是频谱能量围绕质心的分散程度。

提取频谱带宽的代码与提取质心非常相似：

# 为三首音乐计算频谱带宽
spectral_bandwidth_debussy = librosa.feature.spectral_bandwidth(y=debussy,
                                                                sr=sr_debussy,
                                                                n_fft=FRAME_SIZE,
                                                                hop_length=HOP_LENGTH)
spectral_bandwidth_redhot = librosa.feature.spectral_bandwidth(y=redhot,
                                                               sr=sr_redhot,
                                                               n_fft=FRAME_SIZE,
                                                               hop_length=HOP_LENGTH)
spectral_bandwidth_duke = librosa.feature.spectral_bandwidth(y=duke,
                                                             sr=sr_duke,
                                                             n_fft=FRAME_SIZE,
                                                             hop_length=HOP_LENGTH)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/29e74a766bda69812908c49e36510eae_33.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/29e74a766bda69812908c49e36510eae_34.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/29e74a766bda69812908c49e36510eae_36.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/29e74a766bda69812908c49e36510eae_38.png)

# 获取一维时间序列数据
bw_debussy = spectral_bandwidth_debussy[0]
bw_redhot = spectral_bandwidth_redhot[0]
bw_duke = spectral_bandwidth_duke[0]

可视化频谱带宽

同样地，我们可以将频谱带宽的可视化结果绘制出来。

以下是绘制频谱带宽图表的代码：

# 使用相同的时间轴（因为帧数相同）
plt.figure(figsize=(25, 10))

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/29e74a766bda69812908c49e36510eae_46.png)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/29e74a766bda69812908c49e36510eae_47.png)

plt.plot(t, bw_debussy, color=‘b’, label=‘Debussy (Orchestral)’)
plt.plot(t, bw_redhot, color=‘r’, label=‘Red Hot Chili Peppers (Rock)’)
plt.plot(t, bw_duke, color=‘y’, label=‘Duke Ellington (Jazz)’)

![](https://github.com/OpenDocCN/dsai-notes-pt1-zh/raw/master/docs/valerio-ml-aud-sgnproc/img/29e74a766bda69812908c49e36510eae_49.png)

plt.xlabel(‘Time (seconds)’)
plt.ylabel(‘Spectral Bandwidth (Hz)’)
plt.legend()
plt.show()

从图中可以观察到，与频谱质心的趋势类似，古典乐和爵士乐的频谱带宽通常小于摇滚乐或使用大量电子乐器、打击乐的音乐。这是因为前者的频谱能量通常更集中。

总结

本节课中我们一起学习了如何使用Librosa库提取和可视化两个关键的频谱特征：频谱质心和频谱带宽。

我们通过分析三首不同风格的音乐，实践了完整的计算与绘图流程。频谱质心反映了声音的明亮度，而频谱带宽则描述了频谱能量的分散程度。掌握这些特征的提取方法，是进行音频机器学习任务（如音乐分类、情感分析或声音事件检测）的重要基础。

本系列课程至此结束，希望大家已经建立了坚实的音频信号处理知识基础，能够在未来的AI音频或音乐工程项目中灵活运用。

posted @ 2026-03-26 12:20  布客飞龙III  阅读(7)  评论(0)    收藏  举报