【语音识别】作业1：语音特征提取
语音识别原理与应用：语音特征提取（实践）
-
1. 基于标准WAV头文件设计代码以实现PCM WAV文件的读写功能。
-
2. 完成FBank、MFCC及PLP三种声学特性的代码构建。
-
3. 针对MFCC技术展开相关问题分析：
-
- 3.1 探讨采样率、帧长、帧移参数如何影响MFCC向量数量。
-
- 3.2 展述FFT长度与每帧数据点数量之间的相互影响。
-
- 3.3 详细说明Mel尺度频率转换的具体计算流程。
-
- 3.4 深入探讨DCT变换后生成的MFCC静态特征特性。
-
- 3.5 细胞分析一阶和二阶动态特性的计算方法及其意义。
-
4. 深入探讨STFT系列声学特性和CQCC特性在频谱特性分布方面的差异。
- 5. 当语音模拟信号以16kHz进行采样后得到的离散信号中所包含的最大频率值为多少？
-
- 在对一个16kHz的离散信号进行下采样至8kHz的过程中，请说明为什么要先执行低通滤波操作。
-
- 在时域上进行采样的过程（即离散化），为什么会使得频域上出现周期性现象？
-
- 时域上存在的周期性现象为何会导致频域上的离散化现象？
- 9. 实战
参考资料
1. 采用标准的WAV头部，用代码实现PCM WAV文件的读写。
在 Python 标准库 内部，“wave模块”提供了音频 WAV 格式的便捷接口。该模块的功能集能够将原始音频数据写入到诸如对象之类的目标文件，并提取WAV文件的相关属性。 wave 模块 包含 处理 WAV 音频格式 的 便利接口。该模块 不支持 压缩/解压功能 ，但 提供 单声道或立体声的支持。 contextlib 模块 能够 提供 处理 上下文 管理器 和 with 语句 的 实用程序.
#!/usr/bin/env python
# -*- coding:utf-8 -*-
# @FileName :read_write_pcm_wav.py
# @Time :2022/2/15 19:25
# @Author :PangXZ
import wave
import contextlib
import numpy as np
"""
读取test.wav，并将其转换为PCM格式
"""
def wav2pcm(wav_file, pcm_file):
if wav_file is None or pcm_file is None:
return
with contextlib.closing(wave.open(wav_file, 'r')) as f:
wav_params = f.getparams()
## wav_params：_wave_params(nchannels=, sampwidth=, framerate=, nframes=, comptype=,compname=)
print('读取WAV文件：{}'.format(wav_file))
print('声道数：{}'.format(wav_params[0]))
print('量化位数：{}'.format(wav_params[1] * 8))
# 量化位数 1,2,3 分别代表 8 位，16 位，24 位
print('采样率：{}'.format(wav_params[2]))
print('采样总数：{}'.format(wav_params[3]))
print('音频时长：{}'.format(wav_params[3] / wav_params[2]))
print('压缩格式：{}'.format(wav_params[4]))
# 读取采样值，存储为bytes格式
wav_bytes = f.readframes(wav_params[3])
# 转换为array
wave_data = np.frombuffer(wav_bytes, dtype=np.int16)
wave_data.tofile(pcm_file)
print('convert wav to pcm finish.')
""":
读取test.pcm文件，转换为WAV文件"""
def pcm2wav(pcm_file, wav_file):
if pcm_file is None or wav_file is None:
return
f = open(pcm_file, 'rb')
pcm_bytes = f.read()
with contextlib.closing(wave.open(wav_file, 'wb')) as f:
# 加入wav信息头，保存为wav格式
f.setframerate(16000) # 16kHz（采样频率）
f.setsampwidth(2) # 量化位数分为8位，16位，24位
f.setnchannels(1) # 单声道振幅数据位n*1矩阵点，立体声为n*2矩阵点
f.writeframes(pcm_bytes)
print('convert pcm to wav finish.')
if __name__ == "__main__":
wav_name = 'test.wav'
wav2pcm_name = 'test.pcm'
wav2pcm(wav_name, wav2pcm_name)
pcm2wav_name = 'pcm2wav_test.wav'
pcm2wav(wav2pcm_name, pcm2wav_name)2. 完成FBank、MFCC和PLP三种声学特征提取的代码实现
Fbank：(Filter-Bank, 滤波器组)人耳对声音频谱的反应是非线性的性质，在声学工程领域中被广泛研究。其本质是一种模拟人耳前向处理机制的技术手段，在语音信号处理中扮演着关键角色。通过这种算法设计能够显著提升语音识别效能。具体而言，在获取语音信号Fbank特征的过程中一般包含以下几个关键步骤：首先经过预加重处理以增强高频信号能量；其次按照时间间隔截取帧数据；随后施加时域加窗函数减少边缘效应的影响；接着利用短时傅里叶变换方法获取频谱信息；随后计算出功率谱并对其取平方得到幅度平方值；之后应用Mel滤波器组对频率轴进行压缩映射；最后对上述结果取自然对数以降低频谱分布的动态范围差异。完成上述操作后通过对Fbank特征执行离散余弦变换（DCT）能够提取出有效的MFCC特征用于后续的语音分析任务。
用于语音信号处理的MFCC（Mel-frequency cepstral coefficients）是一种广泛使用的特征提取方法。Mel频率是根据人耳对声音的感受特性而提出的，并与赫兹（Hz）频率之间存在非线性对应关系。其通过分析这一关系所得到的Herz频谱特征被称为梅尔倒谱系数（MFCC）。在语音数据分析中被广泛用于特征提取以及降维处理。例如，在实际应用中考虑一帧含有512个采样点的数据时，在经过MFCC处理后通常能够提取出最重要的40个频谱特征并实现了降维效果。
PLP（Perceptual Linear Predictive, 感知线性预测系数）是一种基于Bark听觉模型的特征参数，在语音增强技术中扮演着重要角色。其基本工作原理是在保证语音质量的前提下，在频域中对语音信号进行增强与噪声抑制的同时实现去噪功能。主要包含以下过程：首先进行频谱分析以获取语音信号的频率信息；然后对各临界信道进行能量计算与分配；接着通过乘以预加重因子来增强低频成分；随后将强度信号转化为相应的声压级数值；接着对处理后的数据执行逆离散傅里叶变换（IDFT）以恢复原始时域信号；最后采用线性预测方法实现语音信号的解卷处理以恢复真实的语音信息流
Fbank、MFCC、PLP语音特征提取流程图：
#!/usr/bin/env python
# -*- coding:utf-8 -*-
# @FileName :Fbank_Processing.py
# @Time :2022/2/15 20:44
# @Author :PangXZ
import sys
import wave
import librosa
import contextlib
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from scipy.fftpack import dct
def plot_time(sig, fs):
'''
绘制时域图
:param sig: 语音信号
:param fs: 采样周期
:return:
'''
time = np.arange(0, len(sig)) * (1.0 / fs)
plt.figure(figsize=(20, 5))
plt.plot(time, sig)
plt.xlabel('Times(s)')
plt.ylabel('Amplitude(dB)') # 幅度
plt.grid()
plt.show()
def plot_freq(sig, sample_rate, nfft=512):
'''
绘制嫔御徒
:param sig:
:param sample_rate:
:param nfft:
:return:
'''
freqs = np.linspace(0, sample_rate / 2, nfft // 2 + 1)
xf = np.fft.rfft(sig, nfft) / nfft
xfp = 20 * np.log10(np.clip(np.abs(xf), 1e-20, 1e100)) # 强度
plt.figure(figsize=(20, 5))
plt.plot(freqs, xfp)
plt.xlabel('Freq(Hz)')
plt.ylabel('Amplitude(dB)') # 强度
plt.grid()
plt.show()
def plot_spectrogram(spec, name):
fig = plt.figure(figsize=(20, 5))
heatmap = plt.pcolor(spec)
fig.colorbar(mappable=heatmap)
plt.xlabel('Times(s)')
plt.ylabel('Frequency(Hz)')
# tight_layout会自动调整子图参数，使之填充整个图像区域
plt.tight_layout()
plt.title(name)
plt.show()
# 定义PLP相关的两个公式:bark_change为线性频率坐标转换为Bark坐标，equal_loudness为等响度曲线
def bark_change(x):
return 6 * np.log10(x / (1200 * np.pi) + ((x / 1200 * np.pi)) ** 2 + 1) ** 0.5
def equal_loudness(x):
return ((x ** 2 + 56.8e6) * x ** 4) / ((x ** 2 + 6.3e6) ** 2 * (x ** 2 + 3.8e8))
if __name__ == '__main__':
# 数据准备
wav_file = 'test.wav'
sample_rate, signal = wavfile.read(wav_file)
# sample_rate：是wav文件的采样率， signal是文件的内容，即语音信号
# 查看语音信号的通道数
# print(f"number of channels = {signal.shape[1]}")
# 如果wav为双通道，取第一个通道数据
# signal = signal[:, 0]
signal = signal[0: int(10 * sample_rate)] # # 保留前10s数据
plot_time(signal, sample_rate)
plot_freq(signal, sample_rate)
# 1. 预加重（Preemphasis）
pre_emphasis = 0.97
emphasized_signal = np.append(signal[0], signal[1:] - pre_emphasis * signal[:-1])
plot_time(emphasized_signal, sample_rate)
plot_freq(emphasized_signal, sample_rate)
print('语音信号的维度：{}, 采样率：{}'.format(signal.shape, sample_rate))
# 2. 分帧
frame_size, frame_stride = 0.025, 0.01 # 帧长为25ms，帧移为10ms
frame_length, frame_step = int(round(frame_size * sample_rate)), int(round(frame_stride * sample_rate))
signal_length = len(emphasized_signal)
num_frame = int(np.ceil(np.abs(signal_length - frame_length) / frame_step)) + 1
pad_signal_length = (num_frame - 1) * frame_step + frame_length
pad_zero = np.zeros(int(pad_signal_length - signal_length))
# 分帧后最后一帧点数不足，则补零
pad_signal = np.append(emphasized_signal, pad_zero)
# 每个帧的下表
indices = np.arange(0, frame_length).reshape(1, -1) + np.arange(0, num_frame * frame_step, frame_step).reshape(-1,
1)
frames = pad_signal[indices]
# frames 是二维数组，每一行是一帧，列数是每帧的采样点数，之后的短时 fft 直接在每一列上操作
print('分帧后的维度：{}'.format(frames.shape))
# 3. 加汉明窗
hamming = np.hamming(frame_length)
# hamming = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(0, frame_length) /(frame_length - 1))
# 显示加窗效果
plt.figure(figsize=(20, 5))
plt.plot(hamming)
plt.grid()
plt.xlim(0, 400)
plt.ylim(0, 1)
plt.xlabel('Samples')
plt.ylabel('Amplitude')
frames = frames * hamming
plot_time(frames, sample_rate)
plot_freq(frames.reshape(-1, ), sample_rate)
print('加窗后的维度：{}'.format(frames.reshape(-1, ).shape))
# 4. 快速傅里叶变换（FFT）
NFFT = 512
mag_frames = np.absolute(np.fft.rfft(frames, NFFT)) # 幅度谱
# 5. 功率谱(Power Spectrum)
pow_frames = ((1.0 / NFFT) * (mag_frames ** 2))
# 显示功率谱
plt.figure(figsize=(20, 5))
plt.plot(pow_frames[1])
plt.grid()
plt.show()
# Fbank
low_freq_mel = 0 # 最低mel值
high_freq_mel = 2595 * np.log10(1 + (sample_rate / 2.0) / 700) # 最高mel值，最大信号频率为fs/2
print("Mel滤波器上截止频率：{}，下截止频率：{}".format(low_freq_mel, high_freq_mel))
nfilt = 40 # 滤波器个数
mel_points = np.linspace(low_freq_mel, high_freq_mel, nfilt + 2) # nfilt个mel值均匀分布与最低和最高mel之间
# 所有的 mel 中心点，为了方便后面计算 mel 滤波器组，左右两边各补一个中心点
hz_points = 700 * (10 ** (mel_points / 2595) - 1) # mel值对应回频率点，频率间隔指数化
fbank = np.zeros((nfilt, int(NFFT / 2 + 1))) # 各个mel滤波器在能量谱对应点的取值
bin = (hz_points / (sample_rate / 2)) * (NFFT / 2) # 各个mel滤波器中心点对应FFT的区域编码，找到有值的位置
# 求mel滤波器系数
for i in range(1, nfilt + 1):
left = int(bin[i - 1]) # 左边界点
center = int(bin[i]) # 中心点
right = int(bin[i + 1]) # 右边界点
for j in range(left, center):
fbank[i - 1, j + 1] = (j + 1 - bin[i - 1]) / (bin[i] - bin[i - 1])
for j in range(center, right):
fbank[i - 1, j + 1] = (bin[i + 1] - (j + 1)) / (bin[i + 1] - bin[i])
# mel滤波
filter_banks = np.dot(pow_frames, fbank.T)
filter_banks = np.where(filter_banks == 0, np.finfo(float).eps, filter_banks)
# 取对数
filter_banks = 20 * np.log10(filter_banks) # dB
print('Shape of FBank: {}'.format(filter_banks.shape))
plot_spectrogram(filter_banks.T, 'FBank')
# MFCC
num_ceps = 12 # 保留的倒谱系数的个数
mfcc = dct(filter_banks, type=2, axis=1, norm='ortho')[:, 1:(num_ceps + 1)] # 保持在2~13
print("Shape of MFCC: {}".format(mfcc.shape))
plot_spectrogram(mfcc.T, 'MFCC')
# PLP
N = int(NFFT / 2 + 1)
df = sample_rate / N # 频率分辨率
i = np.arange(N) # 只取大于0部分的频率
freq_hz = i * df # 得到实际频率坐标：0 -> sample_rate
freq_w = 2 * np.pi * np.array(freq_hz) # 转换为角频率
freq_bark = bark_change(freq_w) # 再转换为bark频率
# 选取的临界频带数量一般要大于10，覆盖常用频率范围，这里选取了15个中心频点
point_hz = [250, 350, 450, 570, 700, 840, 1000, 1170, 1370, 1600, 1850, 2150, 2500, 2900, 3400]
point_w = 2 * np.pi * np.array(point_hz) # 转换为角频率
point_bark = bark_change(point_w) # 转换为bark频率
bank = np.zeros((15, N)) # 构造15行257列的矩阵，每一行为一个滤波器向量
filter_data = np.zeros(15) # 构造15维频带能量向量
# -------临界频带分析：构造滤波器组---------
for j in range(15):
for k in range(N):
omg = freq_bark[k] - point_bark[j]
if -1.3 < omg < -0.5:
bank[j, k] = 10 ** (2.5 * (omg + 0.5))
elif -0.5 < omg < 0.5:
bank[j, k] = 1
elif 0.5 < omg < 2.5:
bank[j, k] = 10 ** (-1.0 * (omg - 0.5))
else:
bank[j, k] = 0
bark_filter_banks = np.dot(pow_frames, bank.T)
# 等响度预加重
equal_data = equal_loudness(point_w) * bark_filter_banks
# 强度响度转换，近似模拟声音的强度与人耳感受的响度间的非线性关系
cubic_data = equal_data ** 0.33
# 做30点的ifft（逆傅里叶变换），得到30维PLP向量
plp_data = np.fft.ifft(cubic_data, 30)
# 线性预测
plp = np.zeros((plp_data.shape[0], 12))
for i in range(plp_data.shape[0]):
plp[i] = librosa.lpc(abs(plp_data[i]), order=12)[1:]
# print(plp,plp.shape)
print("shape of PLP", plp.shape)
plot_spectrogram(plp.T, 'PLP')FBank谱图：
MFCC谱图：
PLP特征谱图：
3. 针对MFCC，回答以下问题：
3.1 分析采样率、帧长、帧移与MFCC矢量个数之间的关系。
- 采样频率：定义为模拟声音波形数字化过程中每秒抽取声压级样本的数量；
- 帧长：语音信号呈现准稳态特性，在一定时间段内其特性较为稳定，在处理时将其划分为一系列连续的小时间区间（称为"帧"），通常取值范围为10毫秒至30毫秒；
- 帧移：为了避免前后相邻两"帧"之间变化过于剧烈，在计算相邻"帧"之间的差异时设置一定的"重叠部分"；具体实现方法是采用一个与"帧长"大小相近但稍小的时间间隔作为步进值；
- MFCC矢量个数：即为对每个"帧"长度范围内的语音信号进行分析处理后提取出的特征向量数量
例如：语音时长总计为W_s=33ms的语音信号被采集，在均匀抽样的基础上完成数字处理；其中信号的采样频率设定为f_s=44kHz；根据抽样定理可得总的采样点数计算得出为n_s=W_s \times f_s = 33\times44,000= 2,999,998+个；对于时频分析参数设置方面，则选择合理的频域分辨率f_d=8k并且时间分辨率t_d=5ms以满足分析需求；在频域分析过程中将信号分解成m_t=\lceil (W_s - W_{min})/W_{step} \rceil + 1个子带信号；每个子带信号经过预加重滤波器处理后进行Fourier变换得到频谱数据；随后采用Hamming加窗法对频谱数据进行加窗处理以减少能量泄漏效应；接着利用改进型Cepstral分析方法提取出MFCC特征向量数量确定为\lceil (f_max - f_min)/\Delta_f \rceil + C$个
3.2 分析FFT大小与每帧采样点数的关系
为了确保频谱分析的有效性，在进行快速傅里叶变换（FFT）时要求选取足够的点数以覆盖每帧信号的采样范围（避免采样混叠现象）。通常取值应为2的幂次方以简化算法并提升计算效率。例如每帧采集了400个采样点选择L=512时效果最佳。
3.3 分析Mel频率的计算过程
梅尔频率被视为衡量声音高低的标准单位，在人类听觉系统中具有重要应用价值。具体而言，在1千赫兹以下的声音感知呈现出线性关系，在1千赫兹以上则呈现出非线性的对数关系模式。
3.4 分析DCT变换后得到的MFCC静态特征
离散余弦变换（DCT）是傅里叶变换的一种变形体，在其结果中全部都是实数而非复数。特别地，在离散余弦变换之后，大部分自然信号（如声音和图像）的能量都会集中在低频区域。这一特性使得DCT在压缩数据时表现出色：尤其是对于语音信号，在其频谱中能量主要分布在较低频率的部分，并且由于语音信号的频谱由低频包络和高频细节组成，在实际应用中我们通常只需要关注前13至20个较大的系数就能实现有效的数据压缩效果。为了进一步提取语音特征，在Mel频谱的基础上进行倒谱分析时需要执行以下步骤：首先计算对数值；接着通过离散余弦反变换得到倒谱；其中倒谱的主要成分对应于原始频谱的包络线信息而高频成分则反映了细节信息。最后通过提取倒谱中的第2至第13个系数来获得Mel频率倒谱系数MFCC；这些MFCC值即代表了当前帧语音信号的主要特征信息。需要注意的是在实际计算过程中我们并未进行傅里叶逆变换而是直接对经过三角滤波器处理后的频谱进行DCT运算从而得到MFCC系数序列这一过程虽然简化了计算但其本质与传统方法无异即通过滤波器组获取不同频率区间内的能量分布信息进而利用较少数量的MFCC系数来表征声道特性
3.5 分析一阶和二阶动态特征的计算过程
倒谱参数MFCC主要体现语音信号的静态特征。而语音信号的动力学特性则可以通过这些静态特征的变化来表征。实验结果表明，整合动态与静态特征能够显著提升系统识别性能。对于差分参数而言，其计算公式如下：
d_t=\begin{cases} C_{t+1}-C_{t} & \text{当 t < Q 时} \\ 0 & \text{当 t ≥ Q 时} \end{cases}
其中符号说明如下：
d_t表示第t个时间点的一阶差分值；
C_t表示第t个倒谱系数；
Q为倒谱系数的数量级；
\Delta t = K, 其中K可取值为1或2代表时间步长的变化率。通过将上述计算结果代入公式即可求得二阶差分的相关参数值。
进而可知，在整个构建过程中完整结构由以下几部分组成：N维向量（包含1/3N个原始系数、1/3N个一级差分系数以及1/3N个二级差分系数）以及可选性设置下的帧能量项
4. 对比分析STFT系列的声学特征与CQCC特征在频谱分布上的区别
傅里叶变换采用等间距采样策略，在时频空间实现了等分辨率的时间-频率转换；而常数Q转换（CQT）则将整个频率空间划分为J个音程，并赋予每个音程不同的八度分布特性。这使得该方法能够在低频段提供较高的频域分辨率，在高频段则可获得较高的时域分辨率。基于此方法构建的倒谱分析框架最初主要用于音乐信号处理领域，并后被广泛应用于声纹识别、仿冒攻击检测等应用场景；相较于传统声学特征方法而言
基于STFT的声学特征系列分析涵盖了语谱图、FBank、MFCC以及PLP。这些方法均采用了短时傅里叶变换（STFT），呈现出线性分辨率规律。相比之下，CQCC方法则具有几何级数的比例。FBank和MFCC采用了基于Mel频谱的能量分布模型进行处理；而相比之下，则是采用Bark频谱的人类听觉感知模型来模拟声音特性。因此，在保留更多原始特征方面，FBank表现更为突出；在去相关性方面，则是MFCC提供了较好的解决方案；相比之下，则是PLP具备更强的抗噪声能力。如果考虑相邻帧之间存在信息交互作用的可能性，则能够整合上下文信息以增强...
假设我们对语音模拟信号以16千赫兹的采样率进行采样，则所得离散信号中的最大频率值为多少？
按照奈奎斯特理论，在采样速率超过信号最高频率的两倍时
考虑一个原始信号的采样率是16kHz的情况，在将其抽样率降至8kHz的过程中, 为什么必须在抽样过程中先实施低通滤波器?
基于奈奎斯特定理，在16kHz的采样率下采集到的离散信号可能包含频率不超过8kHz的内容；将其进一步下采样至8kHz后，在该音频中只能包含最高频段不超过4kHz的信息；因此首先要实施低通滤波措施，并对高于4kHz的部分进行消除。
7.时域上的采样（离散化），导致了频域上的周期，为什么？
基于奈奎斯特定理的推导过程表明：在时域中进行采样操作等同于将连续信号与均匀幅值的脉冲序列进行相乘。这种操作实现了对连续信号的时间离散化处理。从频域分析的角度来看，在傅里叶变换下这些均匀幅值的脉冲序列依然保持其原有的形式特性，并且其频率成分是以采样率整数倍间隔分布在整个频谱范围之内。基于时域与频域的对偶特性，在时 Domain 中的乘法操作对应于 Frequency Domain 中的操作特性即卷积运算。因此，在 Time Domain 中完成采样的过程等价于在 Frequency Domain 中对原始信号与其抽样函数进行卷积处理
From the definition of DFT (real), we understand that DFT is a method to represent a discrete-time sequence as a sum of finite cosine and sine waves. For instance, a 32-point time-domain sequence can be synthesized using 17 cosine waves and 17 sine waves. The amplitudes of these cosine and sine waves correspond to the real and imaginary parts of the frequency spectrum, respectively. To extract these specific frequency components from a discrete-time sequence, one should perform cross-correlation calculations. By computing the cross-correlation function, one can obtain the amplitude information at each specified frequency, thereby completing the process of frequency analysis. The formula for DFT is given by X(k)=\sum_{n=0}^{N-1}x(n)e^{-j2\pi kn/N} where x(n) is the time-domain sequence and X(k) represents its frequency-domain counterpart. What is the relationship between DFT and FFT?
- 离散傅里叶变换（DFT）中的实数部分：Re X[k] = \sum_{i=0}^{N-1} x[i] \cdot \cos(2\pi ki/N)
- 离散傅里叶变换（DFT）中的虚数部分：Im X[k] = \sum_{i=0}^{N-1} x[i] \cdot \sin(2\pi ki/N)
这个数学表达式即用于计算时域序列与正余弦序列之间的互相关性，并确定原始序列中是否存在特定频率的正弦和余弦波形。
被分析的离散时域序列由x[i]表示；余弦波形cos(2\pi ki/N)与正弦波形sin(2\pi ki/N)统称为基础函数。同时，在信号处理中常用欧拉公式将这些基础函数表示为复指数形式。这些基础函数具有明确的时间重复特性。通过线性组合这些基波形后的信号同样呈现明显的周期性。
_https://www.zhihu.com/question/26448935/answer/1189713815_
8.时域上的周期，导致了频域上的离散，为什么？
仅限于DFT在时域与频域均呈现离散性与周期性特征
由此可知仅限于DFT能够被计算机所处理
通过分析时频域对称性我们得以推导出
同样地在频域中表示为一个特定的位置
该信号在时域和频域中遵循一定的对称性规律（如连续傅里叶展开式、积分变换以及DTFT和DFT）：
- 时域周期信号与频率特性间的对应关系：当时间轴呈现周期性变化特征且处于连续状态时，在频率轴上表现为对应的序列呈现出均匀采样的特性（即称为傅里叶级数FS）。
- 非平稳信号的时间特性与频率特性间的对应关系：在时间轴上表现出持续不断的状态变化，在频率轴上也呈现出持续不断的变化特征（即称为傅里叶变换FT）。
- 离散时间信号的时间特性与频率特性间的对应关系：当时间变量仅取整数值的情况下，在时间轴上的表现形式与在频率轴上的表现形式均呈现一种均匀分布的特点（即称为离散时间傅里叶变换DTFT）。
- 周期性与有穷性共存的时间序列在进行分析处理后的结果：当输入的时间序列同时具备有限长度和严格重复性的特点，在经过特定的数据处理方法后会得到一种新的具有相同属性的输出结果（即称为离散傅里叶级数DFS）。
9. 实战
请计算音频信号的12维Mel频谱系数以及23维FBank特征，并依赖librosa库进行相关处理。
Librosa是一个功能极为丰富的Python工具包，在音频信号处理方面提供了广泛的支持：它支持基本的时间-频率转换方法、多种时频分析算法、特征提取接口以及声音可视化功能。
#!/usr/bin/env python
# -*- coding:utf-8 -*-
# @FileName :MFCC_test.py
# @Time :2022/2/16 12:10
# @Author :PangXZ
import librosa
import numpy as np
from scipy.fftpack import dct
# Draw spectrogram pictures
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
def plot_spectrogram(spec, note, file_name):
"""
Draw the spectrogram picture
:param spec: a feature_dim by num_frames array (data)
:param note: title of the picture
:param file_name: name of the file
:return:
"""
fig = plt.figure(figsize=(20, 5))
heatmap = plt.pcolor(spec)
fig.colorbar(mappable=heatmap)
plt.xlabel('Time(s')
plt.ylabel(note)
plt.tight_layout()
plt.savefig(file_name)
def preemphasis(signal, coeff=0.97):
"""
语音信号预加重，相当于一个一阶高通滤波器
:param signal: 语音信号
:param coeff: 预加重系数,0表示不加重，默认为0.97
:return: 返回加重后的语音信号
"""
return np.append(signal[0], signal[1:] - coeff * signal[:-1])
def enframe(signal, frame_len=400, frame_shift=160, win=np.hamming(M=400)):
"""
对语音信号进行分帧和加窗（汉明窗）
:param signal: 输入的语音信号
:param frame_len: 帧长，这里用一帧中的采样点数
:param frame_shift: 帧移
:param win: 窗型
:return: 返回加窗之后的语音信号
"""
num_samples = signal.size # 语音信号的大小
num_frames = np.floor((num_samples - frame_len) / frame_shift) + 1 # 帧数
frames = np.zeros((int(num_frames), frame_len))
for i in range(int(num_frames)):
frames[i, :] = signal[i * frame_shift: i * frame_shift + frame_len]
frames[i, :] = frames[i, :] * win
return frames
def get_spectrum(frames, fft_len=512):
"""
使用快速傅里叶变换获得语音信号的频谱
:param frames: 已分帧的语音信号，num_frames by frame_len array
:param fft_len:快速傅里叶变换的长度，默认为512
:return: 返回频谱，a num_frames by fft_len/2+1 array (real)
"""
cFFT = np.fft.fft(frames, n=fft_len)
valid_len = int(fft_len / 2) + 1
spectrum = np.abs(cFFT[:, 0:valid_len])
return spectrum
def fbank(spectrum, num_filter=23, fft_len=512, sample_size=16000):
"""
从频谱中获得mel滤波器组特征
:param spectrum: 频谱
:param num_filter: mel滤波器数目，默认为23
:return: 返回fbank特征
"""
# 1. 获取一组梅尔滤波器组
low_freq_mel = 0
high_freq_mel = 2595 * np.log10(1 + (sample_size / 2) / 700) # 转换到mel尺度
mel_points = np.linspace(low_freq_mel, high_freq_mel, num_filter + 2) # mel空间中线性取点
hz_points = 700 * (np.power(10., (mel_points / 2595)) - 1) # 转回线性谱
bin = np.floor(hz_points / (sample_size / 2) * (fft_len / 2 + 1)) # 把原本的频率对应值缩放到FFT窗长上
# 2. 用滤波器组对每一帧特征滤波，计算特征与滤波器的乘积，使用np.dot
fbank = np.zeros((num_filter, int(np.floor(fft_len / 2 + 1))))
for m in range(1, 1 + num_filter):
f_left = int(bin[m - 1]) # 左边界点
f_center = int(bin[m]) # 中心点
f_right = int(bin[m + 1]) # 右边界点
for k in range(f_left, f_center):
fbank[m - 1, k] = (k - f_left) / (f_center - f_left)
for k in range(f_center, f_right):
fbank[m - 1, k] = (f_right - k) / (f_right - f_center)
filter_banks = np.dot(spectrum, fbank.T)
# finfo函数是根据括号中的类型来获得信息，获得符合这个类型的数型
# eps是取非负的最小值,np.finfo(float).eps的值是2.220446049250313e-16，即最小的正数
filter_banks = np.where(filter_banks == 0, np.finfo(float).eps, filter_banks)
# 3. 取对数操作
filter_banks = 20 * np.log10(filter_banks)
return filter_banks
def mfcc(fbank, num_mfcc=22):
"""
基于Fbank特征获取MFCC特征
:param fbank: fbank特征
:param num_mfcc: MFCC系数数目
:return: MFCC特征
"""
# feats = np.zeros((fbank.shape[0], num_mfcc))
# 1. 从上一步获得的fbank特征，计算获取的mfcc特征，主要是一个离散余弦变换的操作
mfcc = dct(fbank, type=2, axis=1, norm='ortho')[:, 1:(num_mfcc + 1)]
return mfcc
def write_file(feats, file_name):
"""
保存特征到文件中
:param feats: 语音信号特征
:param file_name: 文件名称
"""
f = open(file_name, 'w')
(row, col) = feats.shape
for i in range(row):
f.write('[')
for j in range(col):
f.write(str(feats[i, j]) + ' ')
f.write(']\n')
f.close()
if __name__ == "__main__":
# 预加重系数 0.97
alpha = 0.97
# 帧设置
sample_size = 16000 # fs 16kHz
frame_len = 400 # 一帧有400个采样点，即16000 × 0.025 = 400， 采样率fs=16kHz， 帧长为25ms
frame_shift = 160 # 帧移10ms × 采样率16kHz = 160
fft_len = 512 # 快速傅里叶变换的采样信号长度512 = 2的9次幂
# Mel滤波器组设置
num_filter = 23 # 滤波器个数
num_mfcc = 12 # 倒谱的个数
# 读取WAV格式语音文件
wav, fs = librosa.load('test.wav', sr=None)
signal = preemphasis(wav, coeff=alpha)
frames = enframe(signal, frame_len=frame_len, frame_shift=frame_shift, win=np.hamming(M=frame_len))
spectrum = get_spectrum(frames, fft_len=fft_len)
fbank_feats = fbank(spectrum, num_filter=num_filter, fft_len=fft_len, sample_size=sample_size)
mfcc_feats = mfcc(fbank_feats, num_mfcc=num_mfcc)
plot_spectrogram(fbank_feats, 'Filter Bank', 'fbank.png')
write_file(fbank_feats, 'test.fbank')
plot_spectrogram(mfcc_feats, 'MFCC', 'mfcc.png')
write_file(mfcc_feats, 'test.mfcc')FBank频谱图如下：
MFCC特征图：
参考资料
《语音识别技术及其应用领域》
《语音识别系统的设计与实现》
https://docs.python.org/zh-cn/3/library/wave.html
https://www.jianshu.com/p/94bc38e65fff
https://www.cnblogs.com/yifanrensheng/p/13510742.html
https://www.zhihu.com/question/26448935
https://zhuanlan.zhihu.com/p/449608778
【关于信号在时域和频域中的对称性规律探讨