Python语音分析从绘制好看的波形图和语谱图开始

Posted 2023-03-30 文程公子

目录

1 Python-Librosa库简介

2 音频文件的加载

2.1 返回值与参数

2.2 音频加载示例

2.3 简析返回值y

3 波形图和语谱图的绘制

3.1 绘制波形图

3.1.1 waveshow()方法

3.1.2 波形图绘制示例

3.2 绘制语谱图

3.2.1 specshow()方法

3.2.2 线性频率的语谱图

3.2.3 对数频率的语谱图

3.2.4 梅尔频率的语谱图

4 语谱图颜色的设置（超级好看）

1 Python-Librosa库简介

Librosa是一个强大的用于处理声音信号的第三方库，其官网提供了比较详细的官方文档供使用者学习，主页链接贴在这里：librosa — librosa 0.8.1 documentation，下图是官方文档的主页：

官网提供了多种Librosa库的安装方法，包括pypi、conda和源码安装，最方便的当然是pip安装，而且这种安装方法确保Librosa库满足所有依赖。我们可以使用清华镜像来提升安装速度：

pip install -i https://pypi.tuna.tsinghua.edu.cn/simple librosa

目前官方文档对应的版本是0.8.1，Github（GitHub - librosa/librosa: Python library for audio and music analysis）中也是该版本，应该是当前最高版本，不久的将来会出现0.9.0版本。

2 音频文件的加载

2.1 返回值与参数

将音频文件加载到程序中直接调用 librosa.load()方法即可实现，方法中的参数如下：

librosa.load(path, sr=22050, mono=True, offset=0.0, duration=None, 
dtype=<class 'numpy.float32'>, res_type='kaiser_best')

先交代一下该方法的返回值。返回值有两个：y和sr。y是音频的时间序列，是一个float32的numpy数组。sr是采样率(sampling rate)。

path 是文件路径。可以传入各种格式的文件，对于语音分析来说，最好直接传入.wav格式的音频文件，wav就是波(wave)的意思。

sr 的默认值是22050Hz（人耳识别声音频率的上限就是20000Hz）。如果想保留音频原本的采样率，将sr设为None。也可以根据需求重置采样率。元音和声调的分析11025Hz或16000Hz的采样率就够了，辅音的分析一般设为22050Hz。笔者用22050Hz的采样率录了一个词“文程”，

mono 意为单声道，如果参数值为True，会将不是单声道的音频转换为单声道。一般来说，语音分析用单声道录音就可以了。

offset 和 duration 分别表示音频读取的起始时间点和持续时长，默认值分别是0.0和None，即默认读取整段音频文件。也可以根据需要设置需要读取的音频段。

dtype 和 res_type 分别表示返回值y中元素的数据类型（默认单精度浮点数float32）和重采样(resample)模式，官方文档解释默认的重采样模式kaiser_best是一种高品质的模式，kaiser就是凯撒的意思，听起来就很厉害，应该是不用改的吧。在官方文档提供的示例中，这两个参数都没有改动，估计一般不用改。

2.2 音频加载示例

我用普通话读了一个词“文程”并录音，采样率为22050Hz，音频时长为4.048979591836734秒，音频文件的格式.wav，采用单声道。下面用Librosa加载这个音频。

# 导包
import librosa

# 读取音频
y, sr = librosa.load('文程.wav', offset=1.8, duration=1.7)
print(y)
print(sr)

运行结果：

array([-7.3242188e-04, -7.0190430e-04, -7.3242188e-04, ...,
       -2.4414062e-04, -4.5776367e-04,  3.0517578e-05], dtype=float32)
22050

音频的采样率和声道数跟load()方法参数sr和mono是一样的，这里都保持默认。音频虽然有4秒多，但是前后都有空白，真正有效的大约在1.8秒和3.3秒之间，因此起始点offset设为1.8，持续时间duration设为1.7。

当然也可以把采样率改一下，比如将 sr 调为16000Hz，再来看运行结果：

# 导包
import librosa

# 读取音频
y, sr = librosa.load('文程.wav', sr=16000, offset=1.8, duration=1.7)
print(y)
print(sr)

