音频分析和处理简介

在这个系列的第一部分，我们看了傅里叶变换、短时傅里叶变换、梅尔刻度、滤波器组、梅尔频率表系数的理论和实现，以及光谱特征，如光谱中心点、光谱带宽和光谱对比。

到目前为止，我们还没有尝试去理解任何一种特定的声音。我们对声音分析和任何声音在频域中的样子采取了一种普遍的方法。在这篇文章中，我们将特别研究一下音乐。

把这个项目带入生活

在梯度上运行

简介

我们在上一篇文章中看到了声音是如何在空间中传播的。我们把这些波称为纵波，并把它们的波动称为压缩和稀疏因素。

以吉他为例。在吉他中，这些最终的声波是由不同的吉他弦的振动产生的。所有的物体都有一个自然频率，当敲击、弹奏或受到其他干扰时，它们会振动。

当一个物体被迫以与自然频率共鸣的方式振动时，介质往往会产生我们所说的驻波。驻波是具有恒定振幅轮廓的波。

驻波

峰值在X轴上总是相同的。零振幅的点也是如此。振幅为零的点根本不显示任何运动，被称为结点。

基准频率和谐波

在一个介质中产生驻波的频率被称为谐波。在谐波频率以外的频率，振动是不规则和不重复的。

任何特定乐器产生的最低频率被称为基本频率或乐器的第一次谐波。基频的振动可以表示为一个有两个节点的驻波，一个在琴弦的两端，一个反节点在琴弦的中心。

第一次谐波

第二个谐波的频率是第一个谐波的两倍，而第三个谐波的频率是第一个谐波的三倍。以此类推，不一而足。

谐波

在这里可以找到一个很好的可视化动画，将前四个谐波作为驻波和它们的纵向对应物。

音符、八度和频率

某些频率在一起播放时，我们的耳朵听起来很舒服。这样的频率被称为 "合音"。

在音乐中，音程是衡量不同声音的音高差异的一种方式。音乐中流行的音程可以用以下的整数比来表示。

1. 齐声 --> 1:1
2. 八度 --> 2:1
3. 大六度 --> 5:3
4. 完全五度 --> 3:2
5. 完全四度 --> 4:3
6. 大三度 --> 5:4
7. 小三度 --> 6:5

所以你可以说880Hz是比440Hz高一个八度。或者，如果440Hz是第一个谐波，那么880Hz是第二个谐波。

将两个声音放在一起播放可以创造出有趣的现象，如节拍。当两个频率稍有不同的声音同时播放时，声音的振幅会有一个周期性的变化，其频率是原来两个频率之间的差异。

$f\_{beat} = f\_{1}- f\_{2}$

在音乐中，频率被划分为不同的音高等级。在英语地区，它们由字母表的前7个字母来表示。有些人可能还知道这是Do、Re、Mi、Fa、Sol、La和Ti，而其他人可能知道这是Sa、Re、Ga、Ma、Pa、Dha和Ni。八分音符或八度音符的名称与第一个音符相同，在这之后，这个模式又会重复出现。

今天，大多数乐器都使用12音半音阶，它将一个八度分成12个对数音阶上的等距部分。每个部分据说相隔一个半音，这是音乐中最小的音程。大多数乐器都被调到了等分音阶上，中间八度的A音被调到了440Hz。

建立在C上的12音半音阶在每个八度都有以下12个音。C, C♯/D♭, D, D♯/E♭, E, F, F♯/G♭, G, G♯/A♭, A, A♯/B♭和B。

升读为尖音，降读为平音。尖音意味着高一个半音，平音意味着低一个半音。

半音阶

为了找到中间八度的频率，我们需要找到中间C的频率。中间C或C4的频率是C5的一半。

一旦我们有了中间C，我们可以创建一个大小为13的数组，在C4的对数和C5的对数之间等距排列。我们还将从数组中删除最后一个元素（或本例中的C5）。

代码如下所示。

import numpy as np

C5 = 2**np.linspace(np.log2(440), np.log2(880), 13)[3]
freqs = 2**np.linspace(np.log2(C5/2), np.log2(C5), 13)[:-1]

产生的频率。

