手写Fbank语音特征提取

语音特征-Fbank的绘制

Fbank提取过程如下图所示：

导入需要的包

import numpy as np
import numpy
import scipy.io.wavfile
from scipy.fftpack import dct
import matplotlib.pyplot as plt
import soundfile

1、读取语音信号

#推荐单声道、16k的音频
signal, sample_rate = soundfile.read(r'test.wav')  
signal = signal[0:int(2 * sample_rate)]  # 只取前2s的音频进行展示

2、预加重

# 预加重
pre_emphasis = 0.97
emphasized_signal = numpy.append(signal[0], signal[1:] - pre_emphasis * signal[:-1])

3、分帧

如果采样率为16k，
帧长为25ms，则每一帧有400个点；
帧移为10ms，每次移动160个点。
num_frames 代表的是最终取得帧数，设音频采样率为16k，如果裁剪为2s的音频，则长度为length = num_frames * 160 + 240，这里num_frames 代表的就是编号，240类似于d段

frame_size = 0.025  # 采样率16k时，采样点：400
frame_stride = 0.01  # 采样率16k时，采样点：160
frame_length, frame_step = frame_size * sample_rate, frame_stride * sample_rate  # 从秒转换为采样点
# 信号长度
signal_length = len(emphasized_signal)
frame_length = int(round(frame_length))
frame_step = int(round(frame_step))  # 四舍五入

num_frames = 200 
length = num_frames * 160 + 240
if emphasized_signal.shape[0] <= length:
    shortage = length - emphasized_signal.shape[0]
    audio = numpy.pad(emphasized_signal, (0, shortage), 'wrap')
# 随机选择音频内的2s片段
start_frame = numpy.int64(random.random() * (audio.shape[0] - length))
pad_signal = audio[start_frame:start_frame + length]

4、加窗

indices = numpy.tile(numpy.arange(0, frame_length), (num_frames, 1)) + numpy.tile(numpy.arange(0, num_frames * frame_step, frame_step), (frame_length, 1)).T
# 把每个采样点的值，按照索引的矩阵进行填充
frames = pad_signal[indices.astype(numpy.int32, copy=False)]
frames *= numpy.hamming(frame_length)
print("frames shape:", frames.shape)

np.tile(data,(x,y)) 扩充data 到x行，y列。
indices的每一行代表了一帧的每个采样点的索引，共有多少行就是有多少帧

frames就是通过索引的位置，把每一个点的值填充上

5、FFT

NFFT = 512  # 傅里叶变换的点数，
# 如果 n 小于输入的长度，则裁剪输入。如果它更大，则用零填充输入。如果未给出 n，则使用沿轴指定的轴的输入长度
# 求绝对值，如果是负数，就是求他的模长
`mag_frames = numpy.absolute(numpy.fft.rfft(frames, NFFT))  # fft的幅度(magnitude)`mag_frames = numpy.absolute(numpy.fft.rfft(frames, NFFT))  # fft的幅度(magnitude)

6、功率谱

pow_frames = ((1.0 / NFFT) * ((mag_frames) ** 2))  # 功率谱

7、mel滤波

nfilt = 80  # 滤波器的个数
low_freq_mel = 0
'''
 为什么采样率除以2就是最高频率？
 奈奎斯特采样定理：为了准确地在数字域中表示信号，采样率应至少是信号中最高频率分量的两倍。
'''
high_freq_mel = (2595 * np.log10(1 + (sample_rate / 2) / 700))  # 将Hz转换为Mel

```python
# 需要80个滤波器组，为此需要82个点，这意味着需要low_freq_mel和high_freq_mel之间线性间隔80个点
mel_points = np.linspace(low_freq_mel, high_freq_mel, nfilt + 2)  # 使得Mel scale间距相等
hz_points = (700 * (10 ** (mel_points / 2595) - 1))  # 将Mel转换回Hz
print("hz_points", hz_points)

# 各个mel滤波器中心点对应FFT的区域编码，找到有值的位置
bins = np.floor((NFFT + 1) * hz_points / sample_rate)  # np.floor向下取整
# bins = (hz_points / (sample_rate / 2)) * (NFFT / 2)  # 各个mel滤波器中心点对应FFT的区域编码，找到有值的位置
print("bins: ", bins, bins.shape)

fbank = np.zeros((nfilt, int(np.floor(NFFT / 2 + 1))))
# 此处计算每一个点位的能量值是多少
for m in range(1, nfilt + 1):
    f_m_left = int(bins[m - 1])  # 左
    f_m_center = int(bins[m])  # 中
    f_m_right = int(bins[m + 1])  # 右
    for k in range(f_m_left, f_m_center):
        # 相似三角形
        fbank[m - 1, k] = (k - bins[m - 1]) / (bins[m] - bins[m - 1])
    for k in range(f_m_center, f_m_right):
        fbank[m - 1, k] = (bins[m + 1] - k) / (bins[m + 1] - bins[m])
filter_banks = np.dot(pow_frames, fbank.T)

下面就是绘制谱图：

# np.finfo(float).eps 返回的是 float 类型的最小正浮点数，通常为 2.220446049250313e-16
# 该代码的目的是确保数组中不会存在0值，以避免除以0等运算错误（因为下一行有log运算）。
filter_banks = np.where(filter_banks == 0, np.finfo(float).eps, filter_banks)  # 数值稳定性
filter_banks = 20 * np.log10(filter_banks)  # dB

# 绘制谱图
y = filter_banks.T
plt.imshow(y, cmap='jet', origin='lower')
plt.colorbar()
plt.show()
plt.clf()

plt.imshow(y[::-1, :], cmap='gray')
plt.axis(False)
plt.show()

参考文献：
https://blog.csdn.net/wudibaba21/article/details/108863431
https://blog.csdn.net/chumingqian/article/details/124950613

可能出现的问题

当输入的音频采样率不是16k，不是单声道的时候，会出现绘制出来的谱图能量较大