2021-04-17

特征提取背景

在之前的Task2中，我们初步查看了数据集的信息，并通过音频的声波图和声谱图将声音可视化了。当我们拿到这些音频数据之后，接下来就需要进行特征提取（过滤掉背景噪音等不需要的信息）筛选出我们需要的信息了。

接下来我们将简要介绍以下特征，并详细学习MFCC特征提取知识：

过零率 （Zero Crossing Rate）
频谱质心 （Spectral Centroid）
声谱衰减 (Spectral Roll-off）
梅尔频率倒谱系数 （Mel-frequency cepstral coefficients ，MFCC）
色度频率 （Chroma Frequencies）

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

2.1 过零率

过零率（zero crossing rate）是一个信号符号变化的比率，即，在每帧中，语音信号从正变为负或从负变为正的次数。 这个特征已在语音识别和音乐信息检索领域得到广泛使用，通常对类似金属、摇滚等高冲击性的声音的具有更高的价值。

一般情况下，过零率越大，频率近似越高

x, sr = librosa.load('./train_sample/aloe/24EJ22XBZ5.wav')
#绘制声波信号
plt.figure(figsize=(14, 5))
librosa.display.waveplot(x, sr=sr)#默认采样率（sr
# 放大
n0 = 9000
n1 = 9100
plt.figure(figsize=(14, 5))
plt.plot(x[n0:n1])
plt.grid()

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-zuCHPXyb-1618670530959)(output_3_0.png)]

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-0cHkyKuz-1618670530963)(output_3_1.png)]

# 计算过零率
zero_crossings = librosa.zero_crossings(x[n0:n1], pad=False)
print(sum(zero_crossings))

15

2.2 频谱质心

谱质心(Spectral Centroid)是描述音色属性的重要物理参数之一，是频率成分的重心，是在一定频率范围内通过能量加权平均的频率，其单位是Hz。它是声音信号的频率分布和能量分布的重要信息。在主观感知领域，谱质心描述了声音的明亮度，具有阴暗、低沉品质的声音倾向有较多低频内容，谱质心相对较低，具有明亮、欢快品质的多数集中在高频，谱质心相对较高。该参数常用于对乐器声色的分析研究。

spectral_centroids = librosa.feature.spectral_centroid(x, sr=sr)[0]
print(spectral_centroids.shape)
# (2647,)
# 计算时间变量 
frames = range(len(spectral_centroids))
t = librosa.frames_to_time(frames)
# 归一化频谱质心
def normalize(x, axis=0):
    return sklearn.preprocessing.minmax_scale(x, axis=axis)
#沿波形绘制频谱质心 
librosa.display.waveplot(x, sr=sr, alpha=0.4)
plt.plot(t, normalize(spectral_centroids), color='r')

(140,)





[<matplotlib.lines.Line2D at 0x2511aac7d60>]

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-mCCMLqXg-1618670530965)(output_6_2.png)]

2.3 声谱衰减
它是对声音信号形状（波形图）的一种衡量，表示低于总频谱能量的指定百分比的频率。

spectral_rolloff = librosa.feature.spectral_rolloff(x+0.01, sr=sr)[0]
librosa.display.waveplot(x, sr=sr, alpha=0.4)
plt.plot(t, normalize(spectral_rolloff), color='b')

[<matplotlib.lines.Line2D at 0x2511ab23fd0>]

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-H1U0pLmQ-1618670530968)(output_8_1.png)]

hop_length = 512
chromagram = librosa.feature.chroma_stft(x, sr=sr, hop_length=hop_length)
plt.figure(figsize=(15, 5))
librosa.display.specshow(chromagram, x_axis='time', y_axis='chroma', hop_length=hop_length, cmap='coolwarm')

<matplotlib.collections.QuadMesh at 0x2511ab88520>

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-5D3XtOlk-1618670530970)(output_9_1.png)]

03 MFCC特征提取

人的耳朵在接收信号的时候，不同的频率会引起耳蜗不同部位的震动。耳蜗就像一个频谱仪，自动在做特征提取并进行语音信号的处理。在语音识别领域中MFCC（Mel Frequency Cepstral Coefficents）特征提取是最常用的方法，也是本次音频分类任务中涉及到的特征提取方法。

具体来说，MFCC特征提取的步骤如下：

对语音信号进行分帧处理
用周期图(periodogram)法来进行功率谱(power spectrum)估计
对功率谱用Mel滤波器组进行滤波，计算每个滤波器里的能量
对每个滤波器的能量取log
进行离散余弦变换（DCT）变换
保留DCT的第2-13个系数，去掉其它
其中，前面两步是短时傅里叶变换，后面几步主要涉及梅尔频谱，接下来我们将分别学习这些知识。

3.1 短时傅里叶分析