array([261.6255653 , 277.18263098, 293.66476792, 311.12698372,
       329.62755691, 349.22823143, 369.99442271, 391.99543598,
       415.30469758, 440.        , 466.16376152, 493.88330126])

色度表示法

一个音调可以分成两个部分：音高和色度。

音高是指八度数。例如，中C也被称为C4。

色度特征，也被称为色度图，为音乐的可视化提供了一种强大的方式。色度特征与我们刚才看到的12音阶中的七个音级或12个音符密切相关。色度表征将一个声音在任何特定时间存在的所有频率分成12个桶中的一个，每个桶代表一个音符。

正如我们所讨论的，比A音高一个八度的音符仍然是A音，尽管其频率是前一个A音的两倍。枚举从0到11的色度特征，C音符的索引为0，下面的索引由比前一个索引高一个半音来填充。

短时间傅里叶变换可以通过得到STFT值与色度过滤器的点积，并对输出进行归一化来转换为色度表示。

librosa 滤波器可以用下面的片段进行可视化。

chroma_filters = librosa.filters.chroma(sampling_rate, 2048)
fig, ax = plt.subplots(nrows=6, ncols=2, figsize=(12,9))

for i, row in enumerate(ax):
    for j, col in enumerate(row):
        col.plot(chroma_filters[i + j])
        
plt.show()

这些图看起来就像下面的图片。

色度滤波器

我们将使用与本系列第一部分中相同的音频样本。

example_name = 'nutcracker' 
audio_path = librosa.ex(example_name)

x, sampling_rate = librosa.load(audio_path, sr=None)

你可以用librosa 库得到色度表示，使用下面的代码片段。

S = librosa.stft(x)**2
chroma = librosa.feature.chroma_stft(S=S, sr=sampling_rate)
fig, ax = plt.subplots(figsize=(15,9))
img = librosa.display.specshow(chroma, y_axis='chroma', 
                               x_axis='time', ax=ax)
fig.colorbar(img, ax=ax)
ax.set(title='Chromagram')

色度图

爆发点检测

Onset detection指的是一组方法，让我们通过声音定位音符的起始点。有几种方法可以做到这一点。它们大致可以分为以下几类。

• 基于能量的方法
• 基于音高的方法
• 基于相位的方法
• 监督下的学习

基于能量的方法

在基于能量的方法中，频谱被用来测量时-频域的能量变化。过去，频谱值的一阶差异被用来寻找起始点，但并不精确。根据心理声学原理，在较低的振幅下，类似的频率变化可以被更好地检测出来。一个相对差分的方案可以帮助我们更好地找到峰值。

例如，频谱D可以定义如下。

$D\_{m}(n) = 20log\_{10}(|X\_m(n)|^{2}) - 20log\_{10}(|X\_m(n-1)|^{2})$

其中$X\_{m}(n)$是输入信号的离散STFT。

常用的基于能量的检测方法可以归纳为以下几点。

$$ M=\frac{1}{N}H(D_{m}(n))$$sum_{m=1}^{N} $

其中$H(x) = 夫拉克{x + |x|}{2}$是半波整流函数，N是频谱D中的总频段数，M是检测函数。

探测函数被进一步用移动平均法平滑。一个简单的选峰操作被应用于寻找起始点。数值超过一定阈值的峰被返回作为起始点。

基于相位的方法

STFT可被视为一个带通滤波器，其中$X\_m(n)$表示$m^{th}$滤波器的输出。在只有一个稳定的正弦分量通过$m^{th}$带通滤波器的情况下，$m^{th}$滤波器的输出必须有一个恒定或接近恒定的频率。

因此，连续解包的相位值之间的差异也必须保持近似恒定。

$phi\_{m}(n)- phi\_{m}(n - 1) ´approx phi\_{m}(n - 1) - phi\_{m}(n - 2)$

或