运行结果：

[-6.1596854e-04 -7.5050717e-04 -7.3633762e-04 ...  6.1120314e-05
  9.0007132e-05  9.4523537e-05]
16000

2.3 简析返回值y

最后简单地解析一下返回值y，在官方文档的描述中y是音频时间序列(audio time series)，实际上就是音频的数字表达，即我们常说的模数(A/D)转换过程。

Librosa设定的numpy数组中元素的默认类型是float32，说明采集的样本大小是32bit的，这个采样精度是很高的。32bit的样本一般采用单精度浮点型(float32)存储，而8bit和16bit一般用整型int.

我们可以看一下y的数组长度，即时间序列中有多少个样本点（代码接着2.2按16000Hz采样的数据）：

print(len(y)) # 27200

16000Hz下1.7秒的音频有27200个采样点。所谓采样率就是每一秒样本的数量，1.7秒的样本数应该为16000*1.7=27200.

这里顺便提一下Librosa的get_duration()方法，该方法获取的是音频的长度（单位：秒）。不过传入的参数不是音频文件，而是时间序列y。

print(librosa.get_duration(y,sr=sr))  # 1.7

3 波形图和语谱图的绘制

Librosa.display中提供了多种音频信号可视化的方法，翻阅源码不难发现display的底层是matplotlib，display本身也常常与matplotlib配合使用，官方文档的示例中也体现了一点。

3.1 绘制波形图

3.1.1 waveshow()方法

目前的0.8.1中有两个方法可用于波形图的绘制，分别是waveplot()和waveshow()。查阅官方文档的更新日志，waveshow()是0.8.1版本中新增，waveplot()将会在未来的0.9.0版本中被删除。下面是waveshow()方法的参数：

librosa.display.waveshow(y, sr=22050, max_points=11025, x_axis='time', 
offset=0.0, marker='', where='post', label=None, ax=None, **kwargs)

前两个参数y和sr跟load()方法的返回值是一样的，因此直接与其保持一致就可以了。