声音信号本是一维时域信号（声音信号随时间变化），我们可以通过傅里叶变换将其转换到频域上,但这样又失去了时域信息，无法看出频率分布随时间的变化。短时傅里叶（STFT）就是为了解决这个问题而发明的常用手段。

所谓的短时傅里叶变换，即把一段长信号分帧、加窗，再对每一帧做快速傅里叶变换（FFT），最后把每一帧的结果沿另一个维度堆叠起来，得到类似于一幅图的二维信号形式，也就是我们task2中得到的声谱图。

3.1.1 分帧
语音信号是不稳定的时变信号，但为了便于处理，我们假设在一个很短的时间内，如20-40ms内为一个稳定的系统，也就是1帧。但是我们不能简单平均分割语音，相邻的帧之间需要有一定的重合。我们通常以25ms为1帧，帧移为10ms，因此1秒的信号会有10帧。

3.1.2 对每帧信号进行DFT
我们可以使用离散傅里叶变换(DFT)把每一帧信号变换到时域，公式是：

图片.png

其中s(n)表示时域信号；si(n)是第i帧的数据，其中n的范围是1-400；当我们介绍DFT的时候，Si(k)表示的是第i帧的第k个复系数；Pi(k)是第i帧的功率谱。h(n)是一个N点的窗函数(比如Hamming窗)，K是DFT的长度。有了Si(k)我们就可以估计功率谱：

图片.png

这样，我们就得到了周期图的功率谱估计，即声谱图。以之前task2的芦荟音频为例，它的横坐标是帧下标，纵坐标是不同频率或者叫能量，图中颜色越深(比如红色)，对应频率的能量越大。

图片.png

我们可以利用 librosa库做短时傅立叶变换（STFT），通过STFT能够返回一个复数矩阵，使得：

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

# STFT
y, sr = librosa.load('./train_sample/aloe/24EJ22XBZ5.wav')
S =librosa.stft(y, n_fft=2048, hop_length=None, win_length=None, window='hann', center=True, pad_mode='reflect')
'''
• 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
'''
S = np.abs(S)
print(y.shape)
print(S.shape)

(71661,)
(1025, 140)

04 梅尔频谱和梅尔倒谱
声谱图往往是很大的一张图，且依旧包含了大量无用的信息，所以我们需要通过梅尔标度滤波器组（mel-scale filter banks）将其变为梅尔频谱。

4.1 梅尔尺度
梅尔尺度（Mel Scale）是建立从人类的听觉感知的频率——Pitch到声音实际频率直接的映射。频率的单位是赫兹（Hz），人耳能听到的频率范围是20-20000Hz，但人耳对Hz这种标度单位并不是线性感知关系，例如，若把音调频率从1000Hz提高到2000Hz，我们的耳朵只能觉察到频率似乎提高了一些而不是一倍。但是通过把频率转换成美尔尺度，我们的特征就能够更好的匹配人类的听觉感知效果。从频率到梅尔频率的转换公式如下：

图片.png

我们可以观察一下转换后的映射图，可以发现人耳对于低频声音的分辨率要高于高频的声音，因为赫兹到梅尔是log的关系，所以当频率较小时，mel随Hz变化较快；当频率很大时，mel的上升很缓慢，曲线的斜率很小。这说明了人耳对低频音调的感知较灵敏，在高频时人耳是很迟钝的，梅尔标度滤波器组启发于此。

图片.png

4.2 梅尔滤波器
为了模拟人耳对声音的感知，人们发明的梅尔滤波器组。一组大约20-40(通常26)个三角滤波器组，它会对上一步得到的周期图的功率谱估计进行滤波。而且区间的频率越高，滤波器就越宽(但是如果把它变换到美尔尺度则是一样宽的)。为了计算方便，我们通常把26个滤波器用一个矩阵来表示，这个矩阵有26行，列数就是傅里叶变换的点数。

图片.png

计算过程如下图所示，最后我们会保留这26个滤波器的能量。图(a)是26个滤波器；图(b)是滤波后的信号；图©是其中的第8个滤波器，它只让某一频率范围的信号通过；图(d)通过它的信号的能量；图(e)是第20个滤波器；图(f)是通过它的信号的能量。

图片.png

4.3 梅尔倒谱
在梅尔频谱上做倒谱分析（取对数log，做离散余弦变换（DCT）变换）就得到了梅尔倒谱。

对上面得到的26个点的信号进行DCT，得到26个倒谱系数(Cepstral Coefficents)，最后我们保留2-13这12个数字，这12个数字就叫MFCC特征。对功率谱再做DCT的目的就是为了提取信号的包络。

#梅尔频率倒谱系数 MFCC
mfccs = librosa.feature.mfcc(x, sr)
print (mfccs.shape)
# (20, 151)
#Displaying  the MFCCs:
librosa.display.specshow(mfccs, sr=sr, x_axis='time')

(20, 140)





<matplotlib.collections.QuadMesh at 0x2511adc55b0>

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-0OlLzlwk-1618670530971)(output_14_2.png)]

