语音信号处理库——Librosa

librosa语音信号处理 - 简书 (jianshu.com)这篇文章说的非常详细，但有一些函数已经荒废了我做了一些补充。

librosa — librosa 0.8.1 documentation官方文档

特征提取流程图：

1.读取语音

1	y,sr = librosa.load(path, sr=22050, mono=True, offset=0.0, duration=None)

参数：

path：文件路径；
sr：采样频率，默认22050，可以用参数’None’表示用原语音自身的采样频率；
mono：是否将音频转换为单声道
offset：表示音频读取的位置，以s为单位
duration：表示读取音频的长度，以s为单位

返回值：

y：音频时间序列
sr：音频的采样频率

import librosa
import matplotlib.pyplot as plt
y,sr= librosa.load(path="test.wav",sr=None,mono=False,offset=1,duration=3)
plt.plot(y)

2.读取时长

1	d = librosa.get_duration(y=None, sr=22050, S=None, n_fft=2048, hop_length=512, center=True, filename=None)

参数：

y ：音频时间序列
sr ：y的音频采样率
S ：STFT矩阵或任何STFT衍生的矩阵（例如，色谱图或梅尔频谱图）。根据频谱图输入计算的持续时间仅在达到帧分辨率之前才是准确的。如果需要高精度，则最好直接使用音频时间序列。
n_fft ：S的 FFT窗口大小
hop_length ：S列之间的音频样本数
center ：布尔值
- 如果为True，则S [:, t]的中心为y [t * hop_length]
- 如果为False，则S [:, t]从y[t * hop_length]开始
filename ：如果提供，则所有其他参数都将被忽略，并且持续时间是直接从音频文件中计算得出的。

返回值：

d：音频时长

librosa.get_duration(y,sr)#承接上面的代码
#out：3.0
librosa.get_duration(filename="test.wav")
#out：6.56

3.采样率

1	sr = librosa.get_samplerate(path)

参数：

path：音频文件的路径

sr：音频文件的采样率

1 2	librosa.get_samplerate("test.wav") #out：8000

4.写音频

1 2	#librosa.output.write_wav(path, y, sr, norm=False)在0.7.0版本荒废了现在要用soundfile.write() soundfile.write(file, data, samplerate, subtype=None, endian=None, format=None, closefd=True)

参数：

file：要保存的文件名
data：音频时间序列
samplerate：音频的采样率
其他参数：参考SoundFile — PySoundFile 0.10.3post1-1-g0394588 documentation

import soundfile as sf
sf.write("test2.wav",y,sr)
y2,sr2 = librosa.load("test2.wav",sr=None)
plt.plot(y2)

5.过零率

1	zcr = librosa.feature.zero_crossing_rate(y, frame_length = 2048, hop_length = 512, center = True)

参数：

y ：音频时间序列
frame_length ：帧长
hop_length ：帧移
center：bool，如果为True，则通过填充y的边缘来使帧居中。

返回：

zcr：zcr[0，i]是第i帧中的过零率

librosa.feature.zero_crossing_rate(y)
#out：array([[0.14208984, 0.1796875 , 0.23486328, 0.20800781, 0.21875   ,
#        0.25830078, 0.22314453, 0.21337891, 0.20263672, 0.18115234,
#        0.17480469, 0.16210938, 0.16015625, 0.15722656, 0.22460938,
#        0.23974609, 0.26367188, 0.30615234, 0.30615234, 0.33691406,
#        0.34375   , 0.3125    , 0.29345703, 0.28417969, 0.27539062,
#        0.29345703, 0.26611328, 0.24511719, 0.21386719, 0.17480469,
#        0.17529297, 0.15039062, 0.17236328, 0.16650391, 0.18310547,
#        0.17675781, 0.14208984, 0.22802734, 0.30566406, 0.40576172,
#        0.52880859, 0.55126953, 0.57080078, 0.59765625, 0.58398438,
#        0.54394531, 0.39501953]])

6.波形图

1 2	librosa.display.waveplot(y, sr=22050, x_axis='time', offset=0.0, ax=None) #需要注意，必须import librosa.display，否则会报错“no attribute ‘display’”

参数：

y ：音频时间序列
sr ：y的采样率
x_axis ：str {‘time’，‘off’，‘none’}或None，如果为“时间”，则在x轴上给定时间刻度线。
offset：水平偏移（以秒为单位）开始波形图

1
2