max_points 指的是用于绘图的最大样本数，默认值是采样率的一半。官方文档解释了一通，大概说小于默认值图能画出来，大于默认值就要先降采样再画图。其实这个参数基本是不用动的，原因是奈圭斯特定律(Nyquist's law)。我是这么理解的，声波在时间序列上是一个个周期组成的，每个周期至少两个采样点才能得到完整的信息，那么采样率就是频率的两倍；那么反过来，波形图是对声波的可视化，把一个周期看成一个样本，波形图的样本数就是采样率的一半。

x_axis 顾名思义就是横轴的x轴的标记内容，默认是'time'，x轴单位是秒。除了None以外，有两类六种参数可选。一个是time类，包括'time','s','ms'，x轴单位是秒或毫秒；另一个是lag类，包括'lag','lag_s','lag_ms'，单位也是秒，lag表示滞后时间，从0开始，过了中点为负数。一般采用默认的'time'就可以了。

offset 用于设定音频需要可视化部分的起始时间，一般也不用改。因为在上一步用load()读取音频时就已经把需要处理的音频段设置好了。

后面几个参数都与图像的绘制有关，都是传自matplotlib里面的参数。最后还有一个*kwargs（可变关键字参数），传自matplotlib.pyplot.fill_between 和 matplotlib.pyplot.step，比如我们可以传一个color设置图像的颜色（图像默认是蓝色的）。

3.1.2 波形图绘制示例

代码接着上文2.2按16000Hz的采样率的 load() 方法获取到的 y 和 sr。

在Jupyter Notebook中，只要一行代码就能显示想要的波形图，非常方便：

librosa.display.waveplot(y, sr)

在pycharm等编译器中，需要借助matplotlib.pyplot，否则无法显示图像，示例如下：

plt.figure()
librosa.display.waveplot(y, sr)
plt.show()

我们也可以自由地定制自己喜欢的可视化样式，比如给图改个颜色，再加个标题：

plt.figure()
librosa.display.waveshow(y, sr, color='orange')
plt.title('文程', fontproperties="SimSun")
plt.show()

 改变x轴的范围可以达到局部放大的效果，比如0-0.83秒差不多是“文”这个字音的范围，可以截出来：

plt.figure()
librosa.display.waveshow(y, sr, color='orange')
plt.title('文', fontproperties="SimSun")
plt.xlim(0, 0.83)
plt.show()

 waveshow()方法除了能够显示上面演示的这种包络图(envelope view)，还能显示谐波和脉冲成分(harmonic and percussive components)，官方文档中有相应示例，谐波和脉冲的分析留待后续探索。

3.2 绘制语谱图

语谱图(spectrum)是一种汇集了频率、时间和振幅三种信息的图像，横轴代表时间，纵轴代表频率，频率线的颜色深浅代表能量（反映振幅）。

3.2.1 specshow()方法

librosa.display.specshow(data, x_coords=None, y_coords=None, x_axis=None, y_axis=None, 
sr=22050, hop_length=512, fmin=None, fmax=None, tuning=0.0, bins_per_octave=12, 
key='C:maj', Sa=None, mela=None, thaat=None, auto_aspect=True, htk=False, ax=None, 
**kwargs)

specshow()方法的参数是比较多的，但是一般的需求中要改动的参数并不多。绘制语谱图前要进行频谱分析，即进行时频转换。参数的取值与转换后所得数据的频率尺度(scale)有关。 

频谱分析都是基于短时傅里叶变换(STFT)得到的功率谱(power spectrum)，频率可以采用多种尺度，不管采用哪种尺度的频谱，都要调用librosa.amplitude_to_db()方法将STFT之后得到的幅度值进行对数变换得到分贝值(dB)，在将这个方法的返回值数据传入specshow()方法中，于是就得到了用颜色深浅表示能量大小的显示。

这里演示三种尺度：原尺度（线性linear）、对数尺度(logarithmic scale)和梅尔尺度(mel scale)。另外，笔者在使用中发现采样率在22050Hz以下时得语谱图的图像分辨率并不高，44100Hz的图像效果比较理想，因此重新录了“文程2.wav”，采样率提高到44100Hz。

另外，如果想把多个语谱图画在一起，需要指定specshow()方法中的 ax 参数，这里不做演示。

3.2.2 线性频率的语谱图

将STFT得到的功率谱可视化的实现比较简单，大部分参数直接采用默认值就可以了。其他参数说明见注释。

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

# 读取音频，采样率为44100Hz，持续时间为2秒
y, sr = librosa.load('文程2.wav', offset=0.83, duration=2, sr=None)

# 将 stft 之后的 幅度值的绝对值 转换为 分贝值，将返回值传入specshow()方法中
data = librosa.amplitude_to_db(np.abs(librosa.stft(y)), ref=np.max)

# 绘制图像
fig, ax = plt.subplots(1,1)
# x轴是时间（单位：秒），y轴是由fft窗口和采样率决定的频率值（单位：Hz）
img = librosa.display.specshow(data, sr=sr, x_axis='time', y_axis='linear')
plt.ylim(0, 8500) # 8000Hz以上没有能量显示，因此y轴上限设为8500
plt.title('文程2.wav 线性频率语谱图', fontproperties="SimSun")
fig.colorbar(img, ax=ax, format="%+2.f dB")
plt.show()

3.2.3 对数频率的语谱图

这种语谱图基于对数尺度的功率谱。根据官方文档，首先要把短时傅里叶变换的方法stft()的参数值hop_length(帧移)调大。翻阅stft()方法的源码，知道hop_length指相邻STFT列的音频样本数，在不改变其他参数的情况下是512（stft()方法里的参数关系有点复杂，后续再探索），这里扩大到两倍到1024。

但是得到对数频率下语谱图的关键不在上面的步骤，而是改变y_axis参数，将其设为'log'，如果设为'linear'，得到的图像跟上一节展示的差不多。画出来的图画质简直一言难尽（太糊了，官方文档的图也糊），暂时不知道有什么用，留待今后探索吧。

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

# 读取音频，采样率为44100Hz，持续时间为2秒
y, sr = librosa.load('文程2.wav', offset=0.83, duration=2, sr=None)

# 扩大 帧移 到1024，并传入stft()和specshow()方法中
hop_length = 1024

# 将 stft 之后的 幅度值的绝对值 转换为 分贝值，将返回值传入specshow()方法中
data = librosa.amplitude_to_db(np.abs(librosa.stft(y, hop_length=hop_length)), ref=np.max)

# 绘制图像
fig, ax = plt.subplots(1,1)
# x轴是时间（单位：秒），y轴是对数尺度的频率值（单位：Hz）
img = librosa.display.specshow(data, y_axis='log', sr=sr, hop_length=hop_length, x_axis='time')
plt.title('文程2.wav 对数频率语谱图', fontproperties="SimSun")
fig.colorbar(img, ax=ax, format="%+2.f dB")
plt.show()

3.2.4 梅尔频率的语谱图

梅尔频谱在普通的频谱图加上梅尔滤波函数，相比于普通的频率，梅尔频率更能模拟人耳对声音的感知。在网上查找有关Librosa的资料时，发现不少关于提取MFCC(梅尔倒谱系数)特征提取的文章。因此这里再绘制一个基于梅尔频率的语谱图，看看是什么样子的。

在此之前，要先将功率谱转换为梅尔频谱并将其进行对数变换，具体方法见代码及注释，这些方法的操作细节留待后续探索。

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

# 读取音频，采样率为44100Hz，持续时间为2秒
y, sr = librosa.load('文程2.wav', offset=0.83, duration=2, sr=None)

# 提取梅尔频谱
melspec = librosa.feature.melspectrogram(y, sr)

# 将梅尔频率数据做对数变换
logmelspec = librosa.power_to_db(melspec)

# 绘制图像
fig, ax = plt.subplots(1,1)
# x轴是时间（单位：秒），y轴是梅尔尺度的频率值（单位：Hz）
img = librosa.display.specshow(logmelspec, y_axis='mel', sr=sr, x_axis='time')
plt.title('文程2.wav 梅尔频率语谱图', fontproperties="SimSun")
fig.colorbar(img, ax=ax, format="%+2.f dB")
plt.ylim(0, 9000) # 8192Hz以上没有能量显示，因此y轴上限减为9000
plt.show()

4 语谱图颜色的设置（超级好看）

语谱图的图像颜色是可以改的，除了默认的颜色外，还可改成各种各样的颜色，我们接着3.1.2的线性频谱语谱图的代码绘制几幅不同颜色的图看看效果，要设置的参数是cmap，来自matplotlib。画出来的图真心好看。

import librosa.display
import matplotlib.pyplot as plt
import numpy as np
y, sr = librosa.load('文程2.wav', offset=0.83, duration=2, sr=None)
data = librosa.amplitude_to_db(np.abs(librosa.stft(y)), ref=np.max)
fig = plt.figure(figsize=(12, 9))

ax = fig.add_subplot(221)
spec1 = librosa.display.specshow(data, sr=sr, x_axis='time', 
                                 y_axis='linear',ax=ax, cmap='Greys')
ax.set_title('文程2.wav 灰白色', fontproperties="SimSun")
ax.set_ylim(0, 8500)

ax = fig.add_subplot(222)
spec1 = librosa.display.specshow(data, sr=sr, x_axis='time', 
                                 y_axis='linear',ax=ax, cmap='Blues')
ax.set_title('文程2.wav 蓝色', fontproperties="SimSun")
ax.set_ylim(0, 8500)

ax = fig.add_subplot(223)
spec1 = librosa.display.specshow(data, sr=sr, x_axis='time', 
                                 y_axis='linear',ax=ax, cmap='Greens')
ax.set_title('文程2.wav 绿色', fontproperties="SimSun")
ax.set_ylim(0, 8500)

ax = fig.add_subplot(224)
spec1 = librosa.display.specshow(data, sr=sr, x_axis='time', 
                                 y_axis='linear',ax=ax, cmap='Oranges')
ax.set_title('文程2.wav 橙色', fontproperties="SimSun")
ax.set_ylim(0, 8500)

plt.show()