# mfcc 特征缩放
mfccs = sklearn.preprocessing.scale(mfccs, axis=1)
print(mfccs.mean(axis=1))
print(mfccs.var(axis=1))
librosa.display.specshow(mfccs, sr=sr, x_axis='time')

[ 0.0000000e+00  0.0000000e+00  1.3623919e-08 -3.4059797e-09
  6.8119594e-09  1.3623919e-08  1.3623919e-08  1.3623919e-08
  6.8119594e-09  0.0000000e+00  4.2574748e-09 -2.9802323e-09
  6.8119594e-09 -1.3623919e-08  0.0000000e+00 -8.5149496e-09
  5.1089697e-09 -3.4059797e-09  0.0000000e+00 -6.8119594e-09]
[1.         1.         1.         1.0000001  1.         0.99999976
 1.0000002  1.         1.         0.9999999  1.         1.0000001
 1.         0.9999999  1.         0.9999999  1.         1.
 0.99999976 0.9999999 ]


C:\ProgramData\Anaconda3\lib\site-packages\sklearn\preprocessing\_data.py:194: UserWarning: Numerical issues were encountered when centering the data and might not be solved. Dataset may contain too large values. You may need to prescale your features.
  warnings.warn("Numerical issues were encountered "
C:\ProgramData\Anaconda3\lib\site-packages\sklearn\preprocessing\_data.py:211: UserWarning: Numerical issues were encountered when scaling the data and might not be solved. The standard deviation of the data is probably very close to 0. 
  warnings.warn("Numerical issues were encountered "





<matplotlib.collections.QuadMesh at 0x2511ae40280>

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-5ptHZBCU-1618670530974)(output_15_3.png)]

label = []
# 建立类别标签，不同类别对应不同的数字。
label_dict = {'aloe': 0, 'burger': 1, 'cabbage': 2,'candied_fruits':3, 'carrots': 4, 'chips':5,
                  'chocolate': 6, 'drinks': 7, 'fries': 8, 'grapes': 9, 'gummies': 10, 'ice-cream':11,
                  'jelly': 12, 'noodles': 13, 'pickles': 14, 'pizza': 15, 'ribs': 16, 'salmon':17,
                  'soup': 18, 'wings': 19}
label_dict_inv = {v:k for k,v in label_dict.items()}

from tqdm import tqdm
def extract_features(parent_dir, sub_dirs, max_file=10, file_ext="*.wav"):
    c = 0
    label, feature = [], []
    for sub_dir in sub_dirs:
        for fn in tqdm(glob.glob(os.path.join(parent_dir, sub_dir, file_ext))[:max_file]): # 遍历数据集的所有文件
            
           # segment_log_specgrams, segment_labels = [], []
            #sound_clip,sr = librosa.load(fn)
            #print(fn)
            label_name = fn.split('/')[-2]
            label.extend([label_dict[label_name]])
            X, sample_rate = librosa.load(fn,res_type='kaiser_fast')
            mels = np.mean(librosa.feature.melspectrogram(y=X,sr=sample_rate).T,axis=0) # 计算梅尔频谱(mel spectrogram),并把它作为特征
            feature.extend([mels])
            
    return [feature, label]

temp = np.array(temp)
data = temp.transpose()

---------------------------------------------------------------------------

NameError                                 Traceback (most recent call last)

<ipython-input-17-16088c198d1d> in <module>
----> 1 temp = np.array(temp)
      2 data = temp.transpose()


NameError: name 'temp' is not defined

# 获取特征
X = np.vstack(data[:, 0])

# 获取标签
Y = np.array(data[:, 1])
print('X的特征尺寸是：',X.shape)
print('Y的特征尺寸是：',Y.shape)

---------------------------------------------------------------------------

NameError                                 Traceback (most recent call last)

<ipython-input-18-da33bccf2eb6> in <module>
      1 # 获取特征
----> 2 X = np.vstack(data[:, 0])
      3 
      4 # 获取标签
      5 Y = np.array(data[:, 1])


NameError: name 'data' is not defined

# 在Keras库中：to_categorical就是将类别向量转换为二进制（只有0和1）的矩阵类型表示
Y = to_categorical(Y)

---------------------------------------------------------------------------

NameError                                 Traceback (most recent call last)

<ipython-input-19-7fbd1762e890> in <module>
      1 # 在Keras库中：to_categorical就是将类别向量转换为二进制（只有0和1）的矩阵类型表示
----> 2 Y = to_categorical(Y)


NameError: name 'to_categorical' is not defined

# 最终数据
print(X.shape)
print(Y.shape)

---------------------------------------------------------------------------

NameError                                 Traceback (most recent call last)

<ipython-input-20-aab608e6d1ba> in <module>
      1 # 最终数据
----> 2 print(X.shape)
      3 print(Y.shape)


NameError: name 'X' is not defined