#librosa.display.waveplot在0.9.0会废除，所以这里推荐使用librosa.display.waveshow
librosa.display.waveplot(y, sr=22050, max_points=50000.0, x_axis='time', offset=0.0, max_sr=1000, ax=None, **kwargs)

参数：

max_points：绘制的最大时间点数，如果超过持续时间，则执行下采样
max_sr：可视化的最大采样率‎

import librosa.display
librosa.display.waveplot(y,sr=8000,x_axis="time",offset=0,ax=None)
plt.figure()
librosa.display.waveshow(y,sr=8000,x_axis="ms")

7.短时傅里叶变换

1	D = librosa.stft(y,n_fft=2048,hop_length=None,win_length=None,window='hann',center=True,pad_mode='reflect')

短时傅立叶变换（STFT），返回一个复数矩阵使得D(f,t)

复数的实部：np.abs(D(f,t))频率的振幅
复数的虚部：np.angle(D(f,t))频率的相位

参数：

y：音频时间序列
n_fft：FFT窗口大小，n_fft=hop_length+overlapping
hop_length：帧移，如果未指定，则默认win_length / 4。
win_length：每一帧音频都由window（）加窗。窗长win_length，然后用零填充以匹配N_FFT。默认win_length=n_fft。
window

：字符串，元组，数字，函数 shape =（n_fft, )
- 窗口（字符串，元组或数字）；
- 窗函数，例如scipy.signal.hanning
- 长度为n_fft的向量或数组
center

：bool
- 如果为True，则填充信号y，以使帧 D [:, t]以y [t * hop_length]为中心。
- 如果为False，则D [:, t]从y [t * hop_length]开始
dtype：D的复数值类型。默认值为64-bit complex复数
pad_mode：如果center = True，则在信号的边缘使用填充模式。默认情况下，STFT使用reflection padding。

返回：

STFT矩阵， $shape=(1+\frac{nfft}{2},t)$

librosa.stft(y,center=True)
'''out：array([[-0.00171932+0.0000000e+00j,  0.00592006+0.0000000e+00j,
         0.01431638+0.0000000e+00j, ..., -0.00448615+0.0000000e+00j,
        -0.00437583+0.0000000e+00j, -0.00784964+0.0000000e+00j],
       [ 0.00244478+5.6835911e-18j, -0.00073931+7.5202361e-03j,
        -0.0131455 +9.1534172e-04j, ...,  0.00395171+4.8493536e-04j,
         0.00013135-3.2658517e-04j,  0.00466984-3.1057880e-03j],
       [-0.00165288+1.0003131e-18j, -0.00410656-1.7896685e-03j,
         0.00912124-1.1786177e-03j, ..., -0.00278488-1.6674915e-03j,
         0.00238109-5.2638666e-04j, -0.0019866 +1.9712506e-03j],
       ...,
       [-0.00029123+9.0889913e-19j, -0.0007377 -1.0746513e-03j,
         0.00117161+1.0120005e-03j, ..., -0.00201914+9.7867288e-04j,
        -0.00083203-9.6789573e-04j,  0.00158986+7.2865405e-05j],
       [-0.00081057+3.4067665e-18j, -0.00016138-3.1721860e-04j,
        -0.00199592-1.2107765e-03j, ...,  0.00122393+5.3408562e-04j,
         0.00035368-9.7684760e-04j, -0.0012544 -7.4694009e-04j],
       [ 0.00112953+0.0000000e+00j,  0.00147533+0.0000000e+00j,
         0.0029277 +0.0000000e+00j, ..., -0.00073284+0.0000000e+00j,
        -0.0003154 +0.0000000e+00j,  0.00141847+0.0000000e+00j]],
      dtype=complex64)'''

8.短时傅里叶逆变换

1	y = librosa.istft(stft_matrix, hop_length=None, win_length=None, window='hann', center=True, length=None)

短时傅立叶逆变换（ISTFT），将复数值D(f,t)频谱矩阵转换为时间序列y，窗函数、帧移等参数应与stft相同

参数：

stft_matrix ：经过STFT之后的矩阵
hop_length ：帧移，默认为 $\frac{𝑤𝑖𝑛𝑙𝑒𝑛𝑔𝑡ℎ}{4}$
win_length ：窗长，默认为n_fft
window

：字符串，元组，数字，函数或shape = (n_fft, )
- 窗口（字符串，元组或数字）
- 窗函数，例如scipy.signal.hanning
- 长度为n_fft的向量或数组
center：bool
- 如果为True，则假定D具有居中的帧
- 如果False，则假定D具有左对齐的帧
length：如果提供，则输出y为零填充或剪裁为精确长度音频

返回：

y ：时域信号

9.幅度转dB

1	librosa.amplitude_to_db(S, ref=1.0)

将幅度频谱转换为dB标度频谱。也就是对S取对数。与这个函数相反的是librosa.db_to_amplitude(S)

参数：

S ：输入幅度
ref ：参考值，振幅abs（S）相对于ref进行缩放， $20∗𝑙𝑜𝑔_{10}(\frac{𝑆}{𝑟𝑒𝑓})$

返回：

dB为单位的S

10.功率转dB

1	librosa.core.power_to_db(S, ref=1.0)

将功率谱（幅度平方）转换为分贝（dB）单位，与这个函数相反的是librosa.db_to_power(S)

参数：

S：输入功率
ref ：参考值，振幅abs(S)相对于ref进行缩放， $10∗𝑙𝑜𝑔_{10}(\frac{𝑆}{𝑟𝑒𝑓})$

返回：

dB为单位的S

11.语谱图

1	librosa.display.specshow(data, x_axis=None, y_axis=None, sr=22050, hop_length=512)

参数：

data：要显示的矩阵
sr ：采样率
hop_length ：帧移
x_axis 、y_axis ：x和y轴的范围
频率类型
- ‘linear’，‘fft’，‘hz’：频率范围由FFT窗口和采样率确定
- ‘log’：频谱以对数刻度显示
- ‘mel’：频率由mel标度决定
时间类型
- time：标记以毫秒，秒，分钟或小时显示。值以秒为单位绘制。
- s：标记显示为秒。
- ms：标记以毫秒为单位显示。
所有频率类型均以Hz为单位绘制

#窗长大，窄带语谱图
D = librosa.amplitude_to_db(np.abs(librosa.stft(y)), ref=np.max)
plt.subplot(2, 1, 1)
librosa.display.specshow(D, y_axis='linear')
plt.colorbar(format='%+2.0f dB')
plt.title('Linear-frequency power spectrogram')

plt.subplot(2, 1, 2)
librosa.display.specshow(D, y_axis='log')
plt.colorbar(format='%+2.0f dB')
plt.title('Log-frequency power spectrogram')
plt.show()

#窗长减小，则是宽带语谱图
D = librosa.amplitude_to_db(np.abs(librosa.stft(y,win_length=64)), ref=np.max)
plt.subplot(2, 1, 1)
librosa.display.specshow(D, y_axis='linear')
plt.colorbar(format='%+2.0f dB')
plt.title('Linear-frequency power spectrogram')

plt.subplot(2, 1, 2)
librosa.display.specshow(D, y_axis='log')
plt.colorbar(format='%+2.0f dB')
plt.title('Log-frequency power spectrogram')
plt.show()

12.Mel滤波器组

1	librosa.filters.mel(sr, n_fft, n_mels=128, fmin=0.0, fmax=None, htk=False, norm=1)

创建一个滤波器组矩阵以将FFT合并成Mel频率

参数：

sr ：输入信号的采样率
n_fft ：FFT组件数
n_mels ：产生的梅尔带数
fmin ：最低频率（Hz）
fmax：最高频率（以Hz为单位）。如果为None，则使用fmax = sr / 2.0
norm

：{None，1，np.inf} [标量]
- 如果为1，则将三角mel权重除以mel带的宽度（区域归一化）。否则，保留所有三角形的峰值为1.0

返回：

Mel变换矩阵

melfb = librosa.filters.mel(22050, 2048)
# array([[ 0.   ,  0.016, ...,  0.   ,  0.   ],
#        [ 0.   ,  0.   , ...,  0.   ,  0.   ],
#        ...,
#        [ 0.   ,  0.   , ...,  0.   ,  0.   ],
#        [ 0.   ,  0.   , ...,  0.   ,  0.   ]])
import matplotlib.pyplot as plt
plt.figure()
librosa.display.specshow(melfb, x_axis='linear')
plt.ylabel('Mel filter')
plt.title('Mel filter bank')
plt.colorbar()
plt.tight_layout()
plt.show()

13.计算Mel scaled 频谱

1	librosa.feature.melspectrogram(y=None, sr=22050, S=None, n_fft=2048, hop_length=512, win_length=None, window='hann', center=True, pad_mode='reflect', power=2.0)

如果提供了频谱图输入S，则通过mel_f.dot（S）将其直接映射到mel_f上。