$Bigtriangleup phi\_{m}(n) \`approx \`phi\_{m}(n) - 2 \`phi\_{m}(n - 1) - \`phi\_{m}(n - 2) ´approx 0$

在转换过程中，频率并不恒定。因此，$bigtriangle up phi\_{m}(n)$ 趋向于高。

基于音高的方法

基于音高的方法利用音高特性将声音分成稳态和瞬态两部分，这样就可以只在瞬态部分找到起始点。

在《随意的节奏分组》中，Kristoffer Jensen提出了一个检测功能，称为感知频谱通量，其中频段的差异是由等响度轮廓线来权衡的。

等响度等值线或弗莱彻-蒙森曲线是听者感知稳定音调的恒定振幅时，在频域中的声音轮廓线。

$PFS\_{n} = sum\_{m=1}^{N}W(X\_{m}(n)^{1/3}- X\_{m}(n - 1)^{1/3})$

$X\_{m}$是STFT的幅值，使用汉宁窗获得。$W$是频率加权，用于获得更接近人类响度轮廓的数值，功率函数用于模拟强度-响度幂律。功率函数进一步减少了随机振幅变化。

监督学习

神经网络也被用来检测音频信号中的起始帧。一个典型的基于神经网络的发病检测管道看起来像这样。

来源：eecs.qmul.ac.uk/~josh/docum…

一个音频信号首先被转换为频谱图。这个频谱图被传递给一个神经网络，该网络试图将每一帧分类为发病或不发病。被归类为起始点的，再通过像阈值处理这样的挑峰操作。

有几种神经网络结构可以用于这个问题。最流行的选择是递归网络和卷积神经网络。最近，音频变换器正在成为音乐分类的一个流行选择。

同样的网络也可以被重新利用来进行发声检测。例如，Miyazaki等人使用卷积网络和变压器块的混合体进行发病检测。

librosa ，他们的API中有发病检测功能，可以按以下方式使用。

o_env = librosa.onset.onset_strength(x, sr=sampling_rate)
times = librosa.times_like(o_env, sr=sampling_rate)
onset_frames = librosa.onset.onset_detect(onset_envelope=o_env, sr=sampling_rate)

fig, ax = plt.subplots(nrows=2, sharex=True, figsize=(15,9))
librosa.display.specshow(S_dB, x_axis='time', 
                         y_axis='log', ax=ax[0])

ax[0].set(title='Power spectrogram')
ax[0].label_outer()
ax[1].plot(times, o_env, label='Onset strength')
ax[1].vlines(times[onset_frames], 0, o_env.max(), 
             color='r', alpha=0.9,
             linestyle='--', label='Onsets')
ax[1].legend()

plt.savefig('onsets.png')

可视化将能够指出通过音频播放新音符的地方的估计位置。

起始点检测

前五秒的可视化可以像这样得到。

o_env = librosa.onset.onset_strength(x, sr=sampling_rate)
total_time = librosa.get_duration(x)
five_sec_mark = int(o_env.shape[0]/total_time*5)
o_env = o_env[:five_sec_mark]
times = librosa.times_like(o_env, sr=sampling_rate)
onset_frames = librosa.onset.onset_detect(onset_envelope=o_env, sr=sampling_rate)

fig, ax = plt.subplots(nrows=2, sharex=True, figsize=(15,9))
librosa.display.specshow(S_dB[:, :five_sec_mark], x_axis='time', 
                         y_axis='log', ax=ax[0])

ax[0].set(title='Power spectrogram')
ax[0].label_outer()
ax[1].plot(times, o_env, label='Onset strength')
ax[1].vlines(times[onset_frames], 0, o_env.max(), 
             color='r', alpha=0.9,
             linestyle='--', label='Onsets')
ax[1].legend()

plt.savefig('onsets-5secs.png')

开始检测（前五秒

正如我们在上面看到的，垂直线在发现音符开始的地方正确地排成一排。竖线与幅度曲线作图。频谱图在上面显示了相同的X轴，以便直观地比较音符的起始点。

节拍和节奏

音乐理论中的节拍是一个时间单位，追踪音乐中的周期性元素。你可以把节拍看作是你在听一首歌时拍打你的脚的次数。一首歌曲的节奏是每分钟的拍子数。节奏可能在歌曲的过程中发生变化。

检测节拍和估计歌曲的节奏与检测发病率密切相关。通常情况下，检测歌曲的节拍和节奏分三步进行。

• 发声检测：根据频谱通量计算发声的位置
• 周期性估计：计算起音位置的周期性模式
• 节拍位置估计：根据起始点和周期性信息计算节拍的位置。

在《音乐信号的节奏和节拍估计》中，作者给出了频谱通量的以下定义。

$E(m,n)= sum\_{n} h(n - k)G(m,n)$

其中$h(n)$近似于一个微分器滤波器，并且

$$ G(m,n)= 矩阵cal{F}\{|X(m,n)|}$

其中$X(m,n)$是我们信号的STFT。

变换$mathcal{F}$将信号的STFT通过低通滤波器和使用反双曲正弦函数的非线性压缩算法。

然后，作者使用一个中值滤波器来挑出真正的音符攻击而不是低振幅峰值的起始点。

这个检测功能的结果可以被认为是一个在音符攻击时具有高幅度的脉冲序列。这可以通过两种方式实现。

频谱乘积

这里的想法是，强谐波是在基本频率的整数倍的频率上发现的。

为了找到这个频率，我们计算信号的离散傅里叶变换，得到其所有数值的乘积，从而得到一个强化的频率。

然后，通过挑选出与强化频率值的最大峰值相对应的频率指数，就可以很容易地得到估计的节奏T。节奏被要求在60到200BPM之间，即每分钟的节拍。

自相关

寻找音频信号中的周期性的经典方法是使用自相关。ACF或自相关函数可以定义为

$r(\\tau)=sum\_{k}p(k + \\tau)p(k)$

其中$tau$是滞后值。为了找到估计的节奏，对$r(\\tau)$的三个最大峰值的滞后进行分析。在没有找到倍数关系的情况下，最大峰值的滞后被当作节拍周期。

为了找到节拍位置，我们使用频谱乘积或自相关函数创建了一个节奏为T的人工脉冲序列。然后我们将其与原始检测功能的输出进行交叉关联。

我们把相关度最大的时间指数称为t_{0}$和第一拍。在一个给定的分析窗口中，通过增加一个节拍周期$mathcal{T}$找到第二个和随后的节拍

$t_{i} = t_{i-1} + \mathcal{T}$$

librosa ，我们可以计算出静态节奏、动态节奏，还可以用tempogram直观地看到BPM随时间的变化。

要计算静态节奏，运行以下代码。

onset_env = librosa.onset.onset_strength(x, sr=sampling_rate)
tempo = librosa.beat.tempo(onset_envelope=onset_env, sr=sampling_rate)
print(tempo)

你应该得到一个单元素数组作为输出。

array([107.66601562])

要得到动态节拍。

dtempo = librosa.beat.tempo(onset_envelope=onset_env, sr=sampling_rate,
                            aggregate=None)
print(dtempo)

结果是一个数组。

array([112.34714674, 112.34714674, 112.34714674, ..., 107.66601562, 107.66601562, 107.66601562])

你也可以用tempogram将动态节拍可视化。

fig, ax = plt.subplots(figsize=(15,9))
tg = librosa.feature.tempogram(onset_envelope=onset_env, sr=sampling_rate,
                               hop_length=hop_length)
librosa.display.specshow(tg, x_axis='time', y_axis='tempo', cmap='magma', ax=ax)
ax.plot(librosa.times_like(dtempo), dtempo,
         color='c', linewidth=1.5, label='Tempo estimate (default prior)')
ax.set(title='Dynamic tempo estimation')
ax.legend()
plt.savefig('tempogram.png')

动态节拍估计

谱图分解

频谱图可以提供更多关于一个声音的信息，而不仅仅是关于起音和频率的视觉线索。我们可以使用频谱图来分离一段音乐的不同组成部分。

最常见的使用情况之一是将和声与打击性声音分开，例如钢琴、键盘、吉他和长笛与鼓、邦戈斯等相比。这可以用Harmonic Percussive Source Separation算法或HPSS来完成。

谐波音源分离（HPSS

来源：www.audiolabs-erlangen.de/content/05-…

HPSS算法有4个步骤。

1. 计算频谱图
2. 垂直和水平中值滤波
3. 对过滤后的频谱图进行二进制屏蔽
4. 反转遮蔽后的频谱图

我们知道，给定一个$N$数值的有序列表，如果$N$是奇数，中位数就是$((N-1)/2)^{th}$数值，否则就是$(N/2)^{th}$和$((N+1)/2)^{th}$项之和。

然后，给定一个矩阵$B在R\_{N×K}$中，我们定义谐波和冲击中值滤波器如下。

$medfilt\_{h}(B)(n,k) = median({B(m-l\_{h},k),...,B(m+l\_{h},k) })$

$medfilt\_{p}(B)(n,k) = median({B(n,k-l\_{p}),...,B(n,k+l\_{p})})$

其中$2l\_{h}+ 1$和2l_{p}+ 1$分别为谐波滤波器和冲击滤波器的长度。增强的频谱图被他们计算如下。

$phi\_{h} = medfilt\_{h}(|X(n, k )|^{2})$

$phi\_{p} = medfilt\_{p}(|X(n, k )|^{2})$

其中$X(n,k)$是声音信号的STFT。

滤波后的频谱图采用二进制掩码。对于谐波谱图，如果$phi\_{h}(n,k)$的值大于或等于$phi\_{p}(n,k)$，则元素$M\_{h}(n,k)$的值为1$。否则$M_{h}(n,k)$的值为0。

同样，对于冲击波频谱，如果$phi\_{p}(n,k)$的值大于$phi\_{h}(n,k)$，元素$M\_{p}(n,k)$的值为1$。否则$M_{p}(n,k)$的值为0。

请注意，这种表述意味着$M\_{h}(n,k)+M\_{p}(n,k)=1$，适用于所有$n$和$k$。为了得到打击和谐波元素的频谱图，我们在原始STFT $X$与二元掩码$M\_{p}$和$M\_{h}$之间做一个元素相乘。

由于我们对二元掩码的定义，每个时间点上的每个频仓都被分配给$X\_{h}$或$X\_{p}$。

$X\_{h} = X\\bigodot M\_{h}$

$X\_{p} = X \\bigodot M\_{p}$

一旦你有了这些，你就可以通过应用反STFT得到信号的打击性和谐波性元素。

librosa 为HPSS算法提供了一个开箱即用的实现。

S = librosa.stft(x)
H, P = librosa.decompose.hpss(S)

fig, ax = plt.subplots(nrows=3, sharex=True, sharey=True)
img = librosa.display.specshow(librosa.amplitude_to_db(np.abs(S), ref=np.max),
                         y_axis='log', x_axis='time', ax=ax[0])
ax[0].set(title='Full power spectrogram')
ax[0].label_outer()
librosa.display.specshow(librosa.amplitude_to_db(np.abs(H), ref=np.max(np.abs(S))),
                         y_axis='log', x_axis='time', ax=ax[1])
ax[1].set(title='Harmonic power spectrogram')
ax[1].label_outer()
librosa.display.specshow(librosa.amplitude_to_db(np.abs(P), ref=np.max(np.abs(S))),
                         y_axis='log', x_axis='time', ax=ax[2])
ax[2].set(title='Percussive power spectrogram')
fig.colorbar(img, ax=ax, format='%+2.0f dB')

HPSS（分离声音信号中的击弦和谐波部分

你可以通过以下方式得到分离的信号。

y_harm = librosa.istft(H)
y_perc = librosa.istft(P)

librosa 实现还可以将信号分解成三个元素，而不是两个：谐波、打击乐和残余。你可以通过将margin参数的值设置为大于1的任何值来做到这一点。

H, P = librosa.decompose.hpss(S, margin=3.0)
R = S - (H + P)
y_harm = librosa.istft(H)
y_perc = librosa.istft(P)
y_res = librosa.istft(R)

总结

在这篇文章中，我们利用HPSS算法探索了音符、谐波、八度、色度表示、起音检测方法、节拍、节奏、tempograms和谱图分解。

在本系列的下一部分，我们将尝试用深度学习建立一个文本到语音引擎。正如我们都喜欢说的那样--敬请期待吧

今天就为你的机器学习工作流程增加速度和简单性吧

开始联系销售