如果提供了时间序列输入y，sr，则首先计算其幅值频谱S，然后通过mel_f.dot（S ** power）将其映射到mel scale上。默认情况下，power= 2在功率谱上运行。

参数：

**y **：音频时间序列
**sr **：采样率
**S **：频谱
**n_fft **：FFT窗口的长度
**hop_length **：帧移
**win_length **：窗口的长度为win_length，默认win_length = n_fft
**window **：字符串，元组，数字，函数或shape =（n_fft, )
- 窗口规范（字符串，元组或数字）；看到scipy.signal.get_window
- 窗口函数，例如 scipy.signal.hanning
- 长度为n_fft的向量或数组
center

：bool
- 如果为True，则填充信号y，以使帧 t以y [t * hop_length]为中心。
- 如果为False，则帧t从y [t * hop_length]开始
power：幅度谱的指数。例如1代表能量，2代表功率，等等
n_mels：滤波器组的个数 1288
fmax：最高频率

Mel频谱shape=(n_mels, t)

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

y, sr = librosa.load(librosa.util.example_audio_file())
# 方法一：使用时间序列求Mel频谱
print(librosa.feature.melspectrogram(y=y, sr=sr))
'''[[2.9641967e-07 2.5462259e-03 6.9437072e-02 ... 3.4732797e-09
  5.4255165e-09 1.8523759e-09]
 [1.7560427e-07 1.1623162e-02 1.3323289e+00 ... 2.8872930e-07
  7.6875637e-08 4.4768003e-09]
 [3.8013613e-07 2.3408584e-02 5.0721464e+00 ... 1.6204266e-06
  3.2566246e-07 3.5224299e-09]
 ...
 [2.3962335e-10 2.7136402e-08 2.1219030e-06 ... 2.4645463e-11
  3.0509512e-12 8.6144919e-15]
 [1.4714385e-10 9.5965706e-09 6.5881960e-07 ... 5.2040217e-12
  5.4942682e-13 1.4356046e-15]
 [3.0418751e-12 2.1499022e-10 4.4302659e-08 ... 4.5521724e-13
  9.9850961e-14 5.5009496e-16]]'''

# 方法二：使用stft频谱求Mel频谱
D = np.abs(librosa.stft(y)) ** 2  # stft频谱
S = librosa.feature.melspectrogram(S=D)  # 使用stft频谱求Mel频谱

plt.figure(figsize=(10, 4))
librosa.display.specshow(librosa.power_to_db(S, ref=np.max),
                         y_axis='mel', fmax=8000, x_axis='time')
plt.colorbar(format='%+2.0f dB')
plt.title('Mel spectrogram')
plt.tight_layout()
plt.show()

14.提取Log-Mel Spectrogram 特征(FBank特征)

import librosa

y, sr = librosa.load(librosa.util.example_audio_file(), sr=16000)
# 提取 mel spectrogram feature
melspec = librosa.feature.melspectrogram(y, sr, n_fft=1024, hop_length=512, n_mels=128)
logmelspec = librosa.amplitude_to_db(melspec)        # 转换到对数刻度

print(logmelspec.shape)        # (128, 65)

可见，Log-Mel Spectrogram特征是二维数组的形式，128表示Mel频率的维度（频域），64为时间帧长度（时域），所以Log-Mel Spectrogram特征是音频信号的时频表示特征。其中，n_fft指的是窗的大小，这里为1024；hop_length表示相邻窗之间的距离，这里为512，也就是相邻窗之间有50%的overlap；n_mels为mel bands的数量，这里设为128。

15.提取MFCC系数

1	librosa.feature.mfcc(y=None, sr=22050, S=None, n_mfcc=20, dct_type=2, norm='ortho', **kwargs)

参数：

y：音频数据
sr：采样率
S：np.ndarray，对数功能梅尔谱图
n_mfcc：int>0，要返回的MFCC数量
dct_type：None, or {1, 2, 3} 离散余弦变换（DCT）类型。默认情况下，使用DCT类型2。
norm： None or ‘ortho’ 规范。如果dct_type为2或3，则设置norm =’ortho’使用正交DCT基础。标准化不支持dct_type = 1。

返回：

M： MFCC序列

import librosa

y, sr = librosa.load(librosa.util.example_audio_file(), sr=16000)
# 提取 MFCC feature
mfccs = librosa.feature.mfcc(y=y, sr=sr, n_mfcc=40)

print(mfccs.shape)        # (40, 